`leading_boundary` using `GradientGrassman` fails on multiline MPOs

[leburgel]

leading_boundary(ψ::MPSMultiline, O::MPOMultiline, alg::GradientGrassmann) currently fails for non-trivial multiline MPOs. It seems the cost function definition for the multiline stat mech case is wrong, but it should be fairly straightforward to fix.

For example, when running

using TensorKit, MPSKit, MPSKitModels

# 2D lassical Ising
ψ = InfiniteMPS(ℂ^2, ℂ^10)
O = classical_ising(ComplexF64, Trivial; beta=1.0)

# but multiline
ψmulti = MPSMultiline([ψ, ψ])
Omulti = MPOMultiline([O, O])

alg = GradientGrassmann(; tol=1e-12)
ψ, = leading_boundary(ψmulti, Omulti, alg);

it can go wrong in any of the following ways depending on initialization:

ERROR: DomainError with -5.928428864807194e-17:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt(x::Float64)
   @ Base.Math ./math.jl:677
 [3] optimize(fg::typeof(MPSKit.GrassmannMPS.fg), x::MPSKit.GrassmannMPS.ManifoldPoint{MPSMultiline{InfiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, MPSKit.PerMPOInfEnv{MPOMultiline{DenseMPO{TrivialTensorMap{ComplexSpace, 2, 2, Matrix{ComplexF64}}}}, TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, MPSMultiline{InfiniteMPS{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, KrylovKit.Arnoldi{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 2, 1, Matrix{ComplexF64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{Float64}}, TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}, 2}, Matrix{TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, alg::OptimKit.ConjugateGradient{OptimKit.HagerZhang{Rational{Int64}}, Float64, OptimKit.HagerZhangLineSearch{Rational{Int64}}}; precondition::typeof(OptimKit._precondition), finalize!::typeof(OptimKit._finalize!), retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), transport!::typeof(MPSKit.GrassmannMPS.transport!), scale!::typeof(MPSKit.GrassmannMPS.scale!), add!::typeof(MPSKit.GrassmannMPS.add!), isometrictransport::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/cg.jl:74

or get stuck indefinitely on

┌ Warning: mpo is not hermitian

or get stuck indefinitely on

┌ Warning: Linesearch bracket converged to a point without satisfying Wolfe conditions?
└ @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/linesearches.jl:189

[Jutho]

The whole Grassmann stuff needs a major overhaul.

@maartenvd , if I am correct, the motivation for the current rather complicated implementation is that, if you just use the standard gradient of the left canonical form, and then put the translation to the proper gradient (multiplying with regularised inverse of right reduced matrix), the problem was that the OptimKit methods terminate to soon, because the Euclidean norm of this left canonical gradient is not the physical norm, and the physical norm can still be much bigger.

The best solution would be to just have a better / controllable termination condition in OptimKit. A better choice would probably be inner(gradient, preconditioned gradient) < tol

[lkdvos]

See also #53

[maartenvd]

@maartenvd , if I am correct, the motivation for the current rather complicated implementation is that, if you just use the standard gradient of the left canonical form, and then put the translation to the proper gradient (multiplying with regularised inverse of right reduced matrix), the problem was that the OptimKit methods terminate to soon, because the Euclidean norm of this left canonical gradient is not the physical norm, and the physical norm can still be much bigger.

Essentially, yes. I do recall optimkit's linesearch also terminating too early when you work with the "true" gradient.

The problem is that the gradient with respect to AL is g = v*C', where v is the gradient with respect to AC. Optimkit and its linesearch only look at the norm of the gradient, which converges way faster than the norm of v.

The preconditioned gradient is defined as Pg = v C' pinv(C * C'), so in an ideal situation optimkit could be configured to look at <gradient | preconditioned>, which is approximately the norm of v.

On the other hand, if we don't want to change optimkit, I had to redefine both the derivative and the inner product so that the inner product of the new gradient with itself becomes the norm of v. I defined:

• rho = regularize(C*C')
• the gradient is v C' inverse(rho)
• inner product = trace(x' rho y)

but the inner product then often (always) involves multiplying rho with its exact inverse, so I cobled together PrecGrad, which stores (vC', vC'*inverse(rho),rho). That way I could get away with never having to multiply the inverse of rho with rho, as that is another source of inaccuracy.

There was another complication that made it so that I had to structure the code in a weird way, instead of how it used to be https://github.com/maartenvd/MPSKit.jl/blob/2a9af1f8d3b7c2887608a899a7a6e3e69ea63f37/src/algorithms/grassmann.jl ...

[Jutho]

Thanks for the clarification, this is indeed how I also remember it. But clearly, modifying OptimKit to have a custom termination criterion, with a default that is inner(gradient, preconditioned_gradient)<tol, would be a far easier solution and allow us to restore the simpler/cleaner implementation.