Issues Encountered with IDMRG2 Algorithm

[sunnowrun]

I’m currently using IDMRG2 in MPSKit, but have encountered some unexpected behavior.
Below is the relevant code snippet (I’m using the MPS2MPO function to handle Grassmann MPSs).

using TensorKit, MPSKit

function MPS2MPO(mps::MPSKit.InfiniteMPS)
	mpo = [zeros(scalartype(mps.AL[i]), codomain(mps.AL[i]) ← space(mps.AL[i],2) ⊗ space(mps.AL[i],3)') for i in eachindex(mps.AL)]
	for i in eachindex(mpo)
		for (f1,f2) in fusiontrees(mpo[i])
			f1′ = FusionTree((f1.uncoupled[1], Irrep[ℤ₂](f1.uncoupled[2].n + 1 % 2)), Irrep[ℤ₂](f1.coupled.n+1%2), f1.isdual)
			f2′ = FusionTree(f2.uncoupled[[2]], f2.uncoupled[2], f2.isdual[[2]])
			if f2.uncoupled[1].n == 0
				copyto!(mpo[i][f1,f2] , mps.AL[i][f1,f2′])
			elseif f1.uncoupled[2].n == 1
				copyto!(mpo[i][f1,f2] , mps.AL[i][f1′,f2′])
			end
		end
		for (f1,f2) in fusiontrees(mpo[i])
			if f1.uncoupled[1].n == 1 && f2.uncoupled[1].n == 1
				mpo[i][f1,f2] .*= -1
			end
		end

	end
	return MPSKit.InfiniteMPO(mpo)

end
function test(T = ComplexF64, L = 2)
	Dh = 3
	x = MPSKit.InfiniteMPS(randn, T, fill(Rep[ℤ₂](0=>1, 1=>1), L), fill(Rep[ℤ₂](0=>Dh, 1=>Dh), L))
	y = MPSKit.InfiniteMPS(randn, T, fill(Rep[ℤ₂](0=>1, 1=>1), L), fill(Rep[ℤ₂](0=>Dh, 1=>Dh), L))
	x′ = MPS2MPO(x)

	res0 = x′ * y

	Dh = 18
	z1 = MPSKit.InfiniteMPS(randn, T, fill(Rep[ℤ₂](0=>1, 1=>1), L), fill(Rep[ℤ₂](0=>Dh, 1=>Dh), L))
	res1, _ = approximate(z1, (x′, y), MPSKit.IDMRG(verbosity=1, tol=1e-12, maxiter=100))

	z2 = MPSKit.InfiniteMPS(randn, T, fill(Rep[ℤ₂](0=>1, 1=>1), L), fill(Rep[ℤ₂](0=>Dh, 1=>Dh), L))
	trunc = truncdim(2*Dh)
	res2, _ = approximate(z2, (x′, y), MPSKit.IDMRG2(verbosity=1, tol=1e-12, maxiter=100, trscheme=trunc))

	@assert norm(dot(res0, res0) - 1) < 1e-10
	@assert norm(dot(res1, res1) - 1) < 1e-10
	@assert norm(dot(res2, res2) - 1) < 1e-10

	println(dot(res0, res1))
	println(dot(res0, res2))
	return res1, res2
end
test(Float64, 2)
test(ComplexF64, 2)
test(Float64, 4)
test(ComplexF64, 4)

belwo is the output

┌ Warning: non-unique fixedpoint detected
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/fixedpoint.jl:27
┌ Warning: non-unique fixedpoint detected
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/fixedpoint.jl:27
┌ Warning: IDMRG cancel 100:	err = 1.6503699291e-06	time = 1.04 sec
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/approximate/idmrg.jl:44
0.9999999992086777 - 8.43769498715119e-15im
0.9022513012776303 + 8.326672684688674e-16im
┌ Warning: IDMRG cancel 100:	err = 4.2302533960e-07	time = 0.98 sec
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/approximate/idmrg.jl:44
0.9999999999803788 + 2.1895713171580676e-14im
0.5670804416153772 - 0.10941681962026244im
┌ Warning: non-unique fixedpoint detected
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/fixedpoint.jl:27
┌ Warning: non-unique fixedpoint detected
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/fixedpoint.jl:27
┌ Warning: IDMRG cancel 100:	err = 3.7099543139e-08	time = 0.16 sec
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/approximate/idmrg.jl:44
0.9999999999999138 - 1.4464098719947963e-16im
0.9999999999999991 + 1.0607145446866777e-16im
┌ Warning: IDMRG cancel 100:	err = 1.4863279294e-09	time = 0.21 sec
└ @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/approximate/idmrg.jl:44
0.999999999999998 + 9.922618282587337e-16im
-0.16548138653124678 + 0.7396622075311708im

IDMRG2 returns incorrect results for L = 2, and also for T = ComplexF64, L = 4.

I’m currently using MPSKit v0.13.1 and TensorKit v0.14.6.

Do you know what might be causing this? Thank you in advance!

[lkdvos]

Could you maybe elaborate a bit on what the MPO is that you are trying to construct?

I haven't looked in detail but:
The non-unique fixed point warnings seem to indicate that the state you are trying to approximate is non-injective, which tends to screw up some of the algorithms.
I'm also not sure if we actually respect the phase, ie I would expect that the absolute value of the fidelity is 1, but not necessarily the actual overlap. I'm not entirely sure if we can actually fix this though, I'd have to think that trough a bit more

[sunnowrun]

Thank you for your reply.

Yes, I used MPS to represent fermionic Grassmann tensors, where each physical and auxiliary index carries a Grassmann variable. The MPO is constructed to multiply two Grassmann MPSs, as shown below. (Here, (a) and (b) correspond to the if and elseif branches in the first for (f1,f2) in ... loop, while (c) corresponds to the second for loop, which handles the additional sign introduced by the anticommutation property of Grassmann variables.)

I used transfer_spectrum to compute the spectrum of res0 (the exact result), and the largest eigenvalue appears to be non-degenerate.

To assess the fidelity, I calculated dot(res0, res2). Does the fidelity equal norm(dot(res0, res2)) in this context? If so, the fidelity output still seems to deviate from 1.

[lkdvos]

Indeed, the fidelity would be abs(dot(res0, res2)) == norm(dot(res0, res2)), but in general these computations fail if the MPS is non-injective, since almost every part of MPSKit assumes that the transfer matrix has a unique fixed point. (hence the warnings).

Am I understanding this correctly that you are trying to write an MPO that is the projector on the MPS, ie O = |psi><psi| in a fermionic language?
Is this not simply trying to approximate the MPS itself with a different MPS of given bond dimension?
This could be achieved by approximate(psi1, psi2, ...) directly.

PS, I don't think the drawings line up with the code, keep in mind that our MPO tensors are labeled (left, bottom) <- (top, right), so f1.uncoupled[2] is what I would expect in the if-branch of your construction.

Finally, have you had a look at TensorKit's support for fermionic tensors? Vect[FermionParity] acts similar to Z2Space, but takes care of the anticommutation automatically. MPSKit supports this out of the box, which might be more convenient for you, depending on what you are trying to achieve.

[sunnowrun]

However, the exact result res0 is injective, and the non-unique fixed point warning only appears during the IDMRG calculation, not in IDMRG2.
So I’m curious why IDMRG2 fails in this case, and whether there’s a way to fix it.

No, the MPO is not a projector. My goal is to multiply two Grassmann MPS. To do this, I first convert one Grassmann MPS into an MPO, and then apply it to the other MPS (i.e., MPO * MPS).
Apologies for the inconsistent convention — perhaps the diagram needs to be flipped vertically to match the (left, bottom) ← (top, right) convention.

That sounds great — I’ll give it a try.

[sunnowrun]

Sorry, there's a mistake in the diagram — the index of the MPS in (b) should be 0

[lkdvos]

I don't really understand what you mean with multiplying two MPS.
Given that an MPS is a state (equivalent of a vector), there isn't really a multiplication defined for two states (I think?).
You could take the inner product between two states, which yields a scalar, which is dot, but is the construction with the MPO then to impose the anticommutation relations?
Alternatively you could do an outer product to end up with an MPO, but that also doesn't seem to be what you are trying to get, so I'm slightly confused (this really might just be my bad and I might be missing something!)

[sunnowrun]

I'm implementing the multiplication of two Grassmann MPS as described in this paper, particularly in Sections V.B and VI.

I reformulated their approach as an MPO acting on an MPS, so that I can make use of MPSKit for the computation.

[sunnowrun]

Recently, I have explored the use of the IDMRG2 algorithm for multiplying two infinite MPOs. To achieve this, I first converted one of the MPOs into an MPS representation, which enabled the direct application of IDMRG for the multiplication. The graphical illustration of this process is presented below, where different colors indicate the correspondence between indices:

The corresponding implementation code is as follows:

using TensorKit, MPSKit


function MPO2MPS(mpo::InfiniteMPO)
    phys = physicalspace(mpo)
    fusers = [MPSKit.fuser(scalartype(mpo), s, s') for s in phys]
    res = [@tensor res[-1,-2;-3] := mpo.O[i][-1, 1,2,-3] * fusers[i][-2,1,2] for i in eachindex(fusers)]
    return InfiniteMPS(res)
end
function MPS2MPO(mps::InfiniteMPS, physpace)
    fusers = [MPSKit.fuser(scalartype(mps), physpace, physpace') for i in eachindex(mps)]
    res = [@tensor res[-1,-2;-3,-4] := mps.AL[i][-1,1,-4] * conj(fusers[i][1,-2,-3]) for i in eachindex(fusers)]
    return InfiniteMPO(res)
end
function MPO2MPO(mpo::InfiniteMPO)
    phys = physicalspace(mpo)
    fusers = [MPSKit.fuser(scalartype(mpo), s, s') for s in phys]
    res = [@tensor res[-1,-2;-3,-4] := mpo.O[i][-1, 1,2,-4] * fusers[i][-2,1,3] * conj(fusers[i][-3,2,3]) for i in eachindex(fusers)]
    return InfiniteMPO(res)
end

function multMPO(mpo1::InfiniteMPO, mpo2::InfiniteMPO)
    physpace = physicalspace(mpo1)[1]
    res = MPO2MPO(mpo1) * MPO2MPS(mpo2)
    normalize!(res)
    return MPS2MPO(res, physpace)
end
function multMPO_IDMRG(mpo1::InfiniteMPO, mpo2::InfiniteMPO, chi::Int; alg = IDMRG(verbosity=3, maxiter=1000))
    physpace = physicalspace(mpo1)[1]
    op, psi = MPO2MPO(mpo1), MPO2MPS(mpo2)
    
    ph = physicalspace(psi)[1]
    L = length(mpo1)
    ansatz = InfiniteMPS(randn, ComplexF64, fill(ph,L), fill(ℂ^chi,L))
    res, _ = approximate(ansatz, (op, psi), alg)
    normalize!(res)

    return MPS2MPO(res, physpace)
end

function distance(mpo1::InfiniteMPO, mpo2::InfiniteMPO)
    mps1 = convert(InfiniteMPS, mpo1)
    mps2 = convert(InfiniteMPS, mpo2)
    s1 = dot(mps1, mps1)
    s2 = dot(mps2, mps2)
    s12 = dot(mps1, mps2)
    @assert norm(s1-1) < 1e-10
    @assert norm(s2-1) < 1e-10
    return s1+s2-2*real(s12)
end


function test_multMPO()
    chi1, chi2 = 3, 9
    o1 = InfiniteMPO([randn(ComplexF64, (ℂ^chi1 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^chi1)), randn(ComplexF64, (ℂ^chi1 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^chi1))])
    o2 = InfiniteMPO([randn(ComplexF64, (ℂ^chi1 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^chi1)), randn(ComplexF64, (ℂ^chi1 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^chi1))])

    res1 = o1*o2
    res2 = multMPO(o1, o2)
    res3 = multMPO_IDMRG(o1, o2, chi2)
    res4 = multMPO_IDMRG(o1, o2, chi2, alg=IDMRG2(verbosity=3, maxiter=1000, trscheme=truncdim(chi2)))
    println(distance(res1, res2), " ", distance(res1, res3), " ", distance(res1, res4))
end
test_multMPO()

The output is -4.440892098500626e-16 - 6.245004513516506e-17im 3.552713678800501e-15 + 0.0im 3.1556747937940033 - 1.8735013540549517e-16im, it seems that the IDMRG2 method also producing incorrect results.