# Constant shift in MPS shift-and-invert algo

**URL:** https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174
**Category:** DMRG and Numerical Methods
**Created:** [May 26, 2022, 10:21pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174 "2022-05-26T22:21:24Z")
**Posts on this page:** 11
**Page:** 1

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [May 26, 2022, 10:21pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/1 "2022-05-26T22:21:24Z")

</div>

Hi ITensor,

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

> **[Finding matrix product state representations of highly-excited eigenstates of...](https://arxiv.org/abs/1509.01244)**
>
> A key property of many-body localization, the localization of quantum
> particles in systems with both quenched disorder and interactions, is the area
> law entanglement of even highly excited eigenstates of many-body localized
> Hamiltonians. Matrix...

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

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [May 29, 2022, 8:20pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/2 "2022-05-29T20:20:01Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [May 31, 2022, 9:25am UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/3 "2022-05-31T09:25:51Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [June 8, 2022, 3:54am UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/4 "2022-06-08T03:54:22Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [June 13, 2022, 6:58pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/5 "2022-06-13T18:58:40Z")

</div>

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/](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/](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.

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [June 13, 2022, 7:00pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/6 "2022-06-13T19:00:44Z")

</div>

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](https://github.com/ITensor/ITensors.jl/blob/main/src/mps/projmps.jl)

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

---

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [June 16, 2022, 5:08am UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/7 "2022-06-16T05:08:07Z")

</div>

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

---

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [September 27, 2022, 3:40am UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/8 "2022-09-27T03:40:20Z")

</div>

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

> [@Issues using linsolve (cont'd discussion)](https://itensor.discourse.group/t/issues-using-linsolve-contd-discussion/410):
>
> Continuing the discussion from [Constant shift in MPS shift-and-invert algo](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/6): 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- …

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [September 30, 2022, 2:10am UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/9 "2022-09-30T02:10:23Z")

</div>

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](https://github.com/ITensor/ITensorTDVP.jl)  
(Reference for the DMRG X method: [https://arxiv.org/abs/1509.00483](https://arxiv.org/abs/1509.00483))

Would that give you what you need?

---

<div class="post-metadata">

### Author: ![ray-l](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ray-l/32/76_2.png) [@ray-l](https://itensor.discourse.group/u/ray-l)
#### Post date: [October 7, 2022, 6:32pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/10 "2022-10-07T18:32:01Z")

</div>

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

```auto
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:

```auto
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

```auto
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:

```auto
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

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [October 9, 2022, 6:59pm UTC](https://itensor.discourse.group/t/constant-shift-in-mps-shift-and-invert-algo/174/11 "2022-10-09T18:59:35Z")

</div>

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](mailto: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.
