Issues using linsolve (cont'd discussion)

DMRG and Numerical Methods
1

Continuing the discussion from Constant shift in MPS shift-and-invert algo:

I finally have time to come back to this problem. In short, I’m trying to variationally solve the following linear system:
(H - \lambda) |\psi\rangle = |\phi\rangle

Or that

o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |                                
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o-      -o--o--o--o--o--o |psi>

=
o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o |phi>

(For now let’s assume \lambda=0, since I believe you can always add back the cross terms ( \lambda \langle \psi | H | \psi\rangle , \lambda^2\langle \psi|\psi\rangle, \lambda\langle \psi|\phi\rangle ) which doesn’t change the problem)
(Or there’s a way to add constant terms?)

The way I approach this, is that I defined a new position! function:

function position!(P::AbstractProjMPO, psi::MPS, phi::MPS, pos::Int)
  L = makeL!(P, psi, phi, pos - 1)
  R = makeR!(P, psi, phi, pos + nsite(P))

  # start of the left end
  if L == nothing
    return phi[pos] * phi[pos + 1] * R

  # start of the right end
  elseif R == nothing
    return L * phi[pos] * phi[pos + 1]
  
  # in the middle
  else
    return L * phi[pos] * phi[pos + 1] * R

  end 
end

with appropriate modification to makeL! and makeR!

function _makeR!(P::AbstractProjMPO, psi::MPS, phi::MPS, k::Int)
"""everything were kept the same except"""
R = R * phi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))
"""instead of R = R * psi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))"""
end

Eventually I defined a hamiltonian ‘squared’:

  H2 = contract(H', H; cutoff=1e-12)
  H2 = replaceprime(H2, 2 => 1)

And I calculate the LHS and RHS and a guess

RHS = position!(PH, psi, phi, b)
position!(PH2, psi, b)
guess = psi[b] * psi[b + 1]

Put it through linsolve

          vecs, conv = linsolve(
            PH2,
            RHS,
            guess,  0, 1;
            ishermitian=ishermitian,
            tol=linsolve_tol,
            krylovdim=linsolve_krylovdim,
            maxiter=linsolve_maxiter,
          )

This error initially occurred
LoadError: DimensionMismatch("In scalar(T) or T[], ITensor T is not a scalar

After I removed ‘guess’ from the linsolve arguments, another error

LoadError: The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, this is probably due to an index mismatch.
Common reasons for this error:
(1) You are trying to multiply the ProjMPO with the 2-site wave-function at the wrong position.
(2) `orthogonalize!` was called, changing the MPS without updating the ProjMPO.```
2

Essentially, the LHS (The PH2 here) uses all the unmodified methods of ProjMPS, it’s just that I used a different Hamiltonian, which I believe is done correctly. I can run the LHS Hamiltonian through eigsolve as well.

I guess there’s something wrong in the way I modified the RHS methods, from contraction between the same state \langle \psi | H|\psi\rangle to \langle \psi | H |\phi\rangle, but it’s not immediately clear to me what was the issue.