Non-hermitian DMRG algorithm for steady state

[qkrcksqls135]

Hi. I’m trying to get a most steady state (update MPS in DMRG sweep with eigenvector corresponding the lowest “imaginary part” eigenvalue) in Non-hermitian Hamiltonian with ITensor DMRG.

I already know about ishermitian=false option to use ITensor DMRG code for non-hermitian Hamiltonian, however, I don’t know whether the algorithm choose the ground state (real value) or steady state, if the Hamiltonian is NH.

Is ITensor support to solution of steady state in Non-hermitian Hamiltonians?

Thank you!

[VinceNeede]

Hi @qkrcksqls135, welcome in the community.

The dmrg implementation uses KrylovKit.eigsolve internally, so please give a look to the documentation of that function (Eigenvalue problems · KrylovKit.jl). In particular you can select which eigenvalue to target using the which keyword (this gets passed from dmrg as eigsolve_which_eigenvalue, I just realized this is absent in the documentation, sorry for that, I’ll open an Issues on GitHub about it), you can either use one of the already defined values, or define your custom made EigSorter.

If for Steady State you mean the state with eigenvalue which real part is closer to zero, you could sort the eigenvalues based on the relative weight of the real part:

f(\lambda) = \frac {Re(\lambda)}{\lvert \lambda\rvert}

which would correspond to this:

which = EigSorter(λ->real(λ)/hypot(λ))

As already said in KrylovKit documentation, note that this algorithm is well suited to find the extrema of the spectrum, so if the Steady State is not an extrema, the algorithm may not converge, or require more operations.

[qkrcksqls135]

Thank you for kind answering!

I mean the Steady state that with eigenvalue which imaginary part is minimum, so I can use the option :SI for KrylovKit, then I need change KrylovKit.eigsolv(T = Float64 -> Complex).

Can I change the option T of Krylovkit.eigsolv directly in dmrg ?

or some other process is needed?

[VinceNeede]

Actually for how the call to eigsolve is performed, the passed value of T doesn’t matter, T is actually inferred in eigsovle by applying the MPO to the starting state, so if you have either \psi_0 or the Hamiltonian that is complex, eigsolve will recognize that it has to use complex algebra

[qkrcksqls135]

Thank you @VinceNeede for good comments.

Now, my code well operates, but the result of DMRG is weird.

the sweep result look like oscillation between some eigenstates which are near the lowest imaginary part, and also it has different with the ED results (attached picture is ED result of 10 eigenvalues and the sweep procedure).

However, the result of real-valued case (hermitian) is well-matched with ED(converge to ground state, and also eigenvalue is same with ED).

Is it intrinsic limitation of DMRG for non-hermitian cases?
or if not, would you mind to correct some parts of my codes?

This is my full code,

using ITensors, ITensorMPS


let
  len = 15
  Nbosons = 3
  U = 2-0.5im
  J = 1.

  maxdim = [10, 10, 20, 20, 30, 30, 40, 40, 80, 80]
  nsweeps = 300
  cutoff = 1.e-13

  function constructBoseHubbardHamiltonian(J::Float64, U::Complex, length::Int, pbc::Bool)::OpSum
    os = OpSum()
    for j = 1:length-1
      os += -J, "adag", j, "a", j + 1 
      os += -J, "adag", j + 1, "a", j
    end
    for j = 1:length
      os += U / 2, "n", j, "n", j
      os += -U / 2, "n", j
    end

    if pbc
      os += -J, "adag", length, "a", 1
      os += -J, "adag", 1, "a", length
    end
    return os
  end

  function constructBoseHubbardHamiltonian(J::Float64, U::Complex, length::Int, pbc::Bool, sites::Vector{<:Index})::MPO
    os = constructBoseHubbardHamiltonian(J, U, length, pbc)
    return MPO(os, sites)
  end

  function constructInitialState(length::Int, Nbosons::Int, sites::Vector{<:Index})
    state = [n <= Nbosons ? "1" : "0" for n = 1:length]
    psi0 = random_mps(sites, state)
    return psi0
  end

  sites = siteinds("Boson", len, conserve_qns=true, dim=4)

  psi0 = constructInitialState(len, Nbosons, sites)
  
  H = constructBoseHubbardHamiltonian(J, U, len, false, sites)

  energy, psi = dmrg(H, psi0; nsweeps=nsweeps, maxdim=maxdim, cutoff=cutoff, ishermitian=false, eigsolve_which_eigenvalue = :SI)
  
end

Thank you!

[VinceNeede]

Have you tried to play a little bit with the parameters? Especially eigsolve_krylovdim, eigsolve_maxiter and eigsolve_tol, also consider doing subsequent calls to dmrg changing these values each time and requiring more precision at each step (use the state obtained at the previous call as starting state, maybe adding a random part to prevent the algorithm to get stuck). This should allow you to reach a convergence

[ryanlevy]

As Vince said, check out the list of things here: DMRG FAQs · ITensors.jl
including adding a noise term.

Increasing the eigsolve_krylovdim seems to help, I get ~6.650425306580849 - 1.364722924618177im quickly now

[miles]

Also the way our DMRG is programmed (and probably all DMRG codes out there) makes it perform better if the Hamiltonian is close to being Hermitian. If it is very strongly non-Hermitian (eg if left and right eigenvectors are too different) then it can be difficult to converge and may require an unnecessarily large bond dimension.

But it sounds like using Ryan’s approach it can work in your case.

We are also working on some new algorithms that will include an experimental non-Hermitian DMRG that computes the dominant left and right eigenvectors.

[qkrcksqls135]

Thank you, @VinceNeede , @ryanlevy , and @miles .

yours comments are very helpful for me.

I get perfectly right results now.