Why is a ProjMPO callable but an ITensor not? (Julia)

Question

I am working on certain modifications to ITensors.jl for my own applications, and one thing I want to play with is how to do an optimization step with sites from both ends of the chain (e.g. optimizing sites N and 1 while everything else stays fixed).

To this end, the ProjMPO structure seemed too rigid, so I contracted the other tensors into an ITensor object in the most standard way, as follows:

M = copy(H[1])
for j=2:N-1
    M *= H[j] * psi[j] * dag(prime(psi[j]))
end
M *= H[N]

So far, so good. However, when I then tried to call eigsolve afterward, I got the following error: "ERROR: LoadError: MethodError: objects of type ITensor{8} are not callable".

I looked into the source code to better understand what was going on, but I didn't really get anywhere. To this end, I want to ask the following: what in the code makes a ProjMPS "callable" but an ITensor not, and why was this choice made?

Answer 1

It sounds like you are trying to pass the ITensor M as the operator directly to eigsolve. If that is the case, we have not defined the operation M(v), since it is not clear what the definition of that should be in general (it makes sense if ITensor M has pairs of primed and unprimed indices, but doesn't for other ITensors). I had been considering defining M(v) for certain kinds of ITensors, but we were being conservative about it to start out.

In the meantime, you can define your own callable ITensor wrapper, which here I call an ITensorMap:

using ITensors
using KrylovKit
i = Index(2)
M = randomITensor(i', dag(i))
v0 = randomITensor(i)
# λs, vs = eigsolve(M, v0)
# Gives:
# ERROR: MethodError: objects of type ITensor{2} are not callable
struct ITensorMap
  A::ITensor
end
(M::ITensorMap)(v::ITensor) = noprime(M.A * v)
λs, vs = eigsolve(ITensorMap(M), v0)
@show norm(noprime(M * vs[1]) - vs[1] * λs[1])

This will be a way to test if you are making your custom projected MPO correctly.

However, in general, it is not as efficient to form the operator and then contract it with the ITensor v. Instead, you should define something like the ProjMPO that contains the ITensors of the projected Hamiltonian you are interested in (in your case, H[1], H[N], and then the contraction of the bulk tensors 2...N-1), and define your own function like the product function here:

https://github.com/ITensor/ITensors.jl/blob/master/src/mps/projmpo.jl#L71

which takes the ITensor "vector" and applies the appropriate tensors. In the end the two approaches are equivalent, but the order of contraction is different, which can make a big difference in the efficiency of your algorithm.

Comments

Hi Matt, so this method of making a callable object that contains an ITensor indeed worked, so thank you! However, with the code you sent me above, I was trying it out and found something curious. In particular, you can do something similar to your last line even on completely normal DMRG. Essentially, immediately after the "eigsolve" line in the dmrg function in src/mps/dmrg.jl, you can enter the following lines:

    error = norm(noprime(PH(phi)) - energy * phi)
    println("error = $error")

This error value should be basically zero after every use of eigsolve, i.e. every optimization step within every sweep (please correct me if I am wrong about this), because the resulting phi should be an eigenstate with energy as the corresponding eigenvalue (that's what eigensolve does). However, that is not what I see.  I was doing this using the J1-J2 model (I can send you code for the Hamiltonian, if you want), and for the first few sweeps, I regularly saw nontrivial errors, including as high as 0.2 (the errors were small in later sweeps, but that was just because the state had essentially converged to the ground state by then).

Do you see this as well when you try it out? Am I doing something wrong? And if there is something wrong, how should I fix it?

---

Good question about the accuracy of the eigenvectors. That is expected, because by default we use a very small Krylov space in eigsolve, of just three. So essentially, we are projecting the entire effective Hamiltonian into a matrix of size 3x3. This would seem like a pretty extreme approximation to make, but since at the beginning of DMRG the state is generally not a very good approximation anyway, the effective Hamiltonian isn't a good approximation for the actual Hamiltonian, so there is no point in solving it very accurately. Later on, you should see that the accuracy is better, since the effective Hamiltonian becomes more accurate and the previous 2-site wavefunctions are good approximations for the eigenvector, so a small Krylov size still generally works pretty well (especially if you don't need very high accuracy in your DMRG).

You can use the keyword argument `eigsolve_krylovdim` (https://github.com/ITensor/ITensors.jl/blob/master/src/mps/dmrg.jl#L102 ) to change the maximum Krylov dimension. You can also try changing `eigsolve_maxiter` and `eigsolve_tol`, which also affect the convergence of eigsolve. See the KrylovKit documentation for more information on those parameters (https://jutho.github.io/KrylovKit.jl/stable/ ). By increasing `eigsolve_krylovdim` and `eigsolve_maxiter`, you would be able to improve the accuracy of the initial eigenvector for the initial steps, however the later steps will be slower and in the end you likely will not reach a noticeably different MPS.

---

OK, that makes things a lot clearer. I didn't realize that. Thank you so much!