numpy fill an array with arrays

Question

I want to combine an unspecified (finite) number of matrices under a Kroneckerproduct. In order to do this I want to save the matrices in an array but I don't know how to do this. At the moment I have:

for i in range(LNew-2):
 for j in range(LNew-2):
 Bulk = np.empty(shape=(LNew-1,LNew-1))
 if i == j:
 Bulk[i,j]=H2
 else:
 Bulk[i,j]=idm

Here the H2 and idm are both matrices, which I want to combine under a Kronecker product. But since Bulk is an ndarray object I suppose it wont accept arraylike objects inside it.

edit:

This is the function in which I want to use this idea. I am using it to build a Hamiltonian matrix for a quantum spin chain. So H2 is the Hamiltonian for a two particle chain, H2 is a 4x4 matrix and idm is the 2x2 identity matrix.

and now the three particle chain is np.kron(H2,idm)+np.kron(idm,H2)

and for four particles np.kron(np.kron(H2,idm),idm)+np.kron(idm,np.kron(H2,idm))+np.kron(idm,np.kron(idm,H2)) and so on.

def ExpandHN(LNew):
idm = np.identity(2)
H2 = GetH(2,'N')
HNew = H2
for i in range(LNew-2):
 for j in range(LNew-2):
 Bulk = np.empty(shape=(LNew-1,LNew-1))
 if i == j:
 Bulk[i,j]=H2
 else:
 Bulk[i,j]=idm
i = 0
for i in range(LNew-2):
 for j in range(LNew-3):
 HNew += np.kron(Bulk[i,j],Bulk[i,j+1]) #something like this
return HNew

As you can see the second set of for loops hasn't been worked out.

That being said, if someone has a totaly different but working solution I would be happy with that too.

• user707650: Do you want to store the matrics as individual objects inside the Bulk matrix, or do you want to "blow up" Bulk to accomodate for H2 and idm; if the latter, what's the shape of H2 and idm?
• magnetometer: Maybe you could tell us a bit more about your complete code. As you said that your aim is to calculate a Kronecker product, are you aware of numpy.kron ?
• magnetometer: I don't quite understand the construction rule for the Hamiltonian yet. So for the two particle chain it is H2. For the three particle chain it is H3 = np.kron(H2, idm) + np.kron(idm, H2), so I would naively assume that H4 = np.kron(H3, idm) + np.kron(idm, H3) , but that doesn't seem to be the case. But might it be possible to calculate HNew in a loop without constructing Bulk?
• AppMathPhys: @magnetometer The H2 hamiltonian is the interaction of spins of two particles in the chain. the rest of the chain is kept constant(so up or down in the case of spin 1/2). So I check for each set of two particles what happens. Using np.kron(H3,idm)+np.kron(idm,H3) counts the np.kron(idm,np.kron(H2, idm)) part twice.

Answer 1

If I understand correctly, the your question boils down to how to create arrays of arrays with numpy. I would suggest to use the standard python object dict:

Bulk = dict()
for i in range(LNew-2):
 for j in range(LNew-2):
 if i == j:
 Bulk[(i,j)]=H2
 else:
 Bulk[(i,j)]=idm

The usage of tuples as keys allows you to maintain an array-like indexing of the matrices. Also note, that you should define Bulk outside of the two for loops (in any case).

HTH

• AppMathPhys: thank you very much, with this help I just reduced runtime from 45 minutes to 7.5 seconds.