From d09c02f03107b49f615b91b454e1df88b9bc188c Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 14 Nov 2024 16:53:46 +0100 Subject: [PATCH 1/3] Remove RecursiveVecs and use reallinsolve --- src/algorithms/peps_opt.jl | 2 +- src/operators/derivatives.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index b84467da..ec0b3fe9 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -345,7 +345,7 @@ function fpgrad(∂F∂x, ∂f∂x, ∂f∂A, y₀, alg::ManualIter) end function fpgrad(∂F∂x, ∂f∂x, ∂f∂A, y₀, alg::LinSolver) - y, info = linsolve(∂f∂x, ∂F∂x, y₀, alg.solver, 1, -1) + y, info = reallinsolve(∂f∂x, ∂F∂x, y₀, alg.solver, 1, -1) if alg.solver.verbosity > 0 && info.converged != 1 @warn("gradient fixed-point iteration reached maximal number of iterations:", info) end diff --git a/src/operators/derivatives.jl b/src/operators/derivatives.jl index 36c0d149..0a66cf65 100644 --- a/src/operators/derivatives.jl +++ b/src/operators/derivatives.jl @@ -48,7 +48,7 @@ end (H::PEPS_∂∂C)(x) = MPSKit.∂C(x, H.GL, H.GR) (H::PEPS_∂∂AC)(x) = MPSKit.∂AC(x, (H.top, H.bot), H.GL, H.GR) -function MPSKit.∂AC(x::RecursiveVec, O::Tuple, GL, GR) +function MPSKit.∂AC(x::Vector, O::Tuple, GL, GR) return RecursiveVec( circshift( map((v, O1, O2, l, r) -> ∂AC(v, (O1, O2), l, r), x.vecs, O[1], O[2], GL, GR), 1 @@ -203,7 +203,7 @@ end (H::PEPO_∂∂C)(x) = MPSKit.∂C(x, H.GL, H.GR) (H::PEPO_∂∂AC)(x) = MPSKit.∂AC(x, (H.top, H.bot, H.mid), H.GL, H.GR) -function MPSKit.∂AC(x::RecursiveVec, O::Tuple{T,T,P}, GL, GR) where {T,P} +function MPSKit.∂AC(x::Vector, O::Tuple{T,T,P}, GL, GR) where {T,P} return RecursiveVec( circshift( map( From 07bd8a13131bca0b4759714581ba08993555c57a Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 14 Nov 2024 17:40:10 +0100 Subject: [PATCH 2/3] Add Jacobian real linear test --- test/ctmrg/jacobian_real_linear.jl | 53 ++++++++++++++++++++++++++++++ test/runtests.jl | 3 ++ 2 files changed, 56 insertions(+) create mode 100644 test/ctmrg/jacobian_real_linear.jl diff --git a/test/ctmrg/jacobian_real_linear.jl b/test/ctmrg/jacobian_real_linear.jl new file mode 100644 index 00000000..1fb721d0 --- /dev/null +++ b/test/ctmrg/jacobian_real_linear.jl @@ -0,0 +1,53 @@ +using Test +using Random +using Zygote +using TensorKit, KrylovKit, PEPSKit +using PEPSKit: ctmrg_iter, gauge_fix, fix_relative_phases, fix_global_phases + +iterschemes = [:diffgauge, :fixed] +ctm_alg = CTMRG() +state = InfinitePEPS(2, 2) +envs = leading_boundary(CTMRGEnv(state, ComplexSpace(16)), state, ctm_alg) +@testset "$iterscheme" for iterscheme in iterschemes + # follow code of _rrule + if iterscheme == :fixed + envsconv, info = ctmrg_iter(state, envs, ctm_alg) + envsfix, signs = gauge_fix(envs, envsconv) + Ufix, Vfix = fix_relative_phases(info.U, info.V, signs) + svd_alg_fixed = SVDAdjoint(; + fwd_alg=PEPSKit.FixedSVD(Ufix, info.S, Vfix), + rrule_alg=ctm_alg.projector_alg.svd_alg.rrule_alg, + ) + alg_fixed = CTMRG(; + svd_alg=svd_alg_fixed, trscheme=notrunc(), ctmrgscheme=:simultaneous + ) + + _, env_vjp = pullback(state, envsfix) do A, x + e, = PEPSKit.ctmrg_iter(A, x, alg_fixed) + return PEPSKit.fix_global_phases(x, e) + end + elseif iterscheme == :diffgauge + _, env_vjp = pullback(state, envs) do A, x + return gauge_fix(x, ctmrg_iter(A, x, ctm_alg)[1])[1] + end + end + + # get Jacobians of single iteration + ∂f∂A(x)::typeof(state) = env_vjp(x)[1] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + + # compute real and complex errors + env_in = CTMRGEnv(state, ComplexSpace(16)) + α_real = randn(Float64) + α_complex = randn(ComplexF64) + + real_err_∂A = norm(scale(∂f∂A(env_in), α_real) - ∂f∂A(scale(env_in, α_real))) + real_err_∂x = norm(scale(∂f∂x(env_in), α_real) - ∂f∂x(scale(env_in, α_real))) + complex_err_∂A = norm(scale(∂f∂A(env_in), α_complex) - ∂f∂A(scale(env_in, α_complex))) + complex_err_∂x = norm(scale(∂f∂x(env_in), α_complex) - ∂f∂x(scale(env_in, α_complex))) + + @test real_err_∂A < 1e-10 + @test real_err_∂x < 1e-10 + @test complex_err_∂A > 1e-3 + @test complex_err_∂x > 1e-3 +end diff --git a/test/runtests.jl b/test/runtests.jl index 07a4cd2f..63e92b00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,6 +29,9 @@ end @time @safetestset "CTMRG schemes" begin include("ctmrg/ctmrgschemes.jl") end + @time @safetestset "CTMRG schemes" begin + include("ctmrg/jacobian_real_linear.jl") + end end if GROUP == "ALL" || GROUP == "BOUNDARYMPS" @time @safetestset "VUMPS" begin From 5f5182b5139980190f3e5eb3fa64d7815cf8fd6d Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 12 Dec 2024 17:14:25 +0100 Subject: [PATCH 3/3] Fix Jacobian test --- src/algorithms/peps_opt.jl | 9 +++--- test/ctmrg/jacobian_real_linear.jl | 47 ++++++++++++++++++------------ 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 6316444e..9d7b2056 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -253,13 +253,14 @@ function _rrule( ) @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg) envs = leading_boundary(envinit, state, alg) - envsconv, info = ctmrg_iteration(state, envs, alg) - envs_fixed, signs = gauge_fix(envs, envsconv) + envs_conv, info = ctmrg_iteration(state, envs, alg) + envs_fixed, signs = gauge_fix(envs, envs_conv) # Fix SVD - Ufix, Vfix = fix_relative_phases(info.U, info.V, signs) + U_fixed, V_fixed = fix_relative_phases(info.U, info.V, signs) svd_alg_fixed = SVDAdjoint(; - fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg + fwd_alg=FixedSVD(U_fixed, info.S, V_fixed), + rrule_alg=alg.projector_alg.svd_alg.rrule_alg, ) alg_fixed = @set alg.projector_alg.svd_alg = svd_alg_fixed alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc() diff --git a/test/ctmrg/jacobian_real_linear.jl b/test/ctmrg/jacobian_real_linear.jl index 1fb721d0..d91c4529 100644 --- a/test/ctmrg/jacobian_real_linear.jl +++ b/test/ctmrg/jacobian_real_linear.jl @@ -1,34 +1,45 @@ using Test using Random +using Accessors using Zygote using TensorKit, KrylovKit, PEPSKit -using PEPSKit: ctmrg_iter, gauge_fix, fix_relative_phases, fix_global_phases +using PEPSKit: ctmrg_iteration, gauge_fix, fix_relative_phases, fix_global_phases + +algs = [ + (:fixed, SimultaneousCTMRG(; projector_alg=HalfInfiniteProjector)), + (:diffgauge, SequentialCTMRG(; projector_alg=HalfInfiniteProjector)), + (:diffgauge, SimultaneousCTMRG(; projector_alg=HalfInfiniteProjector)), + # TODO: FullInfiniteProjector errors since even real_err_∂A, real_err_∂x are finite? + # (:fixed, SimultaneousCTMRG(; projector_alg=FullInfiniteProjector)), + # (:diffgauge, SequentialCTMRG(; projector_alg=FullInfiniteProjector)), + # (:diffgauge, SimultaneousCTMRG(; projector_alg=FullInfiniteProjector)), +] +Dbond, χenv = 2, 16 + +@testset "$iterscheme and $ctm_alg" for (iterscheme, ctm_alg) in algs + Random.seed!(123521938519) + state = InfinitePEPS(2, Dbond) + envs = leading_boundary(CTMRGEnv(state, ComplexSpace(χenv)), state, ctm_alg) -iterschemes = [:diffgauge, :fixed] -ctm_alg = CTMRG() -state = InfinitePEPS(2, 2) -envs = leading_boundary(CTMRGEnv(state, ComplexSpace(16)), state, ctm_alg) -@testset "$iterscheme" for iterscheme in iterschemes # follow code of _rrule if iterscheme == :fixed - envsconv, info = ctmrg_iter(state, envs, ctm_alg) - envsfix, signs = gauge_fix(envs, envsconv) - Ufix, Vfix = fix_relative_phases(info.U, info.V, signs) + envs_conv, info = ctmrg_iteration(state, envs, ctm_alg) + envs_fixed, signs = gauge_fix(envs, envs_conv) + U_fixed, V_fixed = fix_relative_phases(info.U, info.V, signs) svd_alg_fixed = SVDAdjoint(; - fwd_alg=PEPSKit.FixedSVD(Ufix, info.S, Vfix), + fwd_alg=PEPSKit.FixedSVD(U_fixed, info.S, V_fixed), rrule_alg=ctm_alg.projector_alg.svd_alg.rrule_alg, ) - alg_fixed = CTMRG(; - svd_alg=svd_alg_fixed, trscheme=notrunc(), ctmrgscheme=:simultaneous - ) + alg_fixed = @set ctm_alg.projector_alg.svd_alg = svd_alg_fixed + alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc() - _, env_vjp = pullback(state, envsfix) do A, x - e, = PEPSKit.ctmrg_iter(A, x, alg_fixed) + _, env_vjp = pullback(state, envs_fixed) do A, x + e, = PEPSKit.ctmrg_iteration(A, x, alg_fixed) return PEPSKit.fix_global_phases(x, e) end elseif iterscheme == :diffgauge _, env_vjp = pullback(state, envs) do A, x - return gauge_fix(x, ctmrg_iter(A, x, ctm_alg)[1])[1] + return gauge_fix(x, ctmrg_iteration(A, x, ctm_alg)[1])[1] end end @@ -46,8 +57,8 @@ envs = leading_boundary(CTMRGEnv(state, ComplexSpace(16)), state, ctm_alg) complex_err_∂A = norm(scale(∂f∂A(env_in), α_complex) - ∂f∂A(scale(env_in, α_complex))) complex_err_∂x = norm(scale(∂f∂x(env_in), α_complex) - ∂f∂x(scale(env_in, α_complex))) - @test real_err_∂A < 1e-10 - @test real_err_∂x < 1e-10 + @test real_err_∂A < 1e-9 + @test real_err_∂x < 1e-9 @test complex_err_∂A > 1e-3 @test complex_err_∂x > 1e-3 end