I’m using the AGP from ITensoragp, at the moment of using the linsolve function I get the following error: ERROR: MethodError: no method matching axpy!(::ComplexF64, ::MPS, ::MPS) . Previously I didn’t have this issue, so I’m not sure if there was an update or something that now doesn’t allow to use the linsolve function with MPS’s. I tried following the error and it leads to the axpy! where indeed it is defined only for type ITensor. I’m not sure what would be the best way around that. Any help you could provide would be appreciated. Thanks.
Please read: Please Read: Make It Easier to Help You for a guide on asking questions on this forum. Can you share a minimal self-contained code that reproduces the error that you are seeing as well as the versions of all of the packages and Julia version you are using?
Hi yes, will try to be more clear: I’m essentially using this:
As,ls_error = agp(HMPOsl[ii], ∂λH,As; error = false,tol=1e-6, maxiter=kwsweeps, krylovdim=kwmaxdim)
where HMPOsl[iii] is an MPO Hamiltonian, also ∂λH is an MPO, and As is some initial ansatz obtained by nested commutator. By using what is written in the ITensorAGP library, I cast this as a linear function.
X = linsolve(A, b, X₀;solver_kwargs..., kwargs...);
these are the types used:
type of A MPO
type of B MPS
type of X₀ MPS.
This eventually leads to a function within ITensors called axpy!, which is only defined for iTensors not MPS so it throws the following error.
ERROR: MethodError: no method matching axpy!(::ComplexF64, ::MPS, ::MPS)
Before I didn’t have this problem, so I’m suspecting is due to an update, but i’m not completely sure. Let me know if this is clear enough if not I could try to be more specific and thanks for the help.
Again, please read: Please Read: Make It Easier to Help You, please format your code and share a single code block that we could copy and paste and run to reproduce the error you are getting, and share all of the versions you are using (the versions of the packages you are using and the Julia version you are using). Instructions for all of that can be found in Please Read: Make It Easier to Help You.
I’m not just trying to be nitpicky or difficult, following those instructions makes it a lot easier for one of us or other users of the forum to help you. Putting in a bit more work on your end when you construct your post lowers the barrier for us to help.
To be more specific, there are a lot of undefined variables in your code, like HMPOsl, ∂λH, etc. Just sharing the types of those objects is not enough, please share a code where you show all of the packages you have loaded, the construction of all of the objects you are using which are needed to run the code, etc.
Hi, yes, I understand. The thing is I use a large number of libraries so it makes it a bit difficult to share all the necessary code. To make it simple, let’s say I literally want to run the example you have on (ITensorAGP.jl/src/agp.jl at main · ITensor/ITensorAGP.jl · GitHub) this:
using ITensors: δ, addtags, apply, replaceinds, scalartype
using KrylovKit
"""
agp(H::MPO, ∂H::MPO, X₀::MPO; use_real = false, init_cutoff = 1e-14, solver_kwargs = (;), kwargs...)
agp(H::MPO, ∂H::MPO; l = 1, use_real = false, init_cutoff = 1e-14, solver_kwargs = (;), kwargs...)
Construct the adiabatic gauge potential (AGP) as a matrix
product operator (MPO) by constructing linear equation
(H^2 \"otimes I + I \"otimes H^2 - 2 H \"otimes H)|x> = |b>
where |x> and |b> = i|∂H H - H ∂H> are matrix product states (MPSs)
in a doubly large Hilbert space compared to that of H. The AGP as an MPO
is constructed by solving for the MPS |x> (with initial ansatz |X₀>)
and converting |X> back into an MPO X.
For H with time reversal symmetry, the AGP and linear equation can be
represented with only real numbers, where solving for X_real in
(H^2 \"otimes I + I \"otimes H^2 - 2 H \"otimes H)|x_real> = |b_real>
gives X_real = iX (|b_real> = i|b> = -|∂H H - H ∂ H>).
Returns:
- `X::MPO` - the AGP as an MPO
- `ls_error::Number` - the `linsolve` error
Optional keyword arguments:
- `use_real` - boolean specifying whether to use only
real numbers/operators (when `H` has time reversal symmetry)
- `init_cutoff` - float specifying the truncation error cutoff
- `solver_kwargs` - a `NamedTuple` containing keyword arguments
that will get forwarded to the local solver, in this case
`KrylovKit.linsolve` which is a GMRES linear solver
"""
function agp(
H::MPO,
∂H::MPO,
X₀::MPO;
use_real=false,
init_cutoff=1e-14,
solver_kwargs=(;),
outputlevel=0,
kwargs...,
)
# solver kwargs
default_solver_kwargs = (; ishermitian=true, rtol=1e-4, maxiter=1, krylovdim=3)
# user input `solver_kwargs` will override `default_solver_kwargs` if duplicates exist
solver_kwargs = (; default_solver_kwargs..., solver_kwargs...)
# length of MPO
L = length(H)
# construct the sites
s = [siteinds(H)[j][2] for j in 1:L]
# construct identity MPO
I_mpo = MPO([δ(scalartype(H), dag(s[j])', s[j]) for j in 1:L])
# convert X₀ from MPO to MPS
X₀ = mpo_to_mps(X₀)
X₀ = MPS([
isodd(j) ? X₀[j] : replaceinds(X₀[j], s[Int(j / 2)]' => addtags(s[Int(j / 2)], "ket"))
for j in 1:(2L)
])
H_bra = interleave(H, addtags(I_mpo, "ket"))
H_ket = interleave(I_mpo, addtags(H, "ket"))
H1 = apply(H_bra, H_bra; cutoff=init_cutoff)
H2 = apply(H_ket, H_ket; cutoff=init_cutoff)
H3 = apply(H_bra, H_ket; cutoff=init_cutoff)
# construct A
A = +(H1, H2, -2 * H3; cutoff=init_cutoff)
# construct b
b = -(
apply(∂H, H; cutoff=init_cutoff), apply(H, ∂H; cutoff=init_cutoff); cutoff=init_cutoff
)
b = use_real ? -b : im * b
# convert b from MPO to MPS
b = mpo_to_mps(b)
b = MPS([
isodd(j) ? b[j] : replaceinds(b[j], s[Int(j / 2)]' => addtags(s[Int(j / 2)], "ket")) for
j in 1:(2L)
])
# solve linear equation Ax = b
outputlevel > 0 && @show linsolve_error(A, b, X₀)
linsolve_time = @elapsed begin
X = linsolve(A, b, X₀; solver_kwargs, kwargs...)
end
outputlevel > 0 && @show linsolve_time
linsolve_error_AXb = linsolve_error(A, b, X)
outputlevel > 0 && @show linsolve_error_AXb
# convert X from MPS to MPO
X = MPO([X[j] * X[j + 1] for j in 1:2:(2L)])
X = MPO([replaceinds(X[j], addtags(s[j], "ket") => s[j]') for j in 1:L])
return X, linsolve_error_AXb
end
hilbert = siteinds("Qubit", 4)
tmpmpo = randomMPO(hilbert)
As,ls_error = agp(tmpmpo, tmpmpo,tmpmpo)
I intentionally deleted the using using ITensorMPS, since I just cannot install it, but I guess that’s another topic. I am however using KrylovKit.linsolve. Using this I got the following problem: ERROR: MethodError: no method matching axpy!(::ComplexF64, ::MPS, ::MPS)
Closest candidates are:
axpy!(::Number, ::ITensor, ::ITensor)
This is my current version of ITensors: ITensors v0.6.11, and i’m using this version of Julia: Julia Version 1.10.0. I understand, and again thanks for all the help, if you need me to be more specific let me know.
Thanks for trying to include a minimal runnable example.
This minimal example works for me and I don’t see the error you are reporting:
using ITensorAGP: agp
using ITensorMPS: MPO, OpSum, siteinds
function main(; n, hz, nsweeps)
s = siteinds("Qubit", n)
os = OpSum()
for j in 1:(n - 1)
os += "X", j, "X", j + 1
end
for j in 1:n
os += hz, "Z", j
end
H = MPO(os, s)
dos = OpSum()
for j in 1:n
dos += "Z", j
end
dH = MPO(dos, s)
return agp(H, dH; nsweeps=2)
end
agp_hz, err = main(; n=4, hz=0.5, nsweeps=2)
Does that example work for you?
I’m using versions:
julia> using Pkg: Pkg; Pkg.status("ITensorMPS")
Status `.../Project.toml`
[0d1a4710] ITensorMPS v0.2.3
julia> versioninfo()
Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin22.4.0)
CPU: 10 × Apple M1 Max
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
and additionally the latest version of ITensorAGP.jl.
Make sure to update to the latest version of ITensorAGP.jl, say by using using Pkg: Pkg; Pkg.update(). I noticed there was an issue that some of the syntax hadn’t been updated for ITensorMPS.jl v0.2 which I just fixed in Update for ITensorMPS v0.2 by mtfishman · Pull Request #7 · ITensor/ITensorAGP.jl · GitHub.
In the future, to create a minimal runnable example, if it requires a function from a package (in this case, the agp(...) function from ITensorAGP) you can just load the package and show how you are using the function from the package, you don’t need to copy and paste the definition of the function into your minimal example. The key point is just that if we download and install the same packages with the same versions that you report and then copy the code from your post verbatim into a Julia REPL we should be able to see the same output that you report, for example the same error message.
I do understand that coming up with a minimal reproducible example can be difficult since an issue may occur inside of a larger project. However, the process of constructing a minimal reproducible example helps us a lot in reproducing your error and debugging, and additionally it can even help you in understanding the bug better and maybe even fixing it yourself. There isn’t much magic in what we do in debugging and answering questions in this forum, we have a bit more insight into what kinds of bugs might occur but in the end we go through the same process of debugging as everyone else (and my first step in debugging my own code and code posted to this forum is always to try to devise a minimal example of the bug), so any help you give us is appreciated and also can be insightful for you.
Looking back at your error message, I believe the issue you are seeing may be caused by the fact that you are not loading ITensorAGP.jl or ITensorMPS.jl.
We overload KrylovKit.linsolve for MPS/MPO objects in ITensorTDVP.jl, which is used by ITensorMPS.jl and in turn by ITensorAGP.jl. If you don’t load one of those packages than a generic version of KrylovKit.linsolve will be called, and you appear to be hitting an error in that generic version of KrylovKit.linsolve. We could make that generic version of KrylovKit.linsolve work for MPS/MPO objects but it would likely be very inefficient and not very accurate.
Instead of copying and pasting the agp function definition into your own project, you should instead install ITensorAGP.jl as instructed here: GitHub - ITensor/ITensorAGP.jl: Adiabatically evolve matrix product states (MPS) using matrix produect operator (MPO) representations of the adiabatic gauge potential (AGP). and then load the package and use the agp function it provides.
Hi Matthew, thanks for all the help. The problem is that I had conflicting versions that involved PastaQ, ITensors, ITensorAGP, and didn’t allow me to install ITensorsMPS. I manage to solve this and now everything works alright. Will try to be more concrete about my question and examples next time. Thanks!
That makes sense and explains the issue, thanks for the update.
Unfortunately right now there are compatibility issues with PastaQ.jl and some other packages that have made it difficult to update to using ITensorMPS.jl (Start using ITensorMPS.jl by mtfishman · Pull Request #309 · GTorlai/PastaQ.jl · GitHub). Out of curiousity, what features of PastaQ.jl are you using? We may be able to direct you on using ITensorMPS.jl directly, such as for gate evolution. For features like tomography, it may be best for us to split it off into a separate package, PastaQ.jl is quite large and has become difficult to maintain (we have limited manpower and our attention is more focused on rewriting NDTensors.jl and developing ITensorNetworks.jl).