Animation Nodes loops have parameters that can be reassigned from within a loop. The option to do that is accessible in the Advanced Node Settings for the Loop Input node.

With reassignable parameters it's possible to keep both a queue of matrices that you want to transform, with each single transformation generating multiple additional queue entries that have to be queued, as well as a list of all matrices that have been calculated.
A first loop is created with three parameters:

- A static list of matrices with the transformations, in this example a list of four matrices that scale an object by .5 and place it into the lower four corners of a cube
- A queue with the matrices that have to be transformed in the current iteration. Once the iteration is finished, the queue will be overwritten with the transformed matrices as a new queue.
- A list of all matrices. Once the transformed matrices are calculated for the current iteration we append them to the (initially empty) list of all matrices and reassign the parameter.
A second loop will do the transfromations for the queue items:

For each transformation matrix in the “transformations loop” we transform the whole list of queued matrices items by that transformation matrix. The Transform Matrix node is vectorized for multiple matrices but not for multiple transformations, so I chose Transformation to be the iterator and Queue Items to be a parameter and not the other way around, even though this approach would seem more intuitive to me. But in principle you can also add another loop and calculate, for each item in the queue, the list of transformed matrices.
Finally, the first loop is invoked like this:

- As the first loop doesn't have a iterator in the Loop Input node, you can set the number of iterations in the Invoke Subprogram node.
- The transformations are four matrices corresponding to the four bottom corners of a cube.
- The queue is initialized with list containing a single standard matrix of location corresponding to the origin and a scale of one.
- The initial list of all matrices is empty, so it doesn't need an input here. You could also feed in a list with the initial standard matrix above to also add the zeroth iteration to the output.
Blend file (using dev version of AN 2.1.2):

(It might be a good idea to disable automatic node execution when working with a larger number of iterations. Just as a word of caution.)