Constant shift in MPS shift-and-invert algo

DMRG and Numerical Methods
1

Hi ITensor,

I’ve been trying to write a code that’s based on this paper:

Specifically, by starting with a random MPS, |\psi\rangle , and by repeatedly solving
\frac{\partial}{\partial \phi} || O|\phi\rangle - |\psi\rangle|| = 0
I can effective obtain the state |\phi\rangle = O^{-1}|\psi\rangle (and eventually O^{-N}|\psi\rangle where N is large), where O = H - \lambda and H is the hamiltonian and \lambda is the target energy above which I want the lowest energy state.

Is there a straightforward way to do this in ITensor ? What I mean is to which extend I can use/modify existing iTensor infra to achieve this, because 1. usually a constant shift in H is trivial) and 2. it involves evaluating things like \langle \phi | O^{\dagger}O|\phi\rangle and \langle \phi |O|\psi\rangle which are not standard in DMRG.

I assume I can potentially do this by separating and adding H |\psi\rangle and \lambda |\psi\rangle, and doing the site-by-site scan completely manually. The system I have in mind is a modified extended Hubbard model

2

Interesting question thanks.

So if I understand correctly, are you wanting to solve a problem of the form

A |\phi\rangle = |\psi\rangle

for some operator A and MPS |\phi\rangle and |\psi\rangle ?

If so then yes the DMRG algorithm can be modified in a nice way to do this very thing. Mostly in the physics literature this has been done in the context of the so-called “correction vector method” but that’s a very specific example and the code one can write to solve this kind of problem can be done in a general way.

Essentially the idea is to take an existing DMRG code, which is solving an eigenvalue problem, and change the core step of the code to instead solve an A x = b type linear algebra problem.

Ok please let me know if that sounds like the kind of thing you want, or if your problem is actually something different, and we could discuss more about how to do it.

3

Hi miles,

Thank you for the response. Yes indeed, I believe what you understood is correct. The goal is to find the state |\phi\rangle, effectively finding A^{-1}|\psi\rangle.

I’ve noticed that ITensor has very efficient and easy to use DMRG code. However in this case the code is rather black-box to me and I’m a bit overwhelmed by the documentation. So it would be incredibly helpful if you could walk me through the modification process a bit. Thanks again in advance.

4

So I’ve looked through some of the source code, particularly at line 219 of the dmrg.jl. I assume this is where you would want to modify the function eigsolve? But I can’t find the source code for that function.

5

Hi, thanks for your patience while I was away on vacation.

Yes you are on the right track and that line would be a key one you’d need to change to make the dmrg function into one that solves an A |x\rangle = |b\rangle type problem instead of an eigenvalue problem.

The eigsolve function is defined in the package KrylovKit.jl. Here is the documentation for it:
https://jutho.github.io/KrylovKit.jl/stable/man/eig/

For an A |x\rangle = |b\rangle problem you’ll want the function linsolve:
https://jutho.github.io/KrylovKit.jl/stable/man/linear/

The other main ingredient you’ll need is to pass the MPS repesenting |b\rangle to your modified dmrg function and to project it into the local MPS basis that you are currently using, similar to how usual DMRG “wraps” pieces of the MPS being optimized around the Hamiltonian. In this case you’d need to use the current MPS basis to both project A and the MPS |b\rangle.

6

I believe the ProjMPS type would have most of what you need for the last part (projecting the MPS |b\rangle): https://github.com/ITensor/ITensors.jl/blob/main/src/mps/projmps.jl

You may need to modify some things though, like the product function.

7

Thanks for the reply! I will give it a try. I might have to bother you again if I’m having trouble

8

Hi miles,

You probably have seen it and working on a solution (or that my new question was formulated badly and if that’s case I do apologize). But I’ve modified the code a bit to do the iterative linear solving, which for the moment does not work.

Just in case, the new question is linked here

9

Something that may be of interest to you, and which I should have thought of sooner is that we have an implementation of the “DMRG X” algorithm available in this package:
https://github.com/ITensor/ITensorTDVP.jl
(Reference for the DMRG X method: https://arxiv.org/abs/1509.00483)

Would that give you what you need?

1 Like
10

Hi Miles,

Thanks for the information. I’m checking it at the moment.

In meantime, I worked on the previous problem a bit more and found several errors. RIght now I have a question about the ProjMPO and ProjMPS class.

The ProjMPO struct is well documented, and it states that it helps with the site-by-site scanning process

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 |psi>

Am I correct in assuming that ProjMPS does something similar:

o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o M

Where M is the ProjMPS object?

Currently, by reading the source code, eigsolve(PH::ProjMPO, phi::ITensor, ...) clearly works, but when I tried to evaluate something similar using linsolve( PH::ProjMPO, RHS::ProjMPS, phi::ITensor, ....) It would give me errors

ERROR: LoadError: MethodError: no method matching iterate(::ProjMPS)
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/range.jl:826
  iterate(::Union{LinRange, StepRangeLen}, ::Integer) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/range.jl:826
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/dict.jl:695
  ...
Stacktrace:
  [1] dot(x::ProjMPS, y::ITensor)
    @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:877
  [2] linsolve(operator::ProjMPO, b::ProjMPS, x₀::ITensor, alg::KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}, a₀::Float64, a₁::Int64)
    @ KrylovKit ~/.julia/packages/KrylovKit/kWdb6/src/linsolve/gmres.jl:4

Or without the initial guess:

ERROR: LoadError: MethodError: no method matching *(::Float64, ::ProjMPS)
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/operators.jl:655
  *(::Union{Float16, Float32, Float64}, ::BigFloat) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/mpfr.jl:406
  *(::Union{Real, Complex}, ::Union{Adjoint{var"#s859", var"#s8591"}, Transpose{var"#s859", var"#s8591"}} where {var"#s859"<:Union{Real, Complex}, var"#s8591"<:(AbstractVector)}, ::AbstractMatrix{<:Union{Real, Complex}}, ::AbstractMatrix{<:Union{Real, Complex}}) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:1172
  ...
Stacktrace:
  [1] linsolve(f::ProjMPO, b::ProjMPS, a₀::Float64, a₁::Int64; kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:ishermitian, :tol, :krylovdim, :maxiter), Tuple{Bool, Float64, Int64, Int64}}})

I’m wondering if I’m missing something here and I have to define the operator behaviors with ProjMPS

11

What is “M”? Is that your ProjMPO_MPS object?

Your drawing of a ProjMPS is correct, in terms of what it stores internally.

The downside of the ProjMPS type at the moment is that its product or contract function (these are two names for the same function) is defined to do a special kind of contracting that is only useful for the “weight penalty DMRG” method for finding excited states. It’s not the definition of contract that is needed for other algorithms like a linear solver. For this reason we are planning to generalize ProjMPS a bit in the near future but haven’t made that change yet.

If you are specifically trying to write a linear solver (i.e. Ax=b and solving for x where A is an MPO and x and b are MPS), then we are also working on a code for that at this very moment. I could share the experimental version of this code with you if you’d like. You can email me at miles@itensor.org. If you’re trying to write a different algorithm – not a linear solver – then the code we have could still help you maybe.