Linear equation solving: A*phi=psi with A hamiltonian and psi MPS gives phi as Itensor object

ITensor Julia Questions
1

Hi everyone,

I am trying to solve A*G=psi for G where A is MPO and psi is MPS but

G, _ = linsolve(
  A, B, ref_state;
  updater         = linsolve_updater,
  updater_kwargs  = (krylovdim=5, tol=1e-6, maxiter=20),
  nsweeps         = 1,
  maxdim          = [20],
  cutoff          = 1e-6,
)

this returns G as Itensor object not MPS. How do I make G MPS such that I can use it for computing inner product.

Here is a minimal full code of what I am trying to do:

using ITensors
using ITensorMPS: linsolve, linsolve_updater

# 1) Small 4-site chain
Nsites = 4
eta    = 0.1
omega  = 1.0
z      = omega + im*eta

sites = siteinds("S=1/2", Nsites; conserve_qns=true)

# 2) Build H MPO
osH = OpSum()
for j in 1:(Nsites-1)
  osH += "Sz", j, "Sz", j+1
  osH += 0.5, "S+", j, "S-", j+1
  osH += 0.5, "S-", j, "S+", j+1
end
H = MPO(osH, sites)

# 3) Ground state via DMRG
psi0 = productMPS(sites, ["Up","Dn","Up","Dn"])
_, ref_state = dmrg(H, psi0; nsweeps=2, maxdim=[10,20], cutoff=[1e-8])
println("typeof(ref_state) = ", typeof(ref_state))  # MPS

# 4) B = Sz_2 |ref_state⟩
r₂ = 2
B   = apply(op(sites, "Sz", r₂), ref_state)
orthogonalize!(B, r₂)

# 5) A = z·I − H MPO
osA = OpSum()
for j in 1:Nsites
  osA += z, "Id", j
end
for j in 1:(Nsites-1)
  osA += -1.0, "Sz", j,    "Sz", j+1
  osA += -0.5, "S+",  j,    "S-",  j+1
  osA += -0.5, "S-",  j,    "S+",  j+1
end
A = MPO(osA, sites)

# 6) Solve (z−H)·G = B
G, _ = linsolve(
  A, B, ref_state;
  updater         = linsolve_updater,
  updater_kwargs  = (krylovdim=5, tol=1e-6, maxiter=20),
  nsweeps         = 1,
  maxdim          = [20],
  cutoff          = 1e-6,
)
println("Type of G = ", typeof(G))  # ITensor
println("Type of A = ", typeof(A))  # MPO
println("Type of B = ", typeof(B))  # MPS

# 7) Attempt inner product
osSz = OpSum(); osSz += 1.0, "Sz", r₂
Msz  = MPO(osSz, sites)
inner(ref_state, Msz, G)  # will throw MethodError since G is not an MPS
2

I think that should just be:

G = linsolve(
  A, B, ref_state;
  updater         = linsolve_updater,
  updater_kwargs  = (krylovdim=5, tol=1e-6, maxiter=20),
  nsweeps         = 1,
  maxdim          = [20],
  cutoff          = 1e-6,
)

linsolve only has one object that gets output, so if you capture the output as G, _, G will be the first tensor of the MPS rather than the entire MPS.