Discrepancy between DMRG energy and <MPS|MPO|MPS>

[shreyasvm3]

Hello,

I am running one-site DMRG, and I find the following peculiarity. Here is a simple reproduction of my code:

dmrg = qtn.DMRG(H_qumpo, max_bond, bsz = 1, cutoffs = compress_cutoff, p0 = guess_mps)
dmrg.solve(tol=conv_tol, bond_dims=max_bond , max_sweeps = sweeps, sweep_sequence = 'R', verbosity = verbosity, suppress_warnings = False, cutoffs = compress_cutoff)
E = dmrg.energy
psi = dmrg.state
H_expt = (psi.H & H_qumpo.apply(psi)).contract()/(psi.H & psi).contract())

I expect E and H_expt to match. They do match for small max_bond, like 2 and 3. However, they differ significantly for max_bond >= 4. Is this expected behaviour? It is very possible that I am misunderstanding something, I am new to both DMRG and quimb.

Thank you!

[jcmgray]

Hi @shreyasvm3, I think they should match - its a bit hard to tell whats going on without the actual hamiltonian and rest of the parameters though, maybe you could post a minimal reproducer?

[shreyasvm3]

@jcmgray Thank you for your response. It is good to know that at a conceptual level the two quantities are supposed to match. It is difficult to convert my code into a small minimal reproduction -- I am buidling MPOs for Liouvillians of spin-boson systems using the manual MPO constructor. For instance, here is a short snippet of the code that constructs the MPO for my Liouvillian:

#Initialize a zero operator
Lvl = qtn.MPO_product_operator([Zero_el2, Zero_phi2])

#Coupled elec-phi
#+1j * -V1/2 * (id ⊗ |1><1|) ⊗ (id ⊗ 1-cos(phi))
Lvl += qtn.MPO_product_operator([np.kron(paulis[0], proj1), np.kron(Id_phi, cos_term)]) * +1j * -V1/2
# Add operators (sigma_k ⊗ |1><1|) ⊗ (B_k ⊗ cos_term) for k in {0,x,y,z}
for pauli, B_op in zip(paulis,B_ops):
Lvl += qtn.MPO_product_operator([np.kron(pauli, proj1), np.kron( B_op, cos_term)])

Here the objects being supplied into qtn.MPO_product_operator are numpy arrays for local operators on the MPO sites. The physical dimensions of the MPO sites are not the same across sites ( since some sites correspond to spin dofs, whereas other correspond to truncated bosonic dofs).
Given that the two quantites E and H_expt are supposed to match conceptually, I'll try to debug further. Thank you again for your comment.

I do have another question.
Let us say L is an MPO constructed as shown above, where I add together MPOs built from local site operators.
I found that L.H doesn't yield the MPO for L^{\dagger}.
I verified this as follows:
L_dense = L.to_dense()
LH_dense = L.H.to_dense()

LH_dense doesn't match L_dense.conj().T!

Do you think this could be an artifact of the way I am constucting the MPO?

[jcmgray]

Is the MPO you construct hermitian?

Regarding .H - see #315, it intentionally only conjugates the network since the 'transpose' is not well defined, but it is a bit confusing.

[shreyasvm3]

Thank you for pointing me to the thread regarding .H.

I am constructing the MPO for M=L^{\dagger}L, which is hermitian, and then perform DMRG on it to obatin the steady state (zero eigenvalue eigenvector). I do this by seperately constructing MPOs for L and L^{\dagger}, compressing them to a chosen cutoff = mpo_cutoff, and then using M_qumpo = Ld_qumpo.apply(L_qumpo, cutoff = mpo_cutoff). For very tight compression tolerances (like mpo_cutoff = 1e-30), I expect M_qumpo to be hermitian within compression errors.

I must note that since the ground state eigenvalue I am targeting is zero, the values for both E and H_expt from my post are very small (order 10^-8 to 10^-12 depending on chosen bond dimension). However, at times, E and H_expt differ by about two orders of magnitude. Is it possible this discrepancy is due to non-hermiticity due to compression?