DMRG state in computational basis when there is no transverse field

[oven8]

What is your issue?

While trying out the quimb.tensor.DMRG2, I constructed the MPO of a 4 site Hamiltonian given as,

$$H=-5\sigma^z_1-3\sigma^z_2-8\sigma^z_3-6\sigma^z_4$$

To do this I tried out an approach analogous to finite state machines (FSM). Given that this can be casted as a classical Hamiltonian (no entangling term) I would expect that the optimized DMRG state should match with brute force results using computational basis which gives me $[1,1,1,1]$ (ground state energy -22), that is all spins are up. I used the following code:

compket = qtn.MPS_computational_state('0000')
compbra = compket.copy(deep=True)
compbra.reindex_({'k0': 'b0','k1': 'b1','k2': 'b2','k3': 'b3'})
(compbra.H & H & compket)^...

where H is the MPO of the mentioned Hamiltonian. Taking the state from DRMG2 gives me $[0,0,0,0]$ as the optimized state with ground state energy -22. I would like to know if this is an issue with DMRG2 in general that it doesn't generate the correct optimized state or is this an implementation issue for my code via FSM?

More broadly speaking, is there any way in Quimb (either via DMRG or quantum circuit ansatz) through which I can get the ground state in computational basis for Hamiltonian with no transverse fields?

[jcmgray]

I think should be $|0000\rangle$ no? Since $\langle 0|\sigma_Z |0 \rangle = 1 $.

[oven8]

Sorry my bad I got confused between notations for Ising computational basis. I checked for a couple of other examples and so far, the basis state results are matching up. Thanks for the clarification. I hope it is okay if I keep this thread open for a few days while I'm testing out a few related things.

[jcmgray]

No worries, it is an easy mistake to make and partly just convention I suppose that $|\uparrow\rangle = |0\rangle$ and $|\downarrow\rangle = |1\rangle$.

[oven8]

Is there any built it method that can calculate the ground state in computational basis for states which are not entangled? So far, the only way I can think of is to calculate the density matrix via the aslinearoperator method for each pair of indices of the bra and ket state. There would be around $\mathcal{O}(2N)$ operations. Do you know if there is any better way to do this?

[jcmgray]

There's no direct function for what you want. But:

• If the groundstate is a product state, DMRG2 at least should find it and output a product state. If it is degenerate it still might be a combination/superposition however.
• You could also supply a series of bond dimensions that decrease to 1 if you wanted to force a product state output, or simply compress the final output to bond dimension 1.
• If you interested in simply sampling the highest weight bit-strings, this is quite easy to do with MPS.

[oven8]

Well, I'm dealing with states which after the DMRG optimization results in a product state. What I seek from this is to find the bit string for the optimized ket. For example, the optimized state for the previously mentioned Hamiltonian generates the bit string 0000.

Now the way I'm thinking about this is roughly on the line of

def ket_comp(i,ket):
    bra = ket.copy(deep=True)
    bra.reindex_({f'k{i}': f'b{i}'})
    BK1 = (bra.H & ket)
    BKAS = BK1.aslinearoperator([f'k{i}'],[f'b{i}'])
    if np.round(BKAS.to_dense()[0,0])==1:
        return 0
    else:
        return 1

np.vectorize(ket_comp)(np.arange(L),dmrg.state)

However, I am not sure if this is the fastest way to do it. For each index I have to separately calculate the traced density matrix. I'm new to quimb so I don't know if there exist functions that can improve the code.

[jcmgray]

Ah ok, well if its a product state then the local state is pure so you can just do, for example:

psi = dmrg.state
# remove the bonds of size 1
psi.squeeze_()
for i in range(psi.L):
    # get the local wavefunction (shorthand for `psi[psi.site_tag(i)].data`
    psi_i = psi[i].data
    # find which bit is set
    xi = np.argmax(psi_i)

[jcmgray]

For completeness, here is how you would get the density matrix for each site (if it wasn't a product state):

cur_orthog = None
for i in range(psi.L):
    # set the orthogonality center, and track it for efficiency
    psi.canonize(i, cur_orthog=cur_orthog)
    cur_orthog = i

    # form the local density matrix as a TN
    kix = psi.site_ind(i)
    bix = kix + "'"
    tn_rho_i = psi[i] & psi[i].reindex({kix: bix}).conj()

    # contract to an array
    rho_i = tn_rho_i.to_dense([kix], [bix])

something like this would be nice to include as built in functionality in quimb.

[oven8]

Ah ok, well if its a product state then the local state is pure so you can just do, for example:

psi = dmrg.state
# remove the bonds of size 1
psi.squeeze_()
for i in range(psi.L):
    # get the local wavefunction (shorthand for `psi[psi.site_tag(i)].data`
    psi_i = psi[i].data
    # find which bit is set
    xi = np.argmax(psi_i)

Some psi_i array elements seem to be -1 resulting in an incorrect bit string. Of course, something like xi = np.absolute(psi_i[1]) can fix that. But why do you think that's the case? I also saw this same thing happening while taking norm of some of the states where I was getting -1.

Thanks for this though. Your implementation is pretty fast for large $L$ states. I will check your other code, thanks for that as well!

[jcmgray]

Ah yes you will need the absolute - it's simply because eigenstates are defined up to an overall arbitrary phase: $$H e^{i\phi} |\psi\rangle = \lambda e^{i\phi} |\psi\rangle$$ which means each tensor in a tensor network wavefunction can have an arbitrary phase - taking the density matrix of course removes this, but simpler to just take the abs.