Add eigsolve-style rrule for CTMRG fixed-point gradient

[leburgel]

Solving the Sylvester problem as an eigenvalue equation, applied to the CTMRG fixed-point gradient linear problem.

Idea after:
https://github.com/Jutho/KrylovKit.jl/blob/fb56bbc7ee952bc8a6af6d278c42e78e553aa62e/ext/KrylovKitChainRulesCoreExt/eigsolve.jl#L167-L290

Actual implementation after:
https://github.com/tangwei94/VUMPSAutoDiff.jl/blob/57afc1cbcb2af9c596197c7d07adfb26ebabcd7a/src/vumps.jl#L161-L203

[codecov]

Codecov Report

Attention: Patch coverage is 86.66667% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...rithms/optimization/fixed_point_differentiation.jl	86.66%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	87.50% <ø> (ø)
...rithms/optimization/fixed_point_differentiation.jl	94.73% <86.66%> (-1.52%)	⬇️

... and 1 file with indirect coverage changes

[leburgel]

I had a go at adding this after a suggestion from @Jutho, for now mainly just so I didn't forget about it again. From what I can tell solving the fixed-point gradient linear problem like this is actually very fast and stable for all the examples. It's also the method of choice for the VUMPS pullback in VUMPSAutoDiff.jl.

However, there seems to be a problem with the gradient tests for iterscheme==:diffgauge. One possibility is that there is actually some residual gauge freedom left which makes it such that the modified eigenvalue problem doesn't work because there's no unique eigenvector with eigenvalue 1. I certainly don't really understand what is going on so I don't know an easy fix, but though I should leave this here as a suggestion.

Do you have an idea why diffgauge would give so much trouble here @pbrehmer?

[pbrehmer]

This is super interesting, thanks for adding this. I will try to take a closer look the next few days.

Regarding :diffgauge: I am kind of surprised by these issues since I thought that :diffgauge is in principle more stable. In that scheme, we differentiate through the gauge fixing itself so I would expect AD to take care of any residual gauge freedom - but perhaps I'm wrong.

[pbrehmer]

I noticed a few things but I can't claim I have really understood what's going on. The problem seems to be that :diffgauge, as you said, doesn't converge properly and outputs an eigenvalue different from 1. When looking at the output, I noticed that the :fixed mode converges after expanding the Krylov subspace to the specified dimension, whereas :diffgauge needs multiple Schur solves before it stops.

So when choosing :diffgauge, it seems that the Krylov dimension really matters. KrylovKit defaults to krylovdim=30 which might be too high in general - reducing the dimension to e.g. krylovdim=8 really speeds up the Arnoldi convergence and it will actually converge to the correct gradient with eigenvalue 1. There seems to be a certain Krylov dimension above which the :diffgauge Arnoldi will break; below that, the Arnoldi convergence is comparable to the :fixed mode.

One can also enable eager=true so that the Krylov subspace is not expanded to the full specified dimension and this also seems to repair the :diffgauge mode and even speeds it up. Perhaps this also speeds up the :fixed mode.

The difference between :fixed and :diffgauge might be a stability thing? After all, with :diffgauge each application of f(X) will differentiate through an eigenvalue problem inside gauge_fix whereas :fixed really just boils down to a fixed SVD.

[pbrehmer]

Thanks again for the addition, this really seems like an efficient approach to differentiating CTMRG! Perhaps this should be the new default? I will try to benchmark against LinSolver and see what seems best in a different PR.

[leburgel]

One can also enable eager=true so that the Krylov subspace is not expanded to the full specified dimension and this also seems to repair the :diffgauge mode and even speeds it up. Perhaps this also speeds up the :fixed mode.

Good catch, I had not thought of that! I added a default gradient_eigsolve which sets eager=true, and enabled eager mode in the tests using iterscheme=:fixed as well.

[leburgel]

Sorry to dismiss the reviews, but I was maybe too quick in un-drafting this since there were two things that still needed to be addressed before merging:

• I expanded the warning given when the norm of the auxiliary component vanishes. This basically means that either the eigsolve is unconverged, or the Jacobian doesn't have a unique leading eigenvalue 1. I was struggling with the tolerance for this check a bit, I picked the easiest thing that seemed sensible but maybe we can do better
• I switched to use realeigsolve for the same reason we switched to reallinsolve before.

If this doesn't break anything, should be good to go then