From ea0ac2023a81448ee3602beef806880e1e42a19e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 12 Oct 2026 16:45:40 +0800 Subject: [PATCH 01/14] [WIP] Split-CTMRG for iPEPS bra-ket double layer networks --- docs/src/man/environments.md | 40 ++- src/PEPSKit.jl | 10 +- .../ctmrg/enlarge_corner_split.jl | 30 +++ src/algorithms/contractions/ctmrg/gaugefix.jl | 8 +- .../contractions/ctmrg/renormalize_split.jl | 30 +++ .../contractions/ctmrg/split_env.jl | 80 ++++++ src/algorithms/ctmrg/ctmrg.jl | 18 +- src/algorithms/ctmrg/gaugefix.jl | 26 +- src/algorithms/ctmrg/gaugefix_split.jl | 80 ++++++ src/algorithms/ctmrg/projectors/projectors.jl | 2 +- src/algorithms/ctmrg/projectors/split.jl | 135 ++++++++++ src/algorithms/ctmrg/split.jl | 150 +++++++++++ .../expectation_value/expectation_value.jl | 3 + .../expectation_value/network_value.jl | 1 + .../optimization/implicit_differentiation.jl | 13 +- .../optimization/peps_optimization.jl | 2 +- .../optimization/preconditioning.jl | 4 + src/algorithms/select_algorithm.jl | 5 + src/environments/split_ctmrg_environments.jl | 250 ++++++++++++++++++ src/networks/infinitesquarenetwork.jl | 19 +- test/ctmrg/split.jl | 7 + test/gradients/split_ctmrg_gradients.jl | 6 + test/testsuite/TestSuite.jl | 13 + test/testsuite/ctmrg/split.jl | 96 +++++++ .../gradients/split_ctmrg_gradients.jl | 84 ++++++ 25 files changed, 1079 insertions(+), 33 deletions(-) create mode 100644 src/algorithms/contractions/ctmrg/enlarge_corner_split.jl create mode 100644 src/algorithms/contractions/ctmrg/renormalize_split.jl create mode 100644 src/algorithms/contractions/ctmrg/split_env.jl create mode 100644 src/algorithms/ctmrg/gaugefix_split.jl create mode 100644 src/algorithms/ctmrg/projectors/split.jl create mode 100644 src/algorithms/ctmrg/split.jl create mode 100644 src/environments/split_ctmrg_environments.jl create mode 100644 test/ctmrg/split.jl create mode 100644 test/gradients/split_ctmrg_gradients.jl create mode 100644 test/testsuite/ctmrg/split.jl create mode 100644 test/testsuite/gradients/split_ctmrg_gradients.jl diff --git a/docs/src/man/environments.md b/docs/src/man/environments.md index 485cd75e4..874b84e3f 100644 --- a/docs/src/man/environments.md +++ b/docs/src/man/environments.md @@ -1,4 +1,42 @@ # Environments !!! note - This section is still under construction. \ No newline at end of file + This section is still under construction. + +## Split CTMRG + +[`SplitCTMRGEnv`](@ref) stores separate ket and bra factors for each CTMRG edge, connected by an interlayer bond. +[`SplitCTMRG`](@ref) absorbs and projects the two layers sequentially using the full projectors of [Naumann et al., Phys. Rev. B 111, 235116 (2025)](https://doi.org/10.1103/PhysRevB.111.235116). +Its default [`SplitProjector`](@ref) constructs the four projector types directly from their split environment halves. +The outer and interlayer spaces can be chosen independently: + +```julia +env₀ = SplitCTMRGEnv(ψ, ComplexSpace(32); interlayer_space = ComplexSpace(24)) +alg = SplitCTMRG(; tol = 1.0e-8, maxiter = 100) +env, info = leading_boundary(env₀, ψ, alg) +energy = expectation_value(ψ, H, CTMRGEnv(env)) +``` + +By default, both spaces remain fixed during iteration. +Use `trunc` for the outer projectors and `interlayer_trunc` for the interlayer projectors to select other truncation strategies. +When omitted, `interlayer_trunc` inherits `trunc`; the default `FixedSpaceTruncation()` preserves the existing spaces of each bond independently. +`SplitCTMRGEnv(env::CTMRGEnv)` factorizes an existing environment exactly, while its `trunc` keyword can reduce the interlayer spaces. +`CTMRGEnv(env::SplitCTMRGEnv)` contracts the interlayer bonds so the existing expectation-value routines can be reused. +This version supports PEPS overlap networks, full projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref). +For comparisons with variPEPS, select its `Projector_Method.FULL` option; its default `FISHMAN` construction reweights the environment halves and uses different projectors. + +Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. +Both sets of bond spaces are held fixed in the differentiated iteration, including when the forward run uses other truncation strategies. +Use a complex environment for fixed-point gradients, as for conventional asymmetric CTMRG. +The first version contracts split edges into `CTMRGEnv` for energy evaluation and local preconditioning, preserving the derivative through both factors: + +```julia +env₀ = SplitCTMRGEnv(randn, ComplexF64, ψ, ComplexSpace(32); interlayer_space = ComplexSpace(24)) +opt_alg = PEPSOptimize(; + boundary_alg = SplitCTMRG(; tol = 1.0e-8), + gradient_alg = FixedPointGradient(; solver_alg = (; alg = :GMRES)), +) +ψ_opt, env_opt, energy, info = fixedpoint(H, ψ, env₀, opt_alg) +``` + +`ImplicitGradient` is not supported for split CTMRG. diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 827a16689..e4e0220c6 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -83,6 +83,7 @@ include("operators/lattices/squarelattice.jl") include("operators/models.jl") include("environments/ctmrg_environments.jl") +include("environments/split_ctmrg_environments.jl") include("environments/vumps_environments.jl") include("environments/suweight.jl") include("environments/bp_environments.jl") @@ -92,11 +93,14 @@ include("algorithms/contractions/ctmrg/types.jl") include("algorithms/contractions/ctmrg/expr.jl") include("algorithms/contractions/ctmrg/enlarge_corner.jl") include("algorithms/contractions/ctmrg/enlarge_corner_column.jl") +include("algorithms/contractions/ctmrg/enlarge_corner_split.jl") include("algorithms/contractions/ctmrg/projector.jl") include("algorithms/contractions/ctmrg/halfinf_env.jl") include("algorithms/contractions/ctmrg/fullinf_env.jl") +include("algorithms/contractions/ctmrg/split_env.jl") include("algorithms/contractions/ctmrg/renormalize_corner.jl") include("algorithms/contractions/ctmrg/renormalize_edge.jl") +include("algorithms/contractions/ctmrg/renormalize_split.jl") include("algorithms/contractions/ctmrg/network_value.jl") include("algorithms/contractions/ctmrg/gaugefix.jl") include("algorithms/contractions/ctmrg/characteristic_equations.jl") @@ -125,11 +129,14 @@ include("algorithms/ctmrg/ctmrg.jl") include("algorithms/ctmrg/projectors/projectors.jl") include("algorithms/ctmrg/projectors/halfinfinite.jl") include("algorithms/ctmrg/projectors/fullinfinite.jl") +include("algorithms/ctmrg/projectors/split.jl") include("algorithms/ctmrg/projectors/c4v_eigh.jl") include("algorithms/ctmrg/projectors/c4v_qr.jl") include("algorithms/ctmrg/simultaneous.jl") include("algorithms/ctmrg/sequential.jl") +include("algorithms/ctmrg/split.jl") include("algorithms/ctmrg/gaugefix.jl") +include("algorithms/ctmrg/gaugefix_split.jl") include("algorithms/ctmrg/c4v.jl") include("algorithms/ctmrg/initialization.jl") @@ -171,11 +178,12 @@ using .Defaults: set_scheduler! export set_scheduler! export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG +export SplitCTMRGEnv, SplitCTMRG export initialize_ctmrg_environment, RandomInitialization, ProductStateInitialization, ApplicationInitialization, IdentityInitialization export corner, edge, setcorner!, setedge! export FixedSpaceTruncation, SiteDependentTruncation -export HalfInfiniteProjector, FullInfiniteProjector +export HalfInfiniteProjector, FullInfiniteProjector, SplitProjector export C4vCTMRG, C4vEighProjector, C4vQRProjector export initialize_random_c4v_env, initialize_singlet_c4v_env export LocalOperator, physicalspace diff --git a/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl b/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl new file mode 100644 index 000000000..ba164dffd --- /dev/null +++ b/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl @@ -0,0 +1,30 @@ +# Corner enlargements +# ------------------- + +"""Project the ket cut of a northwest enlarged corner before contracting a full strip.""" +function enlarge_northwest_corner((Wk, Wb)::Tuple, C, N, A::PEPSSandwich, G) + return @autoopt @tensor Q[χout Dsb; χne Dek Deb] := + G[χout; χsw Dsk] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C[χnw; χnnw] * N[χnnw Dnk Dnb; χne] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) +end + +"""Project the ket cut of a northeast enlarged corner before contracting a full strip.""" +function enlarge_northeast_corner(N, C, (Ek, Eb)::Tuple, A::PEPSSandwich, G) + return @autoopt @tensor Q[χnw Dwk Dwb; χout Dsb] := + N[χnw Dnk Dnb; χnne] * C[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * G[χse Dsk; χout] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) +end + +"""Project the ket cut of a southwest enlarged corner before contracting a full strip.""" +function enlarge_southwest_corner(S, C, (Wk, Wb)::Tuple, A::PEPSSandwich, G) + return @autoopt @tensor Q[χse Dek Deb; χout Dnb] := + S[χse Dsk Dsb; χssw] * C[χssw; χsw] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * G[χnw Dnk; χout] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) +end + +"""Project the ket cut of a southeast enlarged corner before contracting a full strip.""" +function enlarge_southeast_corner((Ek, Eb)::Tuple, C, S, A::PEPSSandwich, G) + return @autoopt @tensor Q[χout Dnb; χsw Dwk Dwb] := + G[χout; χne Dnk] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * C[χse; χsse] * S[χsse Dsk Dsb; χsw] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) +end diff --git a/src/algorithms/contractions/ctmrg/gaugefix.jl b/src/algorithms/contractions/ctmrg/gaugefix.jl index b4d37acbb..9a069b90c 100644 --- a/src/algorithms/contractions/ctmrg/gaugefix.jl +++ b/src/algorithms/contractions/ctmrg/gaugefix.jl @@ -28,7 +28,7 @@ $(SIGNATURES) Apply `fix_gauge_corner` to the northwest corner with appropriate row and column indices. """ -function fix_gauge_northwest_corner((row, col), env::CTMRGEnv, signs) +function fix_gauge_northwest_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs) return fix_gauge_corner( corner(env, NORTHWEST, row, col), signs[WEST, row, col], @@ -41,7 +41,7 @@ $(SIGNATURES) Apply `fix_gauge_corner` to the northeast corner with appropriate row and column indices. """ -function fix_gauge_northeast_corner((row, col), env::CTMRGEnv, signs) +function fix_gauge_northeast_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs) return fix_gauge_corner( corner(env, NORTHEAST, row, col), signs[NORTH, row, col], @@ -54,7 +54,7 @@ $(SIGNATURES) Apply `fix_gauge_corner` to the southeast corner with appropriate row and column indices. """ -function fix_gauge_southeast_corner((row, col), env::CTMRGEnv, signs) +function fix_gauge_southeast_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs) return fix_gauge_corner( corner(env, SOUTHEAST, row, col), signs[EAST, row, col], @@ -67,7 +67,7 @@ $(SIGNATURES) Apply `fix_gauge_corner` to the southwest corner with appropriate row and column indices. """ -function fix_gauge_southwest_corner((row, col), env::CTMRGEnv, signs) +function fix_gauge_southwest_corner((row, col), env::Union{CTMRGEnv, SplitCTMRGEnv}, signs) return fix_gauge_corner( corner(env, SOUTHWEST, row, col), signs[SOUTH, row, col], diff --git a/src/algorithms/contractions/ctmrg/renormalize_split.jl b/src/algorithms/contractions/ctmrg/renormalize_split.jl new file mode 100644 index 000000000..e8081d842 --- /dev/null +++ b/src/algorithms/contractions/ctmrg/renormalize_split.jl @@ -0,0 +1,30 @@ +# Renormalization +# --------------- + +"""Apply the green and teal projectors sequentially when absorbing a north boundary edge.""" +function renormalize_northwest_corner(C::CTMRGCornerTensor, N::CTMRG_PEPS_EdgeTensor, green, teal) + @tensor tmp[χg Dnb; χe] := green[χg; χw Dnk] * C[χw; χn] * N[χn Dnk Dnb; χe] + return @tensor χout[χw; χe] := teal[χw; χg Dnb] * tmp[χg Dnb; χe] +end + +"""Apply the green and teal projectors sequentially when absorbing a south boundary edge.""" +function renormalize_southwest_corner(C::CTMRGCornerTensor, S::CTMRG_PEPS_EdgeTensor, green, teal) + @tensor tmp[χe Dnb; χg] := S[χe Dnk Dnb; χs] * C[χs; χw] * green[χw Dnk; χg] + return @tensor χout[χe; χw] := tmp[χe Dnb; χg] * teal[χg Dnb; χw] +end + +"""Absorb the ket layer and apply both outer and interlayer projector stages to a west edge.""" +function renormalize_west_edge(K::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor) + return @autoopt @tensor E[χinew Dek; χnnew] := + K[χiold Dwk; χnold] * A[d; Dnk Dek Dsk Dwk] * green[χnold Dnk; χng] * + yellow[χiy; χiold Dsk] * teal[χng Dnb; χnnew] * red[χinew; χiy d Dnb] +end + +"""Absorb the bra layer and apply the matching projectors to the other west-edge factor.""" +function renormalize_west_edge(B::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor, ::Val{:bra}) + @autoopt @tensor BY[χsg Dwb; χiy] := + green[χsg; χsold Dsk] * B[χsold Dwb; χiold] * yellow[χiold Dsk; χiy] + @tensor BYA[χsg Dnb Deb Dsb; χiy d] := BY[χsg Dwb; χiy] * conj(A[d; Dnb Deb Dsb Dwb]) + @tensor BYAT[χsnew Dnb Deb; χiy d] := teal[χsnew; χsg Dsb] * BYA[χsg Dnb Deb Dsb; χiy d] + return @tensor E[χsnew Deb; χinew] := BYAT[χsnew Dnb Deb; χiy d] * red[χiy d Dnb; χinew] +end diff --git a/src/algorithms/contractions/ctmrg/split_env.jl b/src/algorithms/contractions/ctmrg/split_env.jl new file mode 100644 index 000000000..c05d031f4 --- /dev/null +++ b/src/algorithms/contractions/ctmrg/split_env.jl @@ -0,0 +1,80 @@ +# Projector environments (green, teal, yellow, red) +# ------------------------------------------------- + +"""Collect the tensors surrounding a site in a split-CTMRG left move.""" +function _split_local_environment((r, c), network, env::SplitCTMRGEnv) + return (; + C1 = corner(env, NORTHWEST, r - 1, c - 1), + C2 = corner(env, NORTHEAST, r - 1, c + 1), + C3 = corner(env, SOUTHEAST, r + 1, c + 1), + C4 = corner(env, SOUTHWEST, r + 1, c - 1), + N = edge(env, NORTH, r - 1, c), S = edge(env, SOUTH, r + 1, c), + Wk = edge(env, :ket, WEST, r, c - 1), Wb = edge(env, :bra, WEST, r, c - 1), + Ek = edge(env, :ket, EAST, r, c + 1), Eb = edge(env, :bra, EAST, r, c + 1), + A = network[r, c], + ) +end + +"""Build the two halves of the environment for the green projectors (Fig. 3).""" +function _split_green_environment((r, c), network, env::SplitCTMRGEnv) + (; C1, C2, C3, C4, N, S, Wk, Wb, Ek, Eb, A) = _split_local_environment((r, c), network, env) + top = @autoopt @tensor Q[χsw Dsk; χse Dsb] := + Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χnne] * + C2[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) + bottom = @autoopt @tensor Q[χne Dnb; χnw Dnk] := + Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C4[χssw; χsw] * S[χsse Dsk Dsb; χssw] * + C3[χse; χsse] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * + ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) + return top, bottom +end + +"""Build the two halves of the environment for the teal projectors (Fig. 7).""" +function _split_teal_environment((r, c), network, env::SplitCTMRGEnv, green_west, green_east) + r = mod1(r, size(env, 2)) + l = _split_local_environment((r, c), network, env) + h = _split_local_environment((r, c + 1), network, env) + nextrow = _next(r, size(env, 2)) + NW = enlarge_northwest_corner((l.Wk, l.Wb), l.C1, l.N, l.A, green_west[2][nextrow]) + NE = enlarge_northeast_corner(h.N, h.C2, (h.Ek, h.Eb), h.A, green_east[1][nextrow]) + SW = enlarge_southwest_corner(l.S, l.C4, (l.Wk, l.Wb), l.A, green_west[1][r]) + SE = enlarge_southeast_corner((h.Ek, h.Eb), h.C3, h.S, h.A, green_east[2][r]) + return NW * NE, SE * SW +end + +"""Build the two halves of the environment for the yellow projectors (Fig. 8).""" +function _split_yellow_environment((r, c), network, env::SplitCTMRGEnv) + (; C1, C2, C3, C4, N, S, Wk, Wb, Ek, Eb, A) = _split_local_environment((r, c), network, env) + top = @autoopt @tensor Q[χiw Dsk; χie Dnb d] := + Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χnne] * C2[χnne; χne] * + Ek[χne Dek; χie] * ket(A)[d; Dnk Dek Dsk Dwk] + bottom = @autoopt @tensor Q[χie Dnb d; χiw Dsk] := + Wb[χsw Dwb; χiw] * C4[χssw; χsw] * S[χsse Dsk Dsb; χssw] * C3[χse; χsse] * + Eb[χie Deb; χse] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) + return top, bottom +end + +"""Build the two halves of the environment for the red projectors (Fig. 9).""" +function _split_red_environment((r, c), network, env::SplitCTMRGEnv, yellow_west, yellow_east) + l = _split_local_environment((r, c), network, env) + h = _split_local_environment((r, c + 1), network, env) + Ywl, Ywr = yellow_west[1][r], yellow_west[2][r] + Yel, Yer = yellow_east[1][r], yellow_east[2][r] + C1, N, Wk, A = l.C1, l.N, l.Wk, l.A + NW = @autoopt @tensor Q[χout d Dnb; χne Dek] := + Ywr[χout; χiw Dsk] * Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χne] * + ket(A)[d; Dnk Dek Dsk Dwk] + C2, N, Ek, A = h.C2, h.N, h.Ek, h.A + NE = @autoopt @tensor Q[χnw Dwk; χout d Dnb] := + N[χnw Dnk Dnb; χnne] * C2[χnne; χne] * Ek[χne Dek; χie] * Yel[χie Dsk; χout] * + ket(A)[d; Dnk Dek Dsk Dwk] + C4, S, Wb, A = l.C4, l.S, l.Wb, l.A + SW = @autoopt @tensor Q[χse Deb; χout d Dnb] := + S[χse Dsk Dsb; χssw] * C4[χssw; χsw] * Wb[χsw Dwb; χiw] * Ywl[χiw Dsk; χout] * + conj(bra(A)[d; Dnb Deb Dsb Dwb]) + C3, S, Eb, A = h.C3, h.S, h.Eb, h.A + SE = @autoopt @tensor Q[χout d Dnb; χsw Dwb] := + Yer[χout; χie Dsk] * Eb[χie Deb; χse] * C3[χse; χsse] * S[χsse Dsk Dsb; χsw] * + conj(bra(A)[d; Dnb Deb Dsb Dwb]) + return NW * NE, SE * SW +end diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index f8cc6e53f..a289d5096 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -21,7 +21,8 @@ function CTMRGAlgorithm(; verbosity = Defaults.ctmrg_verbosity, trunc = (; alg = Defaults.trunc), decomposition_alg = (;), - projector_alg = Defaults.projector_alg, # only allows for Symbol/NamedTuple to expose projector kwargs + projector_alg = alg == :SplitCTMRG ? :SplitProjector : Defaults.projector_alg, # only allows for Symbol/NamedTuple to expose projector kwargs + interlayer_trunc = trunc, ) # replace symbol with projector alg type haskey(CTMRG_SYMBOLS, alg) || throw(ArgumentError("unknown CTMRG algorithm: $alg")) @@ -39,6 +40,11 @@ function CTMRGAlgorithm(; alg = projector_alg, decomposition_alg, trunc, verbosity ) + if alg == :SplitCTMRG + projector_algorithm isa SplitProjector || throw(ArgumentError("SplitCTMRG requires SplitProjector.")) + interlayer_strategy = interlayer_trunc isa NamedTuple ? _TruncationStrategy(; interlayer_trunc...) : interlayer_trunc + return alg_type(tol, maxiter, miniter, verbosity, projector_algorithm, interlayer_strategy) + end return alg_type(tol, maxiter, miniter, verbosity, projector_algorithm) end @@ -74,6 +80,7 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st - `:SimultaneousCTMRG` : Simultaneous expansion and renormalization of all sides. - `:SequentialCTMRG` : Sequential application of left moves and rotations. - `:C4vCTMRG` : CTMRG assuming C₄ᵥ-symmetric PEPS and environment. + - `:SplitCTMRG` : Sequential ket and bra absorption using [`SplitCTMRGEnv`](@ref) and full split projectors. ### Projector algorithm @@ -87,6 +94,7 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st * `projector_alg::Symbol=:$(Defaults.projector_alg)` : Variant of the projector algorithm. See also [`ProjectorAlgorithm`](@ref). - `:HalfInfiniteProjector` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments. - `:FullInfiniteProjector` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments. + - `:SplitProjector` : Projection via SVDs of the split environments, works only for [`SplitCTMRG`](@ref) and is its default. - `:C4vEighProjector` : Projection via `eigh` of the Hermitian enlarged corner, works only for [`C4vCTMRG`](@ref). - `:C4vQRProjector` : Projection via QR decomposition of the lower-rank column-enlarged corner, works only for [`C4vCTMRG`](@ref). * `decomposition_alg::Union{NamedTuple,<:SVDAdjoint,<:EighAdjoint,<:QRAdjoint}` : Tensor @@ -103,12 +111,12 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st - `info.convergence_error::Real` : The final convergence error at the end of the contraction procedure. - `info.contraction_metrics::NamedTuple` : A `NamedTuple` containing metrics which characterize the contraction. The precise contents depend on `alg`. """ -function leading_boundary(env₀::CTMRGEnv, network::InfiniteSquareNetwork; kwargs...) +function leading_boundary(env₀::Union{CTMRGEnv, SplitCTMRGEnv}, network::InfiniteSquareNetwork; kwargs...) alg = select_algorithm(leading_boundary, env₀; kwargs...) return leading_boundary(env₀, network, alg) end function leading_boundary( - env₀::CTMRGEnv, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm + env₀::Union{CTMRGEnv, SplitCTMRGEnv}, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm ) check_input(leading_boundary, network, env₀, alg) log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) @@ -144,7 +152,7 @@ function leading_boundary( return env, info end end -function leading_boundary(env₀::CTMRGEnv, state, args...; kwargs...) +function leading_boundary(env₀::Union{CTMRGEnv, SplitCTMRGEnv}, state, args...; kwargs...) return leading_boundary(env₀, InfiniteSquareNetwork(state), args...; kwargs...) end @@ -215,7 +223,7 @@ convergence_tensors(env::CTMRGEnv, alg) = (env.corners, env.edges) Spectra of the tensors that [`convergence_tensors`](@ref) selects. """ -function convergence_spectra(env::CTMRGEnv, alg) +function convergence_spectra(env::Union{CTMRGEnv, SplitCTMRGEnv}, alg) corners, edges = convergence_tensors(env, alg) return map(C -> corner_spectrum(C, alg), corners), map(T -> edge_spectrum(T, alg), edges) end diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index 6ff44e554..6ffb7d7da 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -80,26 +80,22 @@ function compute_relative_phases( Tsfinal = circshift(reverse(envfinal.edges[dir, :, c]), r) end - # Random MPS of same bond dimension - M = map(Tsfinal) do t - randn(scalartype(t), codomain(t) ← domain(t)) - end - - # Find right fixed points of mixed transfer matrices - eigsolve_alg = Arnoldi() - ρprev = right_transfermatrix_fixedpoint(Tsprev, M, eigsolve_alg) - ρfinal = right_transfermatrix_fixedpoint(Tsfinal, M, eigsolve_alg) - - # Decompose and multiply - Qprev, = left_orth!(ρprev; positive = true) - Qfinal, = left_orth!(ρfinal; positive = true) - - return Qprev * Qfinal' + return _scrambling_gauge(Tsfinal, Tsprev) end return signs end +"""Extract a boundary-bond gauge using the same random MPS for both environments.""" +function _scrambling_gauge(final::AbstractVector, previous::AbstractVector) + M = map(t -> randn(scalartype(t), codomain(t) ← domain(t)), final) + ρprev = right_transfermatrix_fixedpoint(previous, M, Arnoldi()) + ρfinal = right_transfermatrix_fixedpoint(final, M, Arnoldi()) + Qprev, = left_orth!(ρprev; positive = true) + Qfinal, = left_orth!(ρfinal; positive = true) + return Qprev * Qfinal' +end + function compute_relative_phases(envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, ::ScramblingEnvGaugeC4v) where {C, T} # "general" algorithm from https://arxiv.org/abs/2311.11894 Tprev = envprev.edges[1, 1, 1] diff --git a/src/algorithms/ctmrg/gaugefix_split.jl b/src/algorithms/ctmrg/gaugefix_split.jl new file mode 100644 index 000000000..85792f555 --- /dev/null +++ b/src/algorithms/ctmrg/gaugefix_split.jl @@ -0,0 +1,80 @@ +"""List alternating edge factors clockwise, starting at a specified outer bond.""" +function _split_boundary_chain(env::SplitCTMRGEnv, (dir, r, c)) + dr, dc = dir == EAST ? (1, 0) : dir == WEST ? (-1, 0) : dir == NORTH ? (0, 1) : (0, -1) + count = iszero(dr) ? size(env, 3) : size(env, 2) + firstlayer, secondlayer = env.ket_first[dir] ? (:ket, :bra) : (:bra, :ket) + return [edge(env, layer, dir, r + i * dr, c + i * dc) for i in 0:(count - 1) for layer in (firstlayer, secondlayer)] +end + +"""Determine outer and interlayer basis rotations relative to the previous split environment.""" +function compute_relative_phases(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, ::ScramblingEnvGauge) + gauges = map(eachcoordinate(final, 1:4)) do coordinate + current = _split_boundary_chain(final, coordinate) + old = _split_boundary_chain(previous, coordinate) + outer = _scrambling_gauge(current, old) + interlayer = _scrambling_gauge(circshift(current, -1), circshift(old, -1)) + return outer, interlayer + end + return (; outer = map(first, gauges), interlayer = map(last, gauges)) +end + +"""Apply the outer gauges to corners and paired outer and interlayer gauges to edge factors.""" +function fix_relative_phases(env::SplitCTMRGEnv, signs) + corner_fns = (fix_gauge_northwest_corner, fix_gauge_northeast_corner, fix_gauge_southeast_corner, fix_gauge_southwest_corner) + corners = dtmap(eachcoordinate(env, 1:4)) do (dir, r, c) + return corner_fns[dir]((r, c), env, signs.outer) + end + factors = dtmap(eachcoordinate(env, 1:4)) do (dir, r, c) + dr, dc = dir == EAST ? (1, 0) : dir == WEST ? (-1, 0) : dir == NORTH ? (0, 1) : (0, -1) + next = (dir, mod1(r + dr, size(env, 2)), mod1(c + dc, size(env, 3))) + outer, interlayer = signs.outer[dir, r, c], signs.interlayer[dir, r, c] + outer_next = signs.outer[next...] + K, B = edge(env, :ket, dir, r, c), edge(env, :bra, dir, r, c) + return env.ket_first[dir] ? + (fix_gauge_edge(K, outer, interlayer), fix_gauge_edge(B, interlayer, outer_next)) : + (fix_gauge_edge(K, interlayer, outer_next), fix_gauge_edge(B, outer, interlayer)) + end + return SplitCTMRGEnv(corners, map(first, factors), map(last, factors), env.ket_first) +end + +"""Determine overlap phases for split corners and individual ket and bra edge factors.""" +function compute_global_phases(final::SplitCTMRGEnv, previous::SplitCTMRGEnv) + return _split_global_phase.(previous.corners, final.corners), + (; ket = _split_global_phase.(previous.edges_ket, final.edges_ket), bra = _split_global_phase.(previous.edges_bra, final.edges_bra)) +end + +"""Extract an overlap phase without changing the norm of a split factor.""" +function _split_global_phase(previous, current) + overlap = dot(previous, current) + return iszero(overlap) ? one(overlap) : overlap / abs(overlap) +end + +"""Remove global phases from split environment tensors without changing their norms.""" +function fix_global_phases(env::SplitCTMRGEnv, corner_phases, edge_phases) + return SplitCTMRGEnv( + env.corners .* inv.(corner_phases), + env.edges_ket .* inv.(edge_phases.ket), env.edges_bra .* inv.(edge_phases.bra), env.ket_first + ) +end + +"""Check split environment compatibility and determine relative gauges and global phases.""" +function compute_gauge_fix_gauge(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, alg::ScramblingEnvGauge) + final.ket_first == previous.ket_first || throw(ArgumentError("Split-edge orientations must match.")) + for field in (:corners, :edges_ket, :edges_bra) + space.(getproperty(final, field)) == space.(getproperty(previous, field)) || + throw(ArgumentError("Split environment spaces must match for gauge fixing.")) + end + signs = compute_relative_phases(final, previous, alg) + corner_phases, edge_phases = compute_global_phases(fix_relative_phases(final, signs), previous) + return signs, corner_phases, edge_phases +end + +"""Apply relative gauge rotations and global phase corrections to a split environment.""" +function fix_phases(env::SplitCTMRGEnv, signs, corner_phases, edge_phases) + return fix_global_phases(fix_relative_phases(env, signs), corner_phases, edge_phases) +end + +"""Align a split environment with the outer and interlayer gauges of the previous environment.""" +function gauge_fix(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, alg::ScramblingEnvGauge) + return fix_phases(final, compute_gauge_fix_gauge(final, previous, alg)...) +end diff --git a/src/algorithms/ctmrg/projectors/projectors.jl b/src/algorithms/ctmrg/projectors/projectors.jl index 6815669f8..591514556 100644 --- a/src/algorithms/ctmrg/projectors/projectors.jl +++ b/src/algorithms/ctmrg/projectors/projectors.jl @@ -25,7 +25,7 @@ function ProjectorAlgorithm(; # parse SVD forward & rrule algorithm - decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector] + decomposition_algorithm = if alg in [:HalfInfiniteProjector, :FullInfiniteProjector, :SplitProjector] _alg_or_nt(SVDAdjoint, decomposition_alg) elseif alg in [:C4vEighProjector] _alg_or_nt(EighAdjoint, decomposition_alg) diff --git a/src/algorithms/ctmrg/projectors/split.jl b/src/algorithms/ctmrg/projectors/split.jl new file mode 100644 index 000000000..201145111 --- /dev/null +++ b/src/algorithms/ctmrg/projectors/split.jl @@ -0,0 +1,135 @@ +""" +$(TYPEDEF) + +Projector algorithm for [`SplitCTMRG`](@ref), using the green, teal, yellow, and red environments of Naumann et al., Phys. Rev. B **111**, 235116 (2025). +Each projector pair is constructed from a truncated SVD of the product of its environment halves. +The yellow and red halves are QR-compressed before multiplication. + +## Fields + +$(TYPEDFIELDS) + +## Constructors + + SplitProjector(; kwargs...) + +Construct a split projector algorithm with `decomposition_alg`, `trunc`, and `verbosity`, following [`FullInfiniteProjector`](@ref). +`trunc` controls the outer cuts; [`SplitCTMRG`](@ref) supplies `interlayer_trunc` for the yellow and red cuts. +""" +struct SplitProjector{S <: SVDAdjoint, T} <: ProjectorAlgorithm + decomposition_alg::S + trunc::T + verbosity::Int +end +function SplitProjector(; kwargs...) + return ProjectorAlgorithm(; alg = :SplitProjector, kwargs...) +end + +PROJECTOR_SYMBOLS[:SplitProjector] = SplitProjector + +"""Compute a split projector pair, optionally QR-compressing the yellow or red environment halves.""" +function compute_projector( + halves, alg::SplitProjector; + trunc::TruncationStrategy = alg.trunc, compress::Bool = false, + ) + Q, Q_next = halves + # Normalize each half before multiplication. + Q = Q / norm(Q) + Q_next = Q_next / norm(Q_next) + if compress + qr_alg = QRAdjoint() + _, Q = left_orth(Q, qr_alg) + _, Q_next = left_orth(Q_next', qr_alg) + Q_next = Q_next' + end + + # All four colors use an SVD of the product of their prepared halves. + environment = Q * Q_next + svd_alg = decomposition_algorithm(_set_decomposition_truncation(alg, trunc)) + U, S, V, truncation_error = svd_trunc!(environment / norm(environment), svd_alg) + truncation_error = truncation_error / norm(S) + + Zygote.isderiving() && ignore_derivatives() do + if alg.verbosity > 0 && is_degenerate_spectrum(S) + svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S)) + @warn("degenerate singular values detected: ", svals) + end + end + + P_left, P_right = contract_projectors(U, S, V, Q, Q_next) + return (P_left, P_right), (; U, S, V, truncation_error) +end + +"""Resolve fixed-space truncation at an outer or interlayer cut.""" +function _split_truncation(trunc::TruncationStrategy, V::ElementarySpace) + return trunc isa FixedSpaceTruncation ? truncspace(isdual(V) ? flip(V) : V) : trunc +end + +"""Compute green or yellow projector pairs before constructing the teal or red environments.""" +function _split_initial_projectors( + col::Int, network, env::SplitCTMRGEnv, alg::SplitProjector, trunc::TruncationStrategy, + ::Val{interlayer}, ::Val{east} + ) where {interlayer, east} + halves = dtmap(axes(env, 2)) do r + return interlayer ? _split_yellow_environment((r, col), network, env) : + _split_green_environment((r, col), network, env) + end + results = dtmap(axes(env, 2)) do r + if interlayer + top, bottom = halves[r] + if east + top = permute(top, ((1, 4, 5), (3, 2))) + bottom = permute(bottom, ((1, 5), (4, 2, 3))) + end + Q, Q_next = east ? (top, bottom) : (bottom, top) + K = edge(env, :ket, east ? EAST : WEST, r, east ? col + 1 : col) + V = east ? domain(K)[1] : codomain(K)[1] + else + top = halves[_prev(r, size(env, 2))][1] + bottom = halves[r][2] + if east + top = permute(top, ((1, 4), (3, 2))) + bottom = permute(bottom, ((1, 4), (3, 2))) + end + Q, Q_next = east ? (top, bottom) : (bottom, top) + K = edge(env, :ket, east ? EAST : WEST, r, east ? col + 1 : col) + V = east ? codomain(K)[1] : domain(K)[1] + end + return compute_projector((Q, Q_next), alg; trunc = _split_truncation(trunc, V), compress = interlayer) + end + projectors, info = _split_proj_and_info(results) + return projectors, info.truncation_error +end + +"""Construct the green, teal, yellow, and red projector pairs for a split-CTMRG left move.""" +function sequential_projectors( + col::Int, network, env::SplitCTMRGEnv, alg::SplitProjector, interlayer_trunc::TruncationStrategy + ) + # Initial projectors on both sides of the two-site refinement environments. + green, eg = _split_initial_projectors(col, network, env, alg, alg.trunc, Val(false), Val(false)) + green_east, ege = _split_initial_projectors(col + 1, network, env, alg, alg.trunc, Val(false), Val(true)) + yellow, ey = _split_initial_projectors(col, network, env, alg, interlayer_trunc, Val(true), Val(false)) + yellow_east, eye = _split_initial_projectors(col + 1, network, env, alg, interlayer_trunc, Val(true), Val(true)) + + # Teal refines the outer cut using environments with green projectors inserted. + teal_halves = dtmap(r -> _split_teal_environment((r, col), network, env, green, green_east), axes(env, 2)) + teals = dtmap(axes(env, 2)) do r + top = teal_halves[_prev(r, size(env, 2))][1] + bottom = teal_halves[r][2] + V = domain(edge(env, :ket, WEST, r, col))[1] + trunc = _split_truncation(alg.trunc, V) + return compute_projector((bottom, top), alg; trunc) + end + + # Red refines the interlayer cut using environments with yellow projectors inserted. + reds = dtmap(axes(env, 2)) do r + top, bottom = _split_red_environment((r, col), network, env, yellow, yellow_east) + V = codomain(edge(env, :ket, WEST, r, col))[1] + trunc = _split_truncation(interlayer_trunc, V) + return compute_projector((bottom, top), alg; trunc, compress = true) + end + teal, teal_info = _split_proj_and_info(teals) + red, red_info = _split_proj_and_info(reds) + truncation_error = max(eg, ege, ey, eye, teal_info.truncation_error, red_info.truncation_error) + return (; green, teal, yellow, red), (; truncation_error) +end diff --git a/src/algorithms/ctmrg/split.jl b/src/algorithms/ctmrg/split.jl new file mode 100644 index 000000000..e2518fe77 --- /dev/null +++ b/src/algorithms/ctmrg/split.jl @@ -0,0 +1,150 @@ +""" + SplitCTMRG(; kwargs...) + +Sequential split-CTMRG introduced in Naumann et al., Phys. Rev. B **111**, 235116 (2025), Figs. 3 and 7–9. +Uses [`SplitCTMRGEnv`](@ref), with independently truncated outer and interlayer environment bonds. +The common iteration, `trunc`, and `decomposition_alg` keywords follow [`SequentialCTMRG`](@ref). +`interlayer_trunc` controls the interlayer projectors and defaults to `trunc`; specify it separately to override the outer truncation strategy. +Uses [`SplitProjector`](@ref) to construct the green, teal, yellow, and red projector pairs. +The fixed-point gradient aligns the outer and interlayer gauges and holds both environment spaces fixed. +""" +struct SplitCTMRG{P <: SplitProjector, I <: TruncationStrategy} <: CTMRGAlgorithm + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::P + interlayer_trunc::I +end +SplitCTMRG(; kwargs...) = CTMRGAlgorithm(; alg = :SplitCTMRG, kwargs...)::SplitCTMRG +CTMRG_SYMBOLS[:SplitCTMRG] = SplitCTMRG + +function check_input(::typeof(leading_boundary), network, env::SplitCTMRGEnv, alg::CTMRGAlgorithm) + alg isa SplitCTMRG || throw(ArgumentError("SplitCTMRGEnv requires SplitCTMRG.")) + eltype(network) <: PEPSSandwich || throw(ArgumentError("SplitCTMRG currently supports PEPS overlap networks only.")) + size(env)[2:3] == size(network) || throw(DimensionMismatch("Network and environment unit cells must match.")) + env.ket_first == (true, true, false, false) || throw(ArgumentError( + "SplitCTMRG requires ket_first == (true, true, false, false); got $(env.ket_first). " * + "For a rotated standard environment, apply the same rotl90, rotr90, or rot180 transformation " * + "to both the environment and the PEPS state/network to restore this orientation before calling " * + "`leading_boundary`, then rotate the returned environment back." + )) + alg.maxiter > 0 || throw(ArgumentError("maxiter must be positive.")) + return nothing +end +function check_input(::typeof(leading_boundary), network, env::CTMRGEnv, alg::SplitCTMRG) + throw(ArgumentError("SplitCTMRG requires a SplitCTMRGEnv; convert the initial environment first.")) +end + +function convergence_tensors(env::SplitCTMRGEnv, alg) + # The separate factors have an interlayer gauge freedom; compare normalized joined edges. + edges = map(E -> E / norm(E), CTMRGEnv(env).edges) + return env.corners, edges +end +network_value(network::InfiniteSquareNetwork, env::SplitCTMRGEnv) = network_value(network, CTMRGEnv(env)) + +"""Reflect and conjugate a PEPS overlap network to match the split-edge orientation.""" +function _reflect_split_network(network::InfiniteSquareNetwork{<:PEPSSandwich}, diagonal::Bool) + p = diagonal ? (4, 3, 2, 1) : (1, 4, 3, 2) + tensors = dtmap(unitcell(network)) do A + return map(t -> copy(permute(t', ((5,), p))), A) + end + tensors = diagonal ? permutedims(tensors) : reverse(tensors; dims = 2) + return InfiniteSquareNetwork(tensors) +end + +function ChainRulesCore.rrule(::typeof(_reflect_split_network), network::InfiniteSquareNetwork{<:PEPSSandwich}, diagonal::Bool) + reflected = _reflect_split_network(network, diagonal) + function reflection_pullback(Δ_) + Δ = unthunk(Δ_) + Δnetwork = Δ isa AbstractZero ? ZeroTangent() : _reflect_split_network(Δ, diagonal) + return NoTangent(), Δnetwork, NoTangent() + end + return reflected, reflection_pullback +end + +"""Orient a directional split move as a left move without refactorizing its edges.""" +function _orient_split_move(network, env, dir::Int; inverse::Bool = false) + if dir == WEST + return network, env + elseif dir == NORTH || dir == EAST + diagonal = dir == NORTH + return _reflect_split_network(network, diagonal), _reflect_split_env(env, diagonal) + else + rotate = inverse ? rotl90 : rotr90 + return rotate(network), rotate(env) + end +end + +function ctmrg_iteration(network::InfiniteSquareNetwork, env::SplitCTMRGEnv, alg::SplitCTMRG) + truncation_error = zero(real(scalartype(network))) + for dir in (WEST, NORTH, EAST, SOUTH) + oriented_network, oriented_env = _orient_split_move(network, env, dir) + for col in axes(oriented_network, 2) + oriented_env, info = ctmrg_leftmove(col, oriented_network, oriented_env, alg) + truncation_error = max(truncation_error, info.truncation_error) + end + _, env = _orient_split_move(oriented_network, oriented_env, dir; inverse = true) + end + return env, (; contraction_metrics = (; truncation_error)) +end + +function ctmrg_leftmove(col::Int, network, env::SplitCTMRGEnv, alg::SplitCTMRG) + # Infer the contraction graphs for their concrete tensor types, independently of the sweep. + projectors, info = Base.inferencebarrier(sequential_projectors)(col, network, env, alg.projector_alg, alg.interlayer_trunc) + env′ = Base.inferencebarrier(renormalize_sequentially)(col, projectors, network, env)::typeof(env) + return env′, (; truncation_error = real(scalartype(network))(info.truncation_error)) +end + +function renormalize_sequentially(col::Int, projectors, network, env::SplitCTMRGEnv) + (; green, teal, yellow, red) = projectors + updates = dtmap(axes(env, 2)) do r + nextrow = _next(r, size(env, 2)) + NW = renormalize_northwest_corner( + corner(env, NORTHWEST, r, col - 1), edge(env, NORTH, r, col), + green[2][nextrow], teal[2][nextrow], + ) + SW = renormalize_southwest_corner( + corner(env, SOUTHWEST, r, col - 1), edge(env, SOUTH, r, col), + green[1][r], teal[1][r], + ) + K = renormalize_west_edge( + edge(env, :ket, WEST, r, col - 1), green[1][r], teal[1][r], + yellow[2][r], red[2][r], ket(network[r, col]), + ) + B = renormalize_west_edge( + edge(env, :bra, WEST, r, col - 1), green[2][nextrow], teal[2][nextrow], + yellow[1][r], red[1][r], bra(network[r, col]), Val(:bra), + ) + return ( + _normalize_split_tensor(NW, corner(env, NORTHWEST, r, col)), + _normalize_split_tensor(SW, corner(env, SOUTHWEST, r, col)), + _normalize_split_tensor(K, edge(env, :ket, WEST, r, col)), + _normalize_split_tensor(B, edge(env, :bra, WEST, r, col)), + ) + end + corners = map(CartesianIndices(env.corners)) do I + dir, r, c = Tuple(I) + c == col && dir == NORTHWEST && return updates[r][1] + c == col && dir == SOUTHWEST && return updates[r][2] + return env.corners[I] + end + edges_ket = map(CartesianIndices(env.edges_ket)) do I + dir, r, c = Tuple(I) + return c == col && dir == WEST ? updates[r][3] : env.edges_ket[I] + end + edges_bra = map(CartesianIndices(env.edges_bra)) do I + dir, r, c = Tuple(I) + return c == col && dir == WEST ? updates[r][4] : env.edges_bra[I] + end + return SplitCTMRGEnv(corners, edges_ket, edges_bra, env.ket_first) +end + +"""Normalize an updated tensor and align its global phase with the previous factor.""" +function _normalize_split_tensor(t, previous) + t = t / norm(t) + space(t) == space(previous) || return t + overlap = dot(t, previous) + phase = iszero(overlap) ? one(overlap) : overlap / abs(overlap) + return t * phase +end diff --git a/src/algorithms/expectation_value/expectation_value.jl b/src/algorithms/expectation_value/expectation_value.jl index d2ea498f6..07d924386 100644 --- a/src/algorithms/expectation_value/expectation_value.jl +++ b/src/algorithms/expectation_value/expectation_value.jl @@ -20,6 +20,9 @@ function MPSKit.expectation_value( return sum(term_vals) end MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) +function MPSKit.expectation_value(bra::InfinitePEPS, O::LocalOperator, ket::InfinitePEPS, env::SplitCTMRGEnv) + return expectation_value(bra, O, ket, CTMRGEnv(env)) +end function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) checklattice(state, O) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly diff --git a/src/algorithms/expectation_value/network_value.jl b/src/algorithms/expectation_value/network_value.jl index 0bac1e033..2f5152cda 100644 --- a/src/algorithms/expectation_value/network_value.jl +++ b/src/algorithms/expectation_value/network_value.jl @@ -18,6 +18,7 @@ network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state) function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) return network_value(InfiniteSquareNetwork(peps), env) end +LinearAlgebra.norm(peps::InfinitePEPS, env::SplitCTMRGEnv) = network_value(InfiniteSquareNetwork(peps), env) # Local tensor insertions # ----------------------- diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl index d3a0fbfe3..aaa1f13ee 100644 --- a/src/algorithms/optimization/implicit_differentiation.jl +++ b/src/algorithms/optimization/implicit_differentiation.jl @@ -210,16 +210,24 @@ function _set_fixed_truncation(alg::CTMRGAlgorithm) return alg_fixed end +"""Freeze both outer and interlayer spaces in the split fixed-point iteration.""" +function _set_fixed_truncation(alg::SplitCTMRG) + alg_fixed = @set alg.projector_alg = _set_truncation(alg.projector_alg, FixedSpaceTruncation()) + return @set alg_fixed.interlayer_trunc = FixedSpaceTruncation() +end + # compute the CTMRG gradient through fixed-point differentiation function _rrule( gradmode::FixedPointGradient, config::RuleConfig, ::typeof(MPSKit.leading_boundary), - envinit::CTMRGEnv, + envinit::Union{CTMRGEnv, SplitCTMRGEnv}, state, alg::CTMRGAlgorithm, ) _check_algorithm_combination(alg, gradmode) + envinit isa SplitCTMRGEnv && scalartype(envinit) <: Real && + throw(ArgumentError("Split fixed-point gradients require a complex environment; use complex(envinit).")) env, = leading_boundary(envinit, state, alg) @@ -418,6 +426,9 @@ function _check_algorithm_combination(::SequentialCTMRG, ::ImplicitGradient) SimultaneousCTMRG instead to use :ImplicitGradient" throw(ArgumentError(msg)) end +function _check_algorithm_combination(::SplitCTMRG, ::ImplicitGradient) + throw(ArgumentError("SplitCTMRG supports FixedPointGradient; ImplicitGradient is not implemented.")) +end function _check_algorithm_combination(::SimultaneousCTMRG{<:FullInfiniteProjector}, ::ImplicitGradient) msg = "The `:ImplicitGradient` algorithm is not yet compatible with the `FullInfiniteProjector` scheme, \ diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 634292d3e..f39004c81 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -330,7 +330,7 @@ end Check compatibility of an initial PEPS and environment with a specified PEPS optimization algorithm. """ function check_input(::typeof(fixedpoint), peps₀, env₀, alg::PEPSOptimize) - if typeof(parent_alg(alg.boundary_alg)) <: Union{SimultaneousCTMRG, SequentialCTMRG} && + if typeof(parent_alg(alg.boundary_alg)) <: Union{SimultaneousCTMRG, SequentialCTMRG, SplitCTMRG} && parent_alg(alg.gradient_alg) isa FixedPointGradient && scalartype(env₀) <: Real # :fixed mode gauge fixing is incompatible with real environments msg = "the provided real environment is incompatible with :fixed mode \ diff --git a/src/algorithms/optimization/preconditioning.jl b/src/algorithms/optimization/preconditioning.jl index 31494d7d1..71559e347 100644 --- a/src/algorithms/optimization/preconditioning.jl +++ b/src/algorithms/optimization/preconditioning.jl @@ -111,6 +111,10 @@ gradient norm and iteration count tracked in `tol_state`. """ peps_precondition(x, g, tol_state::NamedTuple, alg::Nothing) = g +function peps_precondition(x::Tuple{<:InfinitePEPS, <:SplitCTMRGEnv}, g, tol_state::NamedTuple, alg::LocalPreconditioner) + return peps_precondition((x[1], CTMRGEnv(x[2])), g, tol_state, alg) +end + function peps_precondition(x, g, tol_state::NamedTuple, alg::LocalPreconditioner) peps, env = x δ = alg.regularization * tol_state.gradnorm^2 / max(tol_state.iter, 1) diff --git a/src/algorithms/select_algorithm.jl b/src/algorithms/select_algorithm.jl index c121fb3fd..6324b71a6 100644 --- a/src/algorithms/select_algorithm.jl +++ b/src/algorithms/select_algorithm.jl @@ -158,3 +158,8 @@ function select_algorithm( return CTMRGAlgorithm(; alg, tol, verbosity, decomposition_alg, kwargs...) end + +function select_algorithm(::typeof(leading_boundary), env₀::SplitCTMRGEnv; alg = :SplitCTMRG, kwargs...) + alg === :SplitCTMRG || throw(ArgumentError("SplitCTMRGEnv requires SplitCTMRG.")) + return SplitCTMRG(; kwargs...) +end diff --git a/src/environments/split_ctmrg_environments.jl b/src/environments/split_ctmrg_environments.jl new file mode 100644 index 000000000..574098d0b --- /dev/null +++ b/src/environments/split_ctmrg_environments.jl @@ -0,0 +1,250 @@ +""" + SplitCTMRGEnv(network, virtual_spaces...; interlayer_space=nothing) + SplitCTMRGEnv(env::CTMRGEnv; trunc=notrunc()) + +Corner transfer-matrix environment with separate rank-three ket and bra edge tensors. +The corners and unit-cell coordinates follow [`CTMRGEnv`](@ref). +An additional interlayer bond connects the two factors of each edge. +By default, initialization uses the outer environment space for this bond; conversion from an existing `CTMRGEnv` is exact unless `trunc` is specified. + +Use `edge(env, :ket, dir, r, c)` and `edge(env, :bra, dir, r, c)` to access the factors, or `CTMRGEnv(env)` to contract them for existing observable routines. +This environment supports [`SplitCTMRG`](@ref) and its fixed-point gradient. +""" +struct SplitCTMRGEnv{C, E} + corners::Array{C, 3} + edges_ket::Array{E, 3} + edges_bra::Array{E, 3} + # Whether the ket factor carries the codomain environment leg of the joined edge. + ket_first::NTuple{4, Bool} + function SplitCTMRGEnv( + corners::Array{C, 3}, edges_ket::Array{E, 3}, edges_bra::Array{E, 3}, + ket_first::NTuple{4, Bool} = (true, true, false, false), + ) where {C, E} + size(corners) == size(edges_ket) == size(edges_bra) || + throw(DimensionMismatch("Corner and split-edge arrays must have equal sizes.")) + size(corners, 1) == 4 || throw(DimensionMismatch("Expected four environment directions.")) + return new{C, E}(corners, edges_ket, edges_bra, ket_first) + end +end + +function ChainRulesCore.rrule( + ::Type{SplitCTMRGEnv}, corners::Array{C, 3}, edges_ket::Array{E, 3}, edges_bra::Array{E, 3}, + ket_first::NTuple{4, Bool}, + ) where {C, E} + env = SplitCTMRGEnv(corners, edges_ket, edges_bra, ket_first) + function split_env_pullback(Δenv_) + Δenv = unthunk(Δenv_) + Δenv isa AbstractZero && return NoTangent(), ZeroTangent(), ZeroTangent(), ZeroTangent(), NoTangent() + return NoTangent(), Δenv.corners, Δenv.edges_ket, Δenv.edges_bra, NoTangent() + end + return env, split_env_pullback +end +function ChainRulesCore.rrule( + ::Type{SplitCTMRGEnv}, corners::Array{C, 3}, edges_ket::Array{E, 3}, edges_bra::Array{E, 3}, + ) where {C, E} + env, pullback = rrule(SplitCTMRGEnv, corners, edges_ket, edges_bra, (true, true, false, false)) + return env, Δ -> Base.front(pullback(Δ)) +end +function ChainRulesCore.rrule(::typeof(getproperty), env::SplitCTMRGEnv, name::Symbol) + result = getfield(env, name) + name === :ket_first && return result, _ -> (NoTangent(), ZeroTangent(), NoTangent()) + function split_property_pullback(Δ_) + Δ = unthunk(Δ_) + Δ isa AbstractZero && return NoTangent(), ZeroTangent(), NoTangent() + Δ = _project_split_array(Δ, result) + corners = name === :corners ? Δ : zerovector.(env.corners) + edges_ket = name === :edges_ket ? Δ : zerovector.(env.edges_ket) + edges_bra = name === :edges_bra ? Δ : zerovector.(env.edges_bra) + return NoTangent(), SplitCTMRGEnv(corners, edges_ket, edges_bra, env.ket_first), NoTangent() + end + return result, split_property_pullback +end + +"""Project array cotangents and materialize zero entries in the primal tensor storage.""" +function _project_split_array(Δ, primal) + return map(primal, Δ) do tensor, delta + delta = unthunk(delta) + return delta isa AbstractZero ? zerovector(tensor) : ProjectTo(tensor)(delta) + end +end + +"""Factor an edge across its interlayer bond, distributing the singular values symmetrically.""" +function _split_edge(E::CTMRGEdgeTensor{T, S, 3}, ket_first::Bool, trunc::TruncationStrategy) where {T, S} + p = ket_first ? ((1, 2), (4, 3)) : ((1, 3), (4, 2)) + U, s, V, = svd_trunc(permute(E, p); trunc) + sqs = sdiag_pow(s, 0.5) + left = U * sqs + right = permute(sqs * V, ((1, 3), (2,))) + return ket_first ? (left, right) : (right, left) +end + +function SplitCTMRGEnv(env::CTMRGEnv; trunc::TruncationStrategy = notrunc()) + ket_first = (true, true, false, false) + factors = map(CartesianIndices(env.edges)) do I + return _split_edge(env.edges[I], ket_first[I[1]], trunc) + end + corners = map(C -> one(scalartype(env)) * TensorMap(C), env.corners) + return SplitCTMRGEnv(corners, map(first, factors), map(last, factors), ket_first) +end + +function SplitCTMRGEnv( + f, T, network::InfiniteSquareNetwork{<:PEPSSandwich}, virtual_spaces...; + interlayer_space = nothing, + ) + env = CTMRGEnv(f, T, network, virtual_spaces...) + ket_first = (true, true, false, false) + factors = map(CartesianIndices(env.edges)) do I + V = isnothing(interlayer_space) ? space(env.edges[I], 1) : interlayer_space + return _split_edge(env.edges[I], ket_first[I[1]], truncspace(V)) + end + return SplitCTMRGEnv(env.corners, map(first, factors), map(last, factors), ket_first) +end +function SplitCTMRGEnv(f, T, state::InfinitePEPS, virtual_spaces...; kwargs...) + return SplitCTMRGEnv(f, T, InfiniteSquareNetwork(state), virtual_spaces...; kwargs...) +end +function SplitCTMRGEnv(network::Union{InfinitePEPS, InfiniteSquareNetwork{<:PEPSSandwich}}, virtual_spaces...; kwargs...) + return SplitCTMRGEnv(randn, scalartype(network), network, virtual_spaces...; kwargs...) +end + +Base.size(env::SplitCTMRGEnv, args...) = size(env.corners, args...) +Base.axes(env::SplitCTMRGEnv, args...) = axes(env.corners, args...) +Base.@propagate_inbounds corner(env::SplitCTMRGEnv, I...) = periodic_getindex(env, env.corners, I) +Base.@propagate_inbounds setcorner!(env::SplitCTMRGEnv, v, I...) = periodic_setindex!(env, env.corners, v, I) + +""" + edge(env::SplitCTMRGEnv, layer::Symbol, dir, r, c) + +Return the `:ket` or `:bra` factor of an edge with periodic unit-cell indexing. +Omitting `layer` contracts the interlayer bond and returns a conventional CTMRG edge. +""" +Base.@propagate_inbounds function edge(env::SplitCTMRGEnv, layer::Symbol, I...) + data = layer === :ket ? env.edges_ket : layer === :bra ? env.edges_bra : + throw(ArgumentError("Unknown split edge layer: $layer")) + return periodic_getindex(env, data, I) +end +Base.@propagate_inbounds function setedge!(env::SplitCTMRGEnv, v, layer::Symbol, I...) + data = layer === :ket ? env.edges_ket : layer === :bra ? env.edges_bra : + throw(ArgumentError("Unknown split edge layer: $layer")) + return periodic_setindex!(env, data, v, I) +end + +"""Contract the interlayer bond of a split edge in the conventional index order.""" +function _join_split_edge(K, B, ket_first::Bool) + if ket_first + return @tensor E[l k b; r] := K[l k; i] * B[i b; r] + else + return @tensor E[l k b; r] := B[l b; i] * K[i k; r] + end +end +function edge(env::SplitCTMRGEnv, dir::Int, I...) + return _join_split_edge(edge(env, :ket, dir, I...), edge(env, :bra, dir, I...), env.ket_first[dir]) +end + +""" + CTMRGEnv(env::SplitCTMRGEnv) + +Contract each interlayer bond to obtain an environment for the existing observable routines. +""" +function CTMRGEnv(env::SplitCTMRGEnv) + edges = map(CartesianIndices(env.corners)) do I + return edge(env, Tuple(I)...) + end + return CTMRGEnv(copy.(env.corners), edges) +end +Base.convert(::Type{CTMRGEnv}, env::SplitCTMRGEnv) = CTMRGEnv(env) + +eachcoordinate(env::SplitCTMRGEnv) = collect(Iterators.product(axes(env, 2), axes(env, 3))) +eachcoordinate(env::SplitCTMRGEnv, dirs) = collect(Iterators.product(dirs, axes(env, 2), axes(env, 3))) +cornertype(::Type{SplitCTMRGEnv{C, E}}) where {C, E} = C +edgetype(::Type{SplitCTMRGEnv{C, E}}) where {C, E} = E +cornertype(env::SplitCTMRGEnv) = cornertype(typeof(env)) +edgetype(env::SplitCTMRGEnv) = edgetype(typeof(env)) +VI.scalartype(::Type{SplitCTMRGEnv{C, E}}) where {C, E} = promote_type(scalartype(C), scalartype(E)) +TensorKit.spacetype(::Type{E}) where {E <: SplitCTMRGEnv} = spacetype(cornertype(E)) +TensorKit.storagetype(::Type{E}) where {E <: SplitCTMRGEnv} = TensorKit.promote_storagetype(cornertype(E), edgetype(E)) + +"""Reorder and transform directional tensor arrays without mutating their storage.""" +function _transform_split_array(data, arraymap, tensormap, dirs) + slices = map(1:4) do dir + source = findfirst(==(dir), dirs) + slice = tensormap.(arraymap(data[source, :, :])) + return reshape(slice, 1, size(slice)...) + end + return cat(slices...; dims = 1) +end + +"""Transform unit-cell coordinates, directions, and tensors together without joining the edges.""" +function _transform_split_env(env::SplitCTMRGEnv, arraymap, tensormap, dirs, cornerdirs; reverse_order::Bool = false) + corners = _transform_split_array(env.corners, arraymap, tensormap, cornerdirs) + edges_ket = _transform_split_array(env.edges_ket, arraymap, tensormap, dirs) + edges_bra = _transform_split_array(env.edges_bra, arraymap, tensormap, dirs) + order = ntuple(i -> xor(env.ket_first[findfirst(==(i), dirs)], reverse_order), 4) + return SplitCTMRGEnv(corners, edges_ket, edges_bra, order) +end +Base.rotl90(env::SplitCTMRGEnv) = _transform_split_env(env, rotl90, identity, (4, 1, 2, 3), (4, 1, 2, 3)) +Base.rotr90(env::SplitCTMRGEnv) = _transform_split_env(env, rotr90, identity, (2, 3, 4, 1), (2, 3, 4, 1)) +Base.rot180(env::SplitCTMRGEnv) = rotl90(rotl90(env)) + +"""Reflect and conjugate an environment tensor, preserving nondual outer environment spaces.""" +_reflect_split_tensor(C::CTMRGCornerTensor) = copy(C') +_reflect_split_tensor(E::CTMRGEdgeTensor{<:Any, <:Any, 2}) = copy(permute(E', ((1, 3), (2,)))) + +"""Reflect and conjugate a split environment across the diagonal or the vertical axis.""" +function _reflect_split_env(env::SplitCTMRGEnv, diagonal::Bool) + arraymap = diagonal ? permutedims : A -> reverse(A; dims = 2) + dirs = diagonal ? (4, 3, 2, 1) : (1, 4, 3, 2) + cornerdirs = diagonal ? (1, 4, 3, 2) : (2, 1, 4, 3) + return _transform_split_env(env, arraymap, _reflect_split_tensor, dirs, cornerdirs; reverse_order = true) +end + +"""Map corresponding tensors while preserving the split-edge orientation.""" +function _map_split_env(f, env::SplitCTMRGEnv, others::SplitCTMRGEnv...) + return SplitCTMRGEnv( + map(f, env.corners, map(e -> e.corners, others)...), + map(f, env.edges_ket, map(e -> e.edges_ket, others)...), + map(f, env.edges_bra, map(e -> e.edges_bra, others)...), env.ket_first, + ) +end + +"""Apply an in-place operation to corresponding split-environment tensors.""" +function _foreach_split_env!(f, env::SplitCTMRGEnv, others::SplitCTMRGEnv...) + foreach(f, env.corners, map(e -> e.corners, others)...) + foreach(f, env.edges_ket, map(e -> e.edges_ket, others)...) + foreach(f, env.edges_bra, map(e -> e.edges_bra, others)...) + return env +end + +function update!(env::SplitCTMRGEnv, other::SplitCTMRGEnv) + env.ket_first == other.ket_first || throw(ArgumentError("Split-edge orientations must match.")) + env.corners .= other.corners + env.edges_ket .= other.edges_ket + env.edges_bra .= other.edges_bra + return env +end + +Base.:+(a::SplitCTMRGEnv, b::SplitCTMRGEnv) = VI.add(a, b, 1, 1) +Base.:-(a::SplitCTMRGEnv, b::SplitCTMRGEnv) = VI.add(a, b, -1, 1) +Base.:*(α::Number, env::SplitCTMRGEnv) = VI.scale(env, α) +Base.similar(env::SplitCTMRGEnv) = _map_split_env(similar, env) +Base.real(env::SplitCTMRGEnv) = _map_split_env(real, env) +Base.complex(env::SplitCTMRGEnv) = _map_split_env(complex, env) +function Base.isapprox(a::SplitCTMRGEnv, b::SplitCTMRGEnv; kwargs...) + return a.ket_first == b.ket_first && + isapprox(a.corners, b.corners; kwargs...) && + isapprox(a.edges_ket, b.edges_ket; kwargs...) && + isapprox(a.edges_bra, b.edges_bra; kwargs...) +end + +VI.zerovector(env::SplitCTMRGEnv, ::Type{T}) where {T <: Number} = _map_split_env(Base.Fix2(zerovector, T), env) +VI.zerovector!(env::SplitCTMRGEnv) = _foreach_split_env!(zerovector!, env) +VI.zerovector!!(env::SplitCTMRGEnv) = zerovector!(env) +VI.scale(env::SplitCTMRGEnv, α::Number) = _map_split_env(Base.Fix2(scale, α), env) +VI.scale!(env::SplitCTMRGEnv, α::Number) = _foreach_split_env!(Base.Fix2(scale!, α), env) +VI.scale!(dst::SplitCTMRGEnv, src::SplitCTMRGEnv, α::Number) = _foreach_split_env!((x, y) -> scale!(x, y, α), dst, src) +VI.scale!!(env::SplitCTMRGEnv, α::Number) = scale!(env, α) +VI.scale!!(dst::SplitCTMRGEnv, src::SplitCTMRGEnv, α::Number) = scale!(dst, src, α) +VI.add(a::SplitCTMRGEnv, b::SplitCTMRGEnv, α::Number, β::Number) = _map_split_env((x, y) -> add(x, y, α, β), a, b) +VI.add!(a::SplitCTMRGEnv, b::SplitCTMRGEnv, α::Number, β::Number) = _foreach_split_env!((x, y) -> add!(x, y, α, β), a, b) +VI.add!!(a::SplitCTMRGEnv, b::SplitCTMRGEnv, α::Number, β::Number) = add!(a, b, α, β) +VI.inner(a::SplitCTMRGEnv, b::SplitCTMRGEnv) = inner((a.corners, a.edges_ket, a.edges_bra), (b.corners, b.edges_ket, b.edges_bra)) +VI.norm(env::SplitCTMRGEnv) = norm((env.corners, env.edges_ket, env.edges_bra)) diff --git a/src/networks/infinitesquarenetwork.jl b/src/networks/infinitesquarenetwork.jl index fe62377b3..ccb4e1b57 100644 --- a/src/networks/infinitesquarenetwork.jl +++ b/src/networks/infinitesquarenetwork.jl @@ -117,6 +117,18 @@ end ## Chainrules +"""Materialize unused sandwich layers and project cotangents to their tensor storage.""" +function _project_network_cotangent(delta, primal) + delta = unthunk(delta) + delta isa AbstractZero && return zerovector(primal) + primal isa Tuple || return ProjectTo(primal)(delta) + components = delta isa Tangent ? ChainRulesCore.backing(delta) : delta + return map(primal, components) do tensor, component + component = unthunk(component) + return component isa AbstractZero ? zerovector(tensor) : ProjectTo(tensor)(component) + end +end + # generic implementation function ChainRulesCore.rrule( ::typeof(Base.getindex), network::InfiniteSquareNetwork, r::Int, c::Int @@ -124,10 +136,9 @@ function ChainRulesCore.rrule( O = network[r, c] function getindex_pullback(ΔO_) - ΔO = map(unthunk, ΔO_) - if ΔO isa Tangent - ΔO = ChainRulesCore.construct(typeof(O), ChainRulesCore.backing(ΔO)) - end + ΔO = unthunk(ΔO_) + ΔO isa AbstractZero && return NoTangent(), ZeroTangent(), NoTangent(), NoTangent() + ΔO = _project_network_cotangent(ΔO, O) Δnetwork = zerovector(network) Δnetwork[r, c] = ΔO return NoTangent(), Δnetwork, NoTangent(), NoTangent() diff --git a/test/ctmrg/split.jl b/test/ctmrg/split.jl new file mode 100644 index 000000000..bc80b15e3 --- /dev/null +++ b/test/ctmrg/split.jl @@ -0,0 +1,7 @@ +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +TestSuite.ctmrg_split_projectors(Vector) +TestSuite.ctmrg_split_environment(Vector) +TestSuite.ctmrg_split_unitcell(Vector) +TestSuite.ctmrg_split_expectation_values(Vector) diff --git a/test/gradients/split_ctmrg_gradients.jl b/test/gradients/split_ctmrg_gradients.jl new file mode 100644 index 000000000..5f0bb2611 --- /dev/null +++ b/test/gradients/split_ctmrg_gradients.jl @@ -0,0 +1,6 @@ +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +TestSuite.gradients_split_gauge(Vector) +TestSuite.gradients_split_sweep(Vector) +TestSuite.gradients_split_fixedpoint(Vector) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..fa550601a 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -83,6 +83,13 @@ using .CompressLocal # CTMRG # ----- +module CTMRGSplit + include("ctmrg/split.jl") + export ctmrg_split_projectors + export ctmrg_split_environment, ctmrg_split_unitcell, ctmrg_split_expectation_values +end +using .CTMRGSplit + module CTMRGContractions include("ctmrg/contractions.jl") export ctmrg_contractions_random_cartesian_spaces, ctmrg_contractions_specific_u1_spaces @@ -197,6 +204,12 @@ module GradientsCTMRG end using .GradientsCTMRG +module GradientsSplitCTMRG + include("gradients/split_ctmrg_gradients.jl") + export gradients_split_gauge, gradients_split_sweep, gradients_split_fixedpoint +end +using .GradientsSplitCTMRG + # Time evolution # -------------- module TimeEvolClusterProjectors diff --git a/test/testsuite/ctmrg/split.jl b/test/testsuite/ctmrg/split.jl new file mode 100644 index 000000000..4f9948448 --- /dev/null +++ b/test/testsuite/ctmrg/split.jl @@ -0,0 +1,96 @@ +using Test, Random, LinearAlgebra +using TensorKit, MatrixAlgebraKit, PEPSKit, Adapt +using PEPSKit: ctmrg_iteration, WEST + +function ctmrg_split_projectors(AT) + Random.seed!(37913) + return @testset "Split projector stability ($AT)" begin + alg = SplitCTMRG(; verbosity = 0).projector_alg + trunc = truncspace(ComplexSpace(3)) + # Cover the two-index yellow and three-index red cuts. + for cut in (ComplexSpace(3) ⊗ ComplexSpace(2), ComplexSpace(3) ⊗ ComplexSpace(2) ⊗ ComplexSpace(2)) + Q = adapt(AT, randn(ComplexF64, ComplexSpace(12) ← cut)) + Q_next = adapt(AT, randn(ComplexF64, cut ← ComplexSpace(12))) + reference, refinfo = PEPSKit.compute_projector((Q, Q_next), alg; trunc) + compressed, info = PEPSKit.compute_projector((Q, Q_next), alg; trunc, compress = true) + # QR compression preserves the cut projector and the nonzero singular spectrum. + @test compressed[1] * compressed[2] ≈ reference[1] * reference[2] + @test info.S ≈ refinfo.S + for compress in (false, true), scale in (1.0e200, 1.0e-200) + scaled, = PEPSKit.compute_projector((scale * Q, scale * Q_next), alg; trunc, compress) + @test scaled[1] * scaled[2] ≈ reference[1] * reference[2] + end + end + end +end + +function ctmrg_split_environment(AT) + return @testset "Split environment ($AT)" begin + Random.seed!(47119) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell = (2, 3))) + env = CTMRGEnv(psi, ComplexSpace(3)) + split = SplitCTMRGEnv(env) + @test CTMRGEnv(split) ≈ env + diagonal = CTMRGEnv(DiagonalTensorMap.(real.(env.corners)), copy.(env.edges)) + @test CTMRGEnv(SplitCTMRGEnv(diagonal)) ≈ diagonal + @test CTMRGEnv(rotl90(split)) ≈ rotl90(env) + rescaled = deepcopy(split) + rescaled.edges_ket .*= 2 + rescaled.edges_bra ./= 2 + @test PEPSKit.calc_convergence(rescaled, PEPSKit.convergence_spectra(split, SplitCTMRG())..., SplitCTMRG())[1] < 1.0e-12 + end +end + +function ctmrg_split_unitcell(AT) + Random.seed!(72893) + return @testset "Split sweeps ($AT)" begin + unitcell = (2, 3) + P = ComplexSpace.(rand(2:3, unitcell...)) + N, E = ComplexSpace.(rand(2:3, unitcell...)), ComplexSpace.(rand(2:3, unitcell...)) + chis = ntuple(_ -> ComplexSpace.(rand(3:5, unitcell...)), 4) + psi = adapt(AT, InfinitePEPS(randn, ComplexF64, P, N, E)) + env = SplitCTMRGEnv(psi, chis...; interlayer_space = ComplexSpace(3)) + alg = SplitCTMRG(; verbosity = 0, projector_alg = :SplitProjector) + + # Non-self-dual charge content catches incorrect conjugations in reflected moves. + V = U1Space(-1 => 1, 0 => 1, 1 => 1) + chi = U1Space(-1 => 1, 0 => 2, 1 => 2) + psi_u1 = adapt(AT, InfinitePEPS(randn, ComplexF64, U1Space(0 => 1, 1 => 1), V)) + env_u1 = SplitCTMRGEnv(psi_u1, chi) + @test CTMRGEnv(SplitCTMRGEnv(CTMRGEnv(env_u1))) ≈ CTMRGEnv(env_u1) + for (state, initial) in ((psi, env), (psi_u1, env_u1)) + final, = ctmrg_iteration(InfiniteSquareNetwork(state), initial, alg) + for field in (:corners, :edges_ket, :edges_bra) + @test space.(getproperty(final, field)) == space.(getproperty(initial, field)) + end + end + + # The interlayer and outer refinement parameters can be set independently. + resized, = ctmrg_iteration( + InfiniteSquareNetwork(psi), env, + SplitCTMRG(; verbosity = 0, trunc = truncspace(ComplexSpace(4)), interlayer_trunc = truncspace(ComplexSpace(2))), + ) + @test all(K -> dim(codomain(K)[1]) == 2, resized.edges_ket[WEST, :, :]) + @test all(K -> dim(domain(K)[1]) == 4, resized.edges_ket[WEST, :, :]) + end +end + +function ctmrg_split_expectation_values(AT) + return @testset "Split expectation values ($T, $AT)" for T in (Float64, ComplexF64) + Random.seed!(84109) + psi = adapt(AT, product_peps(randn, T, ComplexSpace(2), ComplexSpace(2); noise_amp = 0.2, state_vector = fill([1.0, 0.0], (1, 1)))) + chi = ComplexSpace(4) + initial = CTMRGEnv(psi, chi) + reference, refinfo = leading_boundary( + initial, psi, SequentialCTMRG(; verbosity = 0, tol = 1.0e-7, maxiter = 40, projector_alg = :FullInfiniteProjector)::SequentialCTMRG, + ) + split_initial = SplitCTMRGEnv(initial; trunc = truncspace(chi)) + split, info = leading_boundary(split_initial, psi; verbosity = 0, tol = 1.0e-7, maxiter = 40) + @test refinfo.converged + @test info.converged + joined = CTMRGEnv(split) + @test abs(norm(psi, joined)) ≈ abs(norm(psi, reference)) rtol = 1.0e-4 + H = adapt(AT, heisenberg_XYZ(InfiniteSquare())) + @test expectation_value(psi, H, joined) ≈ expectation_value(psi, H, reference) atol = 1.0e-4 + end +end diff --git a/test/testsuite/gradients/split_ctmrg_gradients.jl b/test/testsuite/gradients/split_ctmrg_gradients.jl new file mode 100644 index 000000000..63d82e81d --- /dev/null +++ b/test/testsuite/gradients/split_ctmrg_gradients.jl @@ -0,0 +1,84 @@ +using Test, Random, LinearAlgebra +using TensorKit, PEPSKit, Zygote, Adapt +using PEPSKit: ctmrg_iteration, hook_pullback + +function gradients_split_gauge(AT) + return @testset "Split outer and interlayer gauge ($AT)" begin + Random.seed!(32591) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell = (2, 3))) + env = SplitCTMRGEnv(psi, ComplexSpace(2); interlayer_space = ComplexSpace(3)) + pairs = map(PEPSKit.eachcoordinate(env, 1:4)) do (dir, r, c) + firstlayer = env.ket_first[dir] ? :ket : :bra + factor = edge(env, firstlayer, dir, r, c) + outer, interlayer = codomain(factor)[1], domain(factor)[1] + return first(left_orth(adapt(AT, randn(ComplexF64, outer ← outer)))), + first(left_orth(adapt(AT, randn(ComplexF64, interlayer ← interlayer)))) + end + signs = (; outer = map(first, pairs), interlayer = map(last, pairs)) + transformed = PEPSKit.fix_relative_phases(env, signs) + restored = gauge_fix(transformed, env, PEPSKit.ScramblingEnvGauge()) + @test norm(restored - env) / norm(env) < 1.0e-10 + end +end + +"""Prepare a small entangled PEPS and split environment for gradient checks.""" +function _split_gradient_inputs(AT) + Random.seed!(84109) + psi = adapt( + AT, product_peps( + randn, ComplexF64, ComplexSpace(2), ComplexSpace(2); + noise_amp = 0.2, state_vector = fill([1.0, 0.0], (1, 1)) + ) + ) + env = SplitCTMRGEnv(psi, ComplexSpace(2)) + H = adapt(AT, heisenberg_XYZ(InfiniteSquare())) + alg = SplitCTMRG(; verbosity = 0, tol = 1.0e-10, maxiter = 100) + return psi, env, H, alg +end + +function gradients_split_sweep(AT) + return @testset "Split sweep pullback ($AT)" begin + psi, env, H, alg = _split_gradient_inputs(AT) + loss(A, x) = real(expectation_value(A, H, ctmrg_iteration(InfiniteSquareNetwork(A), x, alg)[1])) + _, back = Zygote.pullback(loss, psi, env) + dpsi, denv = back(1.0) + direction = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) + direction_env = SplitCTMRGEnv(psi, ComplexSpace(2)) + epsilon = 1.0e-6 + finite_state = (loss(psi + epsilon * direction, env) - loss(psi - epsilon * direction, env)) / (2 * epsilon) + finite_env = (loss(psi, env + epsilon * direction_env) - loss(psi, env - epsilon * direction_env)) / (2 * epsilon) + @test real(dot(dpsi, direction)) ≈ finite_state atol = 1.0e-5 rtol = 1.0e-4 + @test real(inner(denv, direction_env)) ≈ finite_env atol = 1.0e-5 rtol = 1.0e-4 + end +end + +function gradients_split_fixedpoint(AT) + return @testset "Split fixed-point energy gradient ($AT)" begin + psi, initial, H, alg = _split_gradient_inputs(AT) + env, info = leading_boundary(initial, psi, alg) + @test info.converged + next, = ctmrg_iteration(InfiniteSquareNetwork(psi), env, alg) + aligned = gauge_fix(next, env, PEPSKit.ScramblingEnvGauge()) + @test norm(aligned - env) / norm(env) < 1.0e-7 + gradalg = FixedPointGradient(; solver_alg = (; alg = :GMRES, tol = 1.0e-8, maxiter = 100)) + E, gs = Zygote.withgradient(psi) do A + boundary, = hook_pullback(leading_boundary, env, A, alg; alg_rrule = gradalg) + return cost_function(A, boundary, H) + end + direction = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) + epsilon = 1.0e-5 + energy(A) = cost_function(A, leading_boundary(env, A, alg)[1], H) + finite = (energy(psi + epsilon * direction) - energy(psi - epsilon * direction)) / (2 * epsilon) + @test real(dot(only(gs), direction)) ≈ finite atol = 1.0e-4 rtol = 1.0e-3 + optalg = PEPSOptimize(; + boundary_alg = alg, + optimizer_alg = (; maxiter = 2, verbosity = 0) + ) + optimized_psi, final, optimized, = fixedpoint(H, psi, deepcopy(env), optalg) + @test optimized < E + @test cost_function(optimized_psi, final, H) ≈ optimized + next, = ctmrg_iteration(InfiniteSquareNetwork(optimized_psi), final, alg) + aligned = gauge_fix(next, final, PEPSKit.ScramblingEnvGauge()) + @test norm(aligned - final) / norm(final) < 1.0e-7 + end +end From 70c1517923b5bc963a1eedd488032e8261cf6aee Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 12 Oct 2026 17:08:40 +0800 Subject: [PATCH 02/14] Deduplicate CTMRG gauge fixing --- docs/src/man/environments.md | 1 + src/algorithms/ctmrg/gaugefix.jl | 90 +++++++++++++------------- src/algorithms/ctmrg/gaugefix_split.jl | 39 ++++------- 3 files changed, 57 insertions(+), 73 deletions(-) diff --git a/docs/src/man/environments.md b/docs/src/man/environments.md index 874b84e3f..3118c338c 100644 --- a/docs/src/man/environments.md +++ b/docs/src/man/environments.md @@ -26,6 +26,7 @@ This version supports PEPS overlap networks, full projectors, and automatic diff For comparisons with variPEPS, select its `Projector_Method.FULL` option; its default `FISHMAN` construction reweights the environment halves and uses different projectors. Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. +The ordinary and split environments use the same boundary-MPS gauge method, with each split edge contributing two MPS sites instead of one. Both sets of bond spaces are held fixed in the differentiated iteration, including when the forward run uses other truncation strategies. Use a complex environment for fixed-point gradients, as for conventional asymmetric CTMRG. The first version contracts split edges into `CTMRGEnv` for energy evaluation and local preconditioning, preserving the derivative through both factors: diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index 6ffb7d7da..f5cc9feb9 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -7,6 +7,7 @@ https://arxiv.org/abs/2311.11894. This works by constructing a transfer matrix c of an edge tensor and a random MPS, thus scrambling potential degeneracies, and then performing a QR decomposition to extract the gauge signs. This is adapted accordingly for asymmetric CTMRG algorithms using multi-site unit cell transfer matrices. +Ordinary edges form one MPS site each; split edges form two alternating sites with outer and interlayer bonds. """ struct ScramblingEnvGauge end @@ -27,60 +28,62 @@ Given that the CTMRG run is converged, the returned environment will be element-wise converged to `envprev`. """ function gauge_fix( - envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, alg::G - ) where {C, T, G <: Union{ScramblingEnvGauge, ScramblingEnvGaugeC4v}} + envfinal::Union{CTMRGEnv, SplitCTMRGEnv}, envprev::Union{CTMRGEnv, SplitCTMRGEnv}, + alg::Union{ScramblingEnvGauge, ScramblingEnvGaugeC4v} + ) signs, corner_phases, edge_phases = compute_gauge_fix_gauge(envfinal, envprev, alg) return fix_phases(envfinal, signs, corner_phases, edge_phases) end function compute_gauge_fix_gauge( - envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, alg::G - ) where {C, T, G <: Union{ScramblingEnvGauge, ScramblingEnvGaugeC4v}} - # Check if spaces in envprev and envfinal are the same + envfinal::Union{CTMRGEnv, SplitCTMRGEnv}, envprev::Union{CTMRGEnv, SplitCTMRGEnv}, + alg::Union{ScramblingEnvGauge, ScramblingEnvGaugeC4v} + ) + _check_gauge_compatibility(envfinal, envprev) + signs = compute_relative_phases(envfinal, envprev, alg) + corner_phases, edge_phases = compute_global_phases( + fix_relative_phases(envfinal, signs), envprev + ) + return signs, corner_phases, edge_phases +end + +"""Check that ordinary environments have matching spaces for gauge alignment.""" +function _check_gauge_compatibility(envfinal::CTMRGEnv, envprev::CTMRGEnv) same_spaces = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c) space(envfinal.edges[dir, r, c]) == space(envprev.edges[dir, r, c]) && space(envfinal.corners[dir, r, c]) == space(envprev.corners[dir, r, c]) end all(same_spaces) || throw(ArgumentError("Spaces of envfinal and envprev are not the same")) - # find relative phases - signs = compute_relative_phases(envfinal, envprev, alg) + return nothing +end - # find additional global phases - corner_phases, edge_phases = compute_global_phases( - fix_relative_phases(envfinal, signs), envprev - ) +"""Apply relative gauge rotations followed by global phase corrections.""" +function fix_phases(env::Union{CTMRGEnv, SplitCTMRGEnv}, signs, corner_phases, edge_phases) + return fix_global_phases(fix_relative_phases(env, signs), corner_phases, edge_phases) +end - return signs, corner_phases, edge_phases +"""Return the row and column increments for clockwise boundary traversal.""" +_boundary_step(dir::Int) = dir == EAST ? (1, 0) : dir == WEST ? (-1, 0) : dir == NORTH ? (0, 1) : (0, -1) + +"""List boundary edge coordinates clockwise, starting at the specified outer bond.""" +function _boundary_coordinates(env::Union{CTMRGEnv, SplitCTMRGEnv}, (dir, r, c)) + dr, dc = _boundary_step(dir) + count = iszero(dr) ? size(env, 3) : size(env, 2) + return [(dir, mod1(r + i * dr, size(env, 2)), mod1(c + i * dc, size(env, 3))) for i in 0:(count - 1)] end -function fix_phases(env::CTMRGEnv, signs, corner_phases, edge_phases) - env_fixed_relative = fix_relative_phases(env, signs) - env_fixed_global = fix_global_phases(env_fixed_relative, corner_phases, edge_phases) - return env_fixed_global +"""Represent an ordinary boundary as an MPS with one edge tensor per site.""" +function _boundary_chain(env::CTMRGEnv, coordinate) + return [edge(env, I...) for I in _boundary_coordinates(env, coordinate)] end function compute_relative_phases( envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, ::ScramblingEnvGauge ) where {C, T} - signs = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c) - # Gather edge tensors and pretend they're InfiniteMPSs - if dir == NORTH - Tsprev = circshift(envprev.edges[dir, r, :], 1 - c) - Tsfinal = circshift(envfinal.edges[dir, r, :], 1 - c) - elseif dir == EAST - Tsprev = circshift(envprev.edges[dir, :, c], 1 - r) - Tsfinal = circshift(envfinal.edges[dir, :, c], 1 - r) - elseif dir == SOUTH - Tsprev = circshift(reverse(envprev.edges[dir, r, :]), c) - Tsfinal = circshift(reverse(envfinal.edges[dir, r, :]), c) - elseif dir == WEST - Tsprev = circshift(reverse(envprev.edges[dir, :, c]), r) - Tsfinal = circshift(reverse(envfinal.edges[dir, :, c]), r) - end - - return _scrambling_gauge(Tsfinal, Tsprev) + signs = map(eachcoordinate(envfinal, 1:4)) do coordinate + return _scrambling_gauge(_boundary_chain(envfinal, coordinate), _boundary_chain(envprev, coordinate)) end return signs @@ -160,21 +163,16 @@ function right_transfermatrix_fixedpoint( return first(vecs) end -# Explicit fixing of relative phases (doing this compactly in a loop is annoying) -function fix_relative_phases(envfinal::CTMRGEnv, signs) - corners_fixed = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c) - Cf = if dir == NORTHWEST - fix_gauge_northwest_corner((r, c), envfinal, signs) - elseif dir == NORTHEAST - fix_gauge_northeast_corner((r, c), envfinal, signs) - elseif dir == SOUTHEAST - fix_gauge_southeast_corner((r, c), envfinal, signs) - else # dir == SOUTHWEST - fix_gauge_southwest_corner((r, c), envfinal, signs) - end - return Cf +"""Apply outer boundary gauges to the four corner tensors at each unit-cell site.""" +function _fix_gauge_corners(env::Union{CTMRGEnv, SplitCTMRGEnv}, signs) + corner_fns = (fix_gauge_northwest_corner, fix_gauge_northeast_corner, fix_gauge_southeast_corner, fix_gauge_southwest_corner) + return map(eachcoordinate(env, 1:4)) do (dir, r, c) + return corner_fns[dir]((r, c), env, signs) end +end +function fix_relative_phases(envfinal::CTMRGEnv, signs) + corners_fixed = _fix_gauge_corners(envfinal, signs) edges_fixed = map(eachcoordinate(envfinal, 1:4)) do (dir, r, c) Ef = if dir == NORTHWEST fix_gauge_north_edge((r, c), envfinal, signs) diff --git a/src/algorithms/ctmrg/gaugefix_split.jl b/src/algorithms/ctmrg/gaugefix_split.jl index 85792f555..5f4f00eb2 100644 --- a/src/algorithms/ctmrg/gaugefix_split.jl +++ b/src/algorithms/ctmrg/gaugefix_split.jl @@ -1,17 +1,17 @@ -"""List alternating edge factors clockwise, starting at a specified outer bond.""" -function _split_boundary_chain(env::SplitCTMRGEnv, (dir, r, c)) - dr, dc = dir == EAST ? (1, 0) : dir == WEST ? (-1, 0) : dir == NORTH ? (0, 1) : (0, -1) - count = iszero(dr) ? size(env, 3) : size(env, 2) +"""Represent a split boundary as a doubled MPS chain of alternating ket and bra factors, starting at an outer bond.""" +function _boundary_chain(env::SplitCTMRGEnv, coordinate) + dir = first(coordinate) firstlayer, secondlayer = env.ket_first[dir] ? (:ket, :bra) : (:bra, :ket) - return [edge(env, layer, dir, r + i * dr, c + i * dc) for i in 0:(count - 1) for layer in (firstlayer, secondlayer)] + return [edge(env, layer, I...) for I in _boundary_coordinates(env, coordinate) for layer in (firstlayer, secondlayer)] end """Determine outer and interlayer basis rotations relative to the previous split environment.""" function compute_relative_phases(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, ::ScramblingEnvGauge) gauges = map(eachcoordinate(final, 1:4)) do coordinate - current = _split_boundary_chain(final, coordinate) - old = _split_boundary_chain(previous, coordinate) + current = _boundary_chain(final, coordinate) + old = _boundary_chain(previous, coordinate) outer = _scrambling_gauge(current, old) + # The same method on the doubled chain, shifted to the interlayer cut. interlayer = _scrambling_gauge(circshift(current, -1), circshift(old, -1)) return outer, interlayer end @@ -20,12 +20,9 @@ end """Apply the outer gauges to corners and paired outer and interlayer gauges to edge factors.""" function fix_relative_phases(env::SplitCTMRGEnv, signs) - corner_fns = (fix_gauge_northwest_corner, fix_gauge_northeast_corner, fix_gauge_southeast_corner, fix_gauge_southwest_corner) - corners = dtmap(eachcoordinate(env, 1:4)) do (dir, r, c) - return corner_fns[dir]((r, c), env, signs.outer) - end + corners = _fix_gauge_corners(env, signs.outer) factors = dtmap(eachcoordinate(env, 1:4)) do (dir, r, c) - dr, dc = dir == EAST ? (1, 0) : dir == WEST ? (-1, 0) : dir == NORTH ? (0, 1) : (0, -1) + dr, dc = _boundary_step(dir) next = (dir, mod1(r + dr, size(env, 2)), mod1(c + dc, size(env, 3))) outer, interlayer = signs.outer[dir, r, c], signs.interlayer[dir, r, c] outer_next = signs.outer[next...] @@ -57,24 +54,12 @@ function fix_global_phases(env::SplitCTMRGEnv, corner_phases, edge_phases) ) end -"""Check split environment compatibility and determine relative gauges and global phases.""" -function compute_gauge_fix_gauge(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, alg::ScramblingEnvGauge) +"""Check that split environments have matching orientations and spaces for gauge alignment.""" +function _check_gauge_compatibility(final::SplitCTMRGEnv, previous::SplitCTMRGEnv) final.ket_first == previous.ket_first || throw(ArgumentError("Split-edge orientations must match.")) for field in (:corners, :edges_ket, :edges_bra) space.(getproperty(final, field)) == space.(getproperty(previous, field)) || throw(ArgumentError("Split environment spaces must match for gauge fixing.")) end - signs = compute_relative_phases(final, previous, alg) - corner_phases, edge_phases = compute_global_phases(fix_relative_phases(final, signs), previous) - return signs, corner_phases, edge_phases -end - -"""Apply relative gauge rotations and global phase corrections to a split environment.""" -function fix_phases(env::SplitCTMRGEnv, signs, corner_phases, edge_phases) - return fix_global_phases(fix_relative_phases(env, signs), corner_phases, edge_phases) -end - -"""Align a split environment with the outer and interlayer gauges of the previous environment.""" -function gauge_fix(final::SplitCTMRGEnv, previous::SplitCTMRGEnv, alg::ScramblingEnvGauge) - return fix_phases(final, compute_gauge_fix_gauge(final, previous, alg)...) + return nothing end From 87623cca9e93db502399a1aeb04602f6cf1a1533 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 12 Oct 2026 17:34:24 +0800 Subject: [PATCH 03/14] Fix formatting --- src/algorithms/ctmrg/split.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/algorithms/ctmrg/split.jl b/src/algorithms/ctmrg/split.jl index e2518fe77..31057e577 100644 --- a/src/algorithms/ctmrg/split.jl +++ b/src/algorithms/ctmrg/split.jl @@ -23,12 +23,14 @@ function check_input(::typeof(leading_boundary), network, env::SplitCTMRGEnv, al alg isa SplitCTMRG || throw(ArgumentError("SplitCTMRGEnv requires SplitCTMRG.")) eltype(network) <: PEPSSandwich || throw(ArgumentError("SplitCTMRG currently supports PEPS overlap networks only.")) size(env)[2:3] == size(network) || throw(DimensionMismatch("Network and environment unit cells must match.")) - env.ket_first == (true, true, false, false) || throw(ArgumentError( - "SplitCTMRG requires ket_first == (true, true, false, false); got $(env.ket_first). " * - "For a rotated standard environment, apply the same rotl90, rotr90, or rot180 transformation " * - "to both the environment and the PEPS state/network to restore this orientation before calling " * - "`leading_boundary`, then rotate the returned environment back." - )) + env.ket_first == (true, true, false, false) || throw( + ArgumentError( + "SplitCTMRG requires ket_first == (true, true, false, false); got $(env.ket_first). " * + "For a rotated standard environment, apply the same rotl90, rotr90, or rot180 transformation " * + "to both the environment and the PEPS state/network to restore this orientation before calling " * + "`leading_boundary`, then rotate the returned environment back." + ) + ) alg.maxiter > 0 || throw(ArgumentError("maxiter must be positive.")) return nothing end From 6495885f264cf556b3deec1fb4090e76ef65ae65 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 13 Oct 2026 12:14:57 +0800 Subject: [PATCH 04/14] Forward pass performance improvements --- docs/src/assets/pepskit.bib | 13 ++ docs/src/man/environments.md | 2 +- src/PEPSKit.jl | 3 +- .../ctmrg/enlarge_corner_split.jl | 24 +-- .../contractions/ctmrg/renormalize_split.jl | 10 +- .../contractions/ctmrg/split_env.jl | 14 +- src/algorithms/ctmrg/ctmrg.jl | 4 +- src/algorithms/ctmrg/projectors/split.jl | 169 ++++++++++++------ src/algorithms/ctmrg/split.jl | 6 +- src/utility/lq.jl | 39 ++++ src/utility/qr.jl | 10 +- test/testsuite/TestSuite.jl | 6 + test/testsuite/utility/qr_wrapper.jl | 48 +++++ test/utility/qr_wrapper.jl | 11 ++ 14 files changed, 272 insertions(+), 87 deletions(-) create mode 100644 src/utility/lq.jl create mode 100644 test/testsuite/utility/qr_wrapper.jl create mode 100644 test/utility/qr_wrapper.jl diff --git a/docs/src/assets/pepskit.bib b/docs/src/assets/pepskit.bib index 7eac303f7..0be60cf92 100644 --- a/docs/src/assets/pepskit.bib +++ b/docs/src/assets/pepskit.bib @@ -11,6 +11,19 @@ @article{francuz_stable_2025 doi = {10.1103/PhysRevResearch.7.013237} } +@article{naumann_variationally_2025, + title = {Variationally optimizing infinite projected entangled-pair states at large bond dimensions: A split corner transfer matrix renormalization group approach}, + author = {Naumann, Jan and Weerda, Erik L. and Eisert, Jens and Rizzi, Matteo and Schmoll, Philipp}, + journal = {Phys. Rev. B}, + volume = {111}, + pages = {235116}, + year = {2025}, + month = {Jun}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.111.235116}, + url = {https://doi.org/10.1103/PhysRevB.111.235116} +} + @article{evenbly_gauge_2018, title = {Gauge fixing, canonical forms, and optimal truncations in tensor networks with closed loops}, author = {Evenbly, Glen}, diff --git a/docs/src/man/environments.md b/docs/src/man/environments.md index 3118c338c..498d3631f 100644 --- a/docs/src/man/environments.md +++ b/docs/src/man/environments.md @@ -25,7 +25,7 @@ When omitted, `interlayer_trunc` inherits `trunc`; the default `FixedSpaceTrunca This version supports PEPS overlap networks, full projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref). For comparisons with variPEPS, select its `Projector_Method.FULL` option; its default `FISHMAN` construction reweights the environment halves and uses different projectors. -Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. +Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref PEPSKit.ScramblingEnvGauge) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. The ordinary and split environments use the same boundary-MPS gauge method, with each split edge contributing two MPS sites instead of one. Both sets of bond spaces are held fixed in the differentiated iteration, including when the forward run uses other truncation strategies. Use a complex environment for fixed-point gradients, as for conventional asymmetric CTMRG. diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index e4e0220c6..157b18fc8 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -60,6 +60,7 @@ include("utility/twistdual.jl") include("utility/eigh.jl") include("utility/svd.jl") include("utility/qr.jl") +include("utility/lq.jl") include("utility/rotations.jl") include("utility/hook_pullback.jl") include("utility/autoopt.jl") @@ -176,7 +177,7 @@ include("algorithms/select_algorithm.jl") using .Defaults: set_scheduler! export set_scheduler! -export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint +export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint, LQAdjoint export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG export SplitCTMRGEnv, SplitCTMRG export initialize_ctmrg_environment, diff --git a/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl b/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl index ba164dffd..d58788bee 100644 --- a/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl +++ b/src/algorithms/contractions/ctmrg/enlarge_corner_split.jl @@ -1,30 +1,30 @@ # Corner enlargements # ------------------- -"""Project the ket cut of a northwest enlarged corner before contracting a full strip.""" -function enlarge_northwest_corner((Wk, Wb)::Tuple, C, N, A::PEPSSandwich, G) +"""Build the northwest enlarged corner for the environment of the teal projectors.""" +function enlarge_northwest_corner((Wk, Wb)::Tuple, C, N, A::PEPSSandwich, green) return @autoopt @tensor Q[χout Dsb; χne Dek Deb] := - G[χout; χsw Dsk] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C[χnw; χnnw] * N[χnnw Dnk Dnb; χne] * + green[χout; χsw Dsk] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * C[χnw; χnnw] * N[χnnw Dnk Dnb; χne] * ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) end -"""Project the ket cut of a northeast enlarged corner before contracting a full strip.""" -function enlarge_northeast_corner(N, C, (Ek, Eb)::Tuple, A::PEPSSandwich, G) +"""Build the northeast enlarged corner for the environment of the teal projectors.""" +function enlarge_northeast_corner(N, C, (Ek, Eb)::Tuple, A::PEPSSandwich, green) return @autoopt @tensor Q[χnw Dwk Dwb; χout Dsb] := - N[χnw Dnk Dnb; χnne] * C[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * G[χse Dsk; χout] * + N[χnw Dnk Dnb; χnne] * C[χnne; χne] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * green[χse Dsk; χout] * ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) end -"""Project the ket cut of a southwest enlarged corner before contracting a full strip.""" -function enlarge_southwest_corner(S, C, (Wk, Wb)::Tuple, A::PEPSSandwich, G) +"""Build the southwest enlarged corner for the environment of the teal projectors.""" +function enlarge_southwest_corner(S, C, (Wk, Wb)::Tuple, A::PEPSSandwich, green) return @autoopt @tensor Q[χse Dek Deb; χout Dnb] := - S[χse Dsk Dsb; χssw] * C[χssw; χsw] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * G[χnw Dnk; χout] * + S[χse Dsk Dsb; χssw] * C[χssw; χsw] * Wb[χsw Dwb; χiw] * Wk[χiw Dwk; χnw] * green[χnw Dnk; χout] * ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) end -"""Project the ket cut of a southeast enlarged corner before contracting a full strip.""" -function enlarge_southeast_corner((Ek, Eb)::Tuple, C, S, A::PEPSSandwich, G) +"""Build the southeast enlarged corner for the environment of the teal projectors.""" +function enlarge_southeast_corner((Ek, Eb)::Tuple, C, S, A::PEPSSandwich, green) return @autoopt @tensor Q[χout Dnb; χsw Dwk Dwb] := - G[χout; χne Dnk] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * C[χse; χsse] * S[χsse Dsk Dsb; χsw] * + green[χout; χne Dnk] * Ek[χne Dek; χie] * Eb[χie Deb; χse] * C[χse; χsse] * S[χsse Dsk Dsb; χsw] * ket(A)[d; Dnk Dek Dsk Dwk] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) end diff --git a/src/algorithms/contractions/ctmrg/renormalize_split.jl b/src/algorithms/contractions/ctmrg/renormalize_split.jl index e8081d842..9916978c0 100644 --- a/src/algorithms/contractions/ctmrg/renormalize_split.jl +++ b/src/algorithms/contractions/ctmrg/renormalize_split.jl @@ -13,18 +13,18 @@ function renormalize_southwest_corner(C::CTMRGCornerTensor, S::CTMRG_PEPS_EdgeTe return @tensor χout[χe; χw] := tmp[χe Dnb; χg] * teal[χg Dnb; χw] end -"""Absorb the ket layer and apply both outer and interlayer projector stages to a west edge.""" +"""Absorb the ket layer and apply the split projectors to a west-edge factor.""" function renormalize_west_edge(K::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor) return @autoopt @tensor E[χinew Dek; χnnew] := K[χiold Dwk; χnold] * A[d; Dnk Dek Dsk Dwk] * green[χnold Dnk; χng] * yellow[χiy; χiold Dsk] * teal[χng Dnb; χnnew] * red[χinew; χiy d Dnb] end -"""Absorb the bra layer and apply the matching projectors to the other west-edge factor.""" +"""Absorb the bra layer and apply the split projectors to the other west-edge factor.""" function renormalize_west_edge(B::CTMRGEdgeTensor{<:Any, <:Any, 2}, green, teal, yellow, red, A::PEPSTensor, ::Val{:bra}) @autoopt @tensor BY[χsg Dwb; χiy] := green[χsg; χsold Dsk] * B[χsold Dwb; χiold] * yellow[χiold Dsk; χiy] - @tensor BYA[χsg Dnb Deb Dsb; χiy d] := BY[χsg Dwb; χiy] * conj(A[d; Dnb Deb Dsb Dwb]) - @tensor BYAT[χsnew Dnb Deb; χiy d] := teal[χsnew; χsg Dsb] * BYA[χsg Dnb Deb Dsb; χiy d] - return @tensor E[χsnew Deb; χinew] := BYAT[χsnew Dnb Deb; χiy d] * red[χiy d Dnb; χinew] + @tensor BYT[χsnew Dwb Dsb; χiy] := teal[χsnew; χsg Dsb] * BY[χsg Dwb; χiy] + @tensor BYTA[χsnew Dnb Deb; χiy d] := BYT[χsnew Dwb Dsb; χiy] * conj(A[d; Dnb Deb Dsb Dwb]) + return @tensor E[χsnew Deb; χinew] := BYTA[χsnew Dnb Deb; χiy d] * red[χiy d Dnb; χinew] end diff --git a/src/algorithms/contractions/ctmrg/split_env.jl b/src/algorithms/contractions/ctmrg/split_env.jl index c05d031f4..e4e2c4ede 100644 --- a/src/algorithms/contractions/ctmrg/split_env.jl +++ b/src/algorithms/contractions/ctmrg/split_env.jl @@ -54,27 +54,27 @@ function _split_yellow_environment((r, c), network, env::SplitCTMRGEnv) return top, bottom end -"""Build the two halves of the environment for the red projectors (Fig. 9).""" -function _split_red_environment((r, c), network, env::SplitCTMRGEnv, yellow_west, yellow_east) +"""Build the quarter-environments for the red projectors before multiplying their factors (Fig. 9).""" +function _split_red_quarters((r, c), network, env::SplitCTMRGEnv, yellow_west, yellow_east) l = _split_local_environment((r, c), network, env) h = _split_local_environment((r, c + 1), network, env) Ywl, Ywr = yellow_west[1][r], yellow_west[2][r] Yel, Yer = yellow_east[1][r], yellow_east[2][r] C1, N, Wk, A = l.C1, l.N, l.Wk, l.A - NW = @autoopt @tensor Q[χout d Dnb; χne Dek] := + @autoopt @tensor NW[χout d Dnb; χne Dek] := Ywr[χout; χiw Dsk] * Wk[χiw Dwk; χnw] * C1[χnw; χnnw] * N[χnnw Dnk Dnb; χne] * ket(A)[d; Dnk Dek Dsk Dwk] C2, N, Ek, A = h.C2, h.N, h.Ek, h.A - NE = @autoopt @tensor Q[χnw Dwk; χout d Dnb] := + @autoopt @tensor NE[χnw Dwk; χout d Dnb] := N[χnw Dnk Dnb; χnne] * C2[χnne; χne] * Ek[χne Dek; χie] * Yel[χie Dsk; χout] * ket(A)[d; Dnk Dek Dsk Dwk] C4, S, Wb, A = l.C4, l.S, l.Wb, l.A - SW = @autoopt @tensor Q[χse Deb; χout d Dnb] := + @autoopt @tensor SW[χse Deb; χout d Dnb] := S[χse Dsk Dsb; χssw] * C4[χssw; χsw] * Wb[χsw Dwb; χiw] * Ywl[χiw Dsk; χout] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) C3, S, Eb, A = h.C3, h.S, h.Eb, h.A - SE = @autoopt @tensor Q[χout d Dnb; χsw Dwb] := + @autoopt @tensor SE[χout d Dnb; χsw Dwb] := Yer[χout; χie Dsk] * Eb[χie Deb; χse] * C3[χse; χsse] * S[χsse Dsk Dsb; χsw] * conj(bra(A)[d; Dnb Deb Dsb Dwb]) - return NW * NE, SE * SW + return (; NW, NE, SE, SW) end diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index a289d5096..abda74d78 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -80,7 +80,7 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st - `:SimultaneousCTMRG` : Simultaneous expansion and renormalization of all sides. - `:SequentialCTMRG` : Sequential application of left moves and rotations. - `:C4vCTMRG` : CTMRG assuming C₄ᵥ-symmetric PEPS and environment. - - `:SplitCTMRG` : Sequential ket and bra absorption using [`SplitCTMRGEnv`](@ref) and full split projectors. + - `:SplitCTMRG` : Sequential ket and bra absorption using [`SplitCTMRGEnv`](@ref) and the split projectors of [`SplitProjector`](@ref). ### Projector algorithm @@ -94,7 +94,7 @@ supplied via the keyword arguments or directly as an [`CTMRGAlgorithm`](@ref) st * `projector_alg::Symbol=:$(Defaults.projector_alg)` : Variant of the projector algorithm. See also [`ProjectorAlgorithm`](@ref). - `:HalfInfiniteProjector` : Projection via SVDs of half-infinite (two enlarged corners) CTMRG environments. - `:FullInfiniteProjector` : Projection via SVDs of full-infinite (all four enlarged corners) CTMRG environments. - - `:SplitProjector` : Projection via SVDs of the split environments, works only for [`SplitCTMRG`](@ref) and is its default. + - `:SplitProjector` : Split projector pairs constructed via SVDs of their environments, works only for [`SplitCTMRG`](@ref) and is its default. - `:C4vEighProjector` : Projection via `eigh` of the Hermitian enlarged corner, works only for [`C4vCTMRG`](@ref). - `:C4vQRProjector` : Projection via QR decomposition of the lower-rank column-enlarged corner, works only for [`C4vCTMRG`](@ref). * `decomposition_alg::Union{NamedTuple,<:SVDAdjoint,<:EighAdjoint,<:QRAdjoint}` : Tensor diff --git a/src/algorithms/ctmrg/projectors/split.jl b/src/algorithms/ctmrg/projectors/split.jl index 201145111..664ecb573 100644 --- a/src/algorithms/ctmrg/projectors/split.jl +++ b/src/algorithms/ctmrg/projectors/split.jl @@ -1,9 +1,9 @@ """ $(TYPEDEF) -Projector algorithm for [`SplitCTMRG`](@ref), using the green, teal, yellow, and red environments of Naumann et al., Phys. Rev. B **111**, 235116 (2025). +Projector algorithm for [`SplitCTMRG`](@ref), constructing the green (Fig. 3), teal (Fig. 7), yellow (Fig. 8), and red (Fig. 9) projector pairs of [Naumann et al.](@cite naumann_variationally_2025). Each projector pair is constructed from a truncated SVD of the product of its environment halves. -The yellow and red halves are QR-compressed before multiplication. +QR compression is applied to the environment halves for the yellow projectors and, when the truncation permits, the quarter-environments for the red projectors before multiplication. ## Fields @@ -14,7 +14,7 @@ $(TYPEDFIELDS) SplitProjector(; kwargs...) Construct a split projector algorithm with `decomposition_alg`, `trunc`, and `verbosity`, following [`FullInfiniteProjector`](@ref). -`trunc` controls the outer cuts; [`SplitCTMRG`](@ref) supplies `interlayer_trunc` for the yellow and red cuts. +`trunc` controls the green and teal projectors; [`SplitCTMRG`](@ref) supplies `interlayer_trunc` for the yellow and red projectors. """ struct SplitProjector{S <: SVDAdjoint, T} <: ProjectorAlgorithm decomposition_alg::S @@ -27,7 +27,7 @@ end PROJECTOR_SYMBOLS[:SplitProjector] = SplitProjector -"""Compute a split projector pair, optionally QR-compressing the yellow or red environment halves.""" +"""Compute a split projector pair, optionally QR-compressing the environment halves for the yellow or red projectors.""" function compute_projector( halves, alg::SplitProjector; trunc::TruncationStrategy = alg.trunc, compress::Bool = false, @@ -37,10 +37,8 @@ function compute_projector( Q = Q / norm(Q) Q_next = Q_next / norm(Q_next) if compress - qr_alg = QRAdjoint() - _, Q = left_orth(Q, qr_alg) - _, Q_next = left_orth(Q_next', qr_alg) - Q_next = Q_next' + _, Q = left_orth(Q, QRAdjoint()) + Q_next, _ = right_orth(Q_next, LQAdjoint()) end # All four colors use an SVD of the product of their prepared halves. @@ -60,73 +58,138 @@ function compute_projector( return (P_left, P_right), (; U, S, V, truncation_error) end -"""Resolve fixed-space truncation at an outer or interlayer cut.""" +"""Resolve fixed-space truncation for the green and teal projectors or the yellow and red projectors.""" function _split_truncation(trunc::TruncationStrategy, V::ElementarySpace) return trunc isa FixedSpaceTruncation ? truncspace(isdual(V) ? flip(V) : V) : trunc end -"""Compute green or yellow projector pairs before constructing the teal or red environments.""" -function _split_initial_projectors( - col::Int, network, env::SplitCTMRGEnv, alg::SplitProjector, trunc::TruncationStrategy, - ::Val{interlayer}, ::Val{east} - ) where {interlayer, east} - halves = dtmap(axes(env, 2)) do r +"""Construct the half-environments for the green or yellow projectors for each row of a directional move.""" +function _split_initial_halves( + col::Int, network, env::SplitCTMRGEnv, ::Val{interlayer} + ) where {interlayer} + return dtmap(axes(env, 2)) do r return interlayer ? _split_yellow_environment((r, col), network, env) : _split_green_environment((r, col), network, env) end +end + +"""Compute a green projector pair at one cut from prepared half-environments.""" +function _split_green_projector( + (r, col), halves, env::SplitCTMRGEnv, alg::SplitProjector, ::Val{E} = Val(false); + trunc::TruncationStrategy = alg.trunc, + ) where {E} + top = halves[_prev(r, size(env, 2))][1] + bottom = halves[r][2] + if E + top = permute(top, ((1, 4), (3, 2))) + bottom = permute(bottom, ((1, 4), (3, 2))) + end + K = edge(env, :ket, E ? EAST : WEST, r, E ? col + 1 : col) + V = E ? codomain(K)[1] : domain(K)[1] + Q, Q_next = E ? (top, bottom) : (bottom, top) + return compute_projector((Q, Q_next), alg; trunc = _split_truncation(trunc, V)) +end + +"""Compute a yellow projector pair at one cut using QR/LQ compression of prepared half-environments.""" +function _split_yellow_projector( + (r, col), halves, env::SplitCTMRGEnv, alg::SplitProjector, ::Val{E} = Val(false); + trunc::TruncationStrategy = alg.trunc, + ) where {E} + top, bottom = halves[r] + if E + top = permute(top, ((1, 4, 5), (3, 2))) + bottom = permute(bottom, ((1, 5), (4, 2, 3))) + end + K = edge(env, :ket, E ? EAST : WEST, r, E ? col + 1 : col) + V = E ? domain(K)[1] : codomain(K)[1] + Q, Q_next = E ? (top, bottom) : (bottom, top) + return compute_projector((Q, Q_next), alg; trunc = _split_truncation(trunc, V), compress = true) +end + +"""Collect green or yellow projector pairs for one side of a column from prepared half-environments.""" +function _split_projectors_from_halves( + col::Int, halves, env::SplitCTMRGEnv, alg::SplitProjector, trunc::TruncationStrategy, + ::Val{interlayer}, east::Val, + ) where {interlayer} results = dtmap(axes(env, 2)) do r - if interlayer - top, bottom = halves[r] - if east - top = permute(top, ((1, 4, 5), (3, 2))) - bottom = permute(bottom, ((1, 5), (4, 2, 3))) - end - Q, Q_next = east ? (top, bottom) : (bottom, top) - K = edge(env, :ket, east ? EAST : WEST, r, east ? col + 1 : col) - V = east ? domain(K)[1] : codomain(K)[1] - else - top = halves[_prev(r, size(env, 2))][1] - bottom = halves[r][2] - if east - top = permute(top, ((1, 4), (3, 2))) - bottom = permute(bottom, ((1, 4), (3, 2))) - end - Q, Q_next = east ? (top, bottom) : (bottom, top) - K = edge(env, :ket, east ? EAST : WEST, r, east ? col + 1 : col) - V = east ? codomain(K)[1] : domain(K)[1] - end - return compute_projector((Q, Q_next), alg; trunc = _split_truncation(trunc, V), compress = interlayer) + return interlayer ? _split_yellow_projector((r, col), halves, env, alg, east; trunc) : + _split_green_projector((r, col), halves, env, alg, east; trunc) end projectors, info = _split_proj_and_info(results) return projectors, info.truncation_error end -"""Construct the green, teal, yellow, and red projector pairs for a split-CTMRG left move.""" +"""Compute a teal projector pair from adjacent-row green-projected halves using the outer truncation strategy.""" +function _split_teal_projector( + (r, col), halves, env::SplitCTMRGEnv, alg::SplitProjector; + trunc::TruncationStrategy = alg.trunc, + ) + top = halves[_prev(r, size(env, 2))][1] + bottom = halves[r][2] + V = domain(edge(env, :ket, WEST, r, col))[1] + trunc = _split_truncation(trunc, V) + return compute_projector((bottom, top), alg; trunc) +end + +"""Check whether the reduced environment spectrum for the red projectors can supply the requested largest singular values.""" +function _split_reduce_red(quarters::NamedTuple, trunc::TruncationStrategy) + trunc isa TruncationSpace && trunc.by === abs && trunc.rev || return false + (; NW, NE, SE, SW) = quarters + # Respect both the reduced quarter ranks and the original half-QR dimensions. + return all(blocksectors(trunc.space)) do c + return blockdim(trunc.space, c) <= min( + blockdim(codomain(SE), c), blockdim(domain(SE), c), + blockdim(domain(NE), c), blockdim(codomain(NE), c), + blockdim(domain(SW), c), blockdim(codomain(NW), c), + ) + end +end + +"""Compute the red projector pair using quarter QR factors when its truncation fits the reduced spectrum.""" +function _split_red_projector( + (r, col), quarters::NamedTuple, env::SplitCTMRGEnv, alg::SplitProjector; + trunc::TruncationStrategy = alg.trunc, + ) + V = codomain(edge(env, :ket, WEST, r, col))[1] + trunc = _split_truncation(trunc, V) + (; NW, NE, SE, SW) = quarters + reduce = ignore_derivatives() do + _split_reduce_red(quarters, trunc) + end + if reduce + # SE = Qb * Rb and NE = Lt * Qt with isometric Qb and Qt. + # Their factors retain the norms and nonzero spectrum of (SE * SW) * (NW * NE). + _, Rb = left_orth(SE, QRAdjoint()) + Lt, _ = right_orth(NE, LQAdjoint()) + return compute_projector((Rb * SW, NW * Lt), alg; trunc) + end + # Preserve zero-spectrum dimensions and custom truncation ordering on the full halves. + return compute_projector((SE * SW, NW * NE), alg; trunc, compress = true) +end + +"""Construct the split projector pairs for a split-CTMRG left move.""" function sequential_projectors( col::Int, network, env::SplitCTMRGEnv, alg::SplitProjector, interlayer_trunc::TruncationStrategy ) - # Initial projectors on both sides of the two-site refinement environments. - green, eg = _split_initial_projectors(col, network, env, alg, alg.trunc, Val(false), Val(false)) - green_east, ege = _split_initial_projectors(col + 1, network, env, alg, alg.trunc, Val(false), Val(true)) - yellow, ey = _split_initial_projectors(col, network, env, alg, interlayer_trunc, Val(true), Val(false)) - yellow_east, eye = _split_initial_projectors(col + 1, network, env, alg, interlayer_trunc, Val(true), Val(true)) + # Prepare each side once; a one-column cell shares the same periodic halves. + single_column = size(env, 3) == 1 + green_halves = _split_initial_halves(col, network, env, Val(false)) + green, eg = _split_projectors_from_halves(col, green_halves, env, alg, alg.trunc, Val(false), Val(false)) + green_east_halves = single_column ? green_halves : _split_initial_halves(col + 1, network, env, Val(false)) + green_east, ege = _split_projectors_from_halves(col + 1, green_east_halves, env, alg, alg.trunc, Val(false), Val(true)) + yellow_halves = _split_initial_halves(col, network, env, Val(true)) + yellow, ey = _split_projectors_from_halves(col, yellow_halves, env, alg, interlayer_trunc, Val(true), Val(false)) + yellow_east_halves = single_column ? yellow_halves : _split_initial_halves(col + 1, network, env, Val(true)) + yellow_east, eye = _split_projectors_from_halves(col + 1, yellow_east_halves, env, alg, interlayer_trunc, Val(true), Val(true)) # Teal refines the outer cut using environments with green projectors inserted. teal_halves = dtmap(r -> _split_teal_environment((r, col), network, env, green, green_east), axes(env, 2)) - teals = dtmap(axes(env, 2)) do r - top = teal_halves[_prev(r, size(env, 2))][1] - bottom = teal_halves[r][2] - V = domain(edge(env, :ket, WEST, r, col))[1] - trunc = _split_truncation(alg.trunc, V) - return compute_projector((bottom, top), alg; trunc) - end + teals = dtmap(r -> _split_teal_projector((r, col), teal_halves, env, alg), axes(env, 2)) # Red refines the interlayer cut using environments with yellow projectors inserted. reds = dtmap(axes(env, 2)) do r - top, bottom = _split_red_environment((r, col), network, env, yellow, yellow_east) - V = codomain(edge(env, :ket, WEST, r, col))[1] - trunc = _split_truncation(interlayer_trunc, V) - return compute_projector((bottom, top), alg; trunc, compress = true) + quarters = _split_red_quarters((r, col), network, env, yellow, yellow_east) + return _split_red_projector((r, col), quarters, env, alg; trunc = interlayer_trunc) end teal, teal_info = _split_proj_and_info(teals) red, red_info = _split_proj_and_info(reds) diff --git a/src/algorithms/ctmrg/split.jl b/src/algorithms/ctmrg/split.jl index 31057e577..582237243 100644 --- a/src/algorithms/ctmrg/split.jl +++ b/src/algorithms/ctmrg/split.jl @@ -1,11 +1,11 @@ """ SplitCTMRG(; kwargs...) -Sequential split-CTMRG introduced in Naumann et al., Phys. Rev. B **111**, 235116 (2025), Figs. 3 and 7–9. +Sequential split-CTMRG introduced in [Naumann et al.](@cite naumann_variationally_2025). Uses [`SplitCTMRGEnv`](@ref), with independently truncated outer and interlayer environment bonds. The common iteration, `trunc`, and `decomposition_alg` keywords follow [`SequentialCTMRG`](@ref). -`interlayer_trunc` controls the interlayer projectors and defaults to `trunc`; specify it separately to override the outer truncation strategy. -Uses [`SplitProjector`](@ref) to construct the green, teal, yellow, and red projector pairs. +`trunc` controls the green and teal projectors; `interlayer_trunc` controls the yellow and red projectors and defaults to `trunc`. +Uses [`SplitProjector`](@ref) to construct the split projector pairs. The fixed-point gradient aligns the outer and interlayer gauges and holds both environment spaces fixed. """ struct SplitCTMRG{P <: SplitProjector, I <: TruncationStrategy} <: CTMRGAlgorithm diff --git a/src/utility/lq.jl b/src/utility/lq.jl new file mode 100644 index 000000000..37396c194 --- /dev/null +++ b/src/utility/lq.jl @@ -0,0 +1,39 @@ +""" +$(TYPEDEF) + +Wrapper for an LQ decomposition algorithm with a defined reverse rule, used by `right_orth(!)`. +See [`QRAdjoint`](@ref) for details of the constructor options and defaults. +""" +struct LQAdjoint{F, R} + fwd_alg::F + rrule_alg::R +end + +LQAdjoint(; kwargs...) = LQAdjoint(_orthogonalization_algorithms(; kwargs...)...) + +"""Compute an LQ decomposition using the forward algorithm and reverse rule configured by `LQAdjoint`.""" +MatrixAlgebraKit.right_orth(t, alg::LQAdjoint) = right_orth!(copy(t), alg) +"""Factorize the input in place using the LQ forward algorithm configured by `LQAdjoint`.""" +MatrixAlgebraKit.right_orth!(t, alg::LQAdjoint) = right_orth!(t, alg.fwd_alg) + +# right_orth! rrule wrapping MatrixAlgebraKit's lq_pullback! +function ChainRulesCore.rrule( + ::typeof(right_orth!), + t::AbstractTensorMap, + alg::LQAdjoint{F, R}, + ) where {F <: MatrixAlgebraKit.Algorithm, R <: FullPullback} + LQ = right_orth(t, alg) + gtol = _get_pullback_gauge_tol(alg.rrule_alg.verbosity) + + """Propagate LQ cotangents with the configured gauge-sensitivity tolerance.""" + function right_orth!_pullback(ΔLQ) + Δt = zeros(scalartype(t), space(t)) + MatrixAlgebraKit.lq_pullback!(Δt, t, LQ, unthunk.(ΔLQ); gauge_atol = gtol(ΔLQ)) + return NoTangent(), Δt, NoTangent() + end + function right_orth!_pullback(::Tuple{ZeroTangent, ZeroTangent}) + return NoTangent(), ZeroTangent(), NoTangent() + end + + return LQ, right_orth!_pullback +end diff --git a/src/utility/qr.jl b/src/utility/qr.jl index a512821c8..194a5f5f6 100644 --- a/src/utility/qr.jl +++ b/src/utility/qr.jl @@ -19,7 +19,8 @@ Construct a `QRAdjoint` algorithm struct based on the following keyword argument The available algorithms are provided through MatrixAlgebraKit and include: - `:DefaultAlgorithm` : MatrixAlgebraKit's [default QR algorithm](@extref MatrixAlgebraKit.DefaultAlgorithm) for a given matrix type. - `:Householder` : MatrixAlgebraKit's [`Householder`](@extref MatrixAlgebraKit.Householder) -* `rrule_alg::Union{Algorithm,NamedTuple}=(; alg::Symbol=$(Defaults.qr_rrule_alg))`: Reverse-rule algorithm for differentiating the eigenvalue decomposition. Can be supplied by an `Algorithm` instance directly or as a `NamedTuple` where `alg` is one of the following: +* `rrule_alg::Union{Algorithm,NamedTuple}=(; alg::Symbol=$(Defaults.qr_rrule_alg))`: Reverse-rule algorithm for differentiating the QR decomposition. + Can be supplied by an `Algorithm` instance directly or as a `NamedTuple` where `alg` is one of the following: - `:FullPullback` : MatrixAlgebraKit's [`qr_pullback!`](@extref MatrixAlgebraKit.qr_pullback!) !!! note @@ -39,7 +40,10 @@ const QR_RRULE_SYMBOLS = IdDict{Symbol, Type{<:Any}}( :FullPullback => FullPullback ) -function QRAdjoint(; fwd_alg = (;), rrule_alg = (;)) +QRAdjoint(; kwargs...) = QRAdjoint(_orthogonalization_algorithms(; kwargs...)...) + +"""Parse the shared forward and reverse algorithm options for QR and LQ decompositions.""" +function _orthogonalization_algorithms(; fwd_alg = (;), rrule_alg = (;)) # parse forward algorithm fwd_algorithm = if fwd_alg isa NamedTuple fwd_kwargs = (; @@ -72,7 +76,7 @@ function QRAdjoint(; fwd_alg = (;), rrule_alg = (;)) rrule_alg end - return QRAdjoint(fwd_algorithm, rrule_algorithm) + return fwd_algorithm, rrule_algorithm end """ diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index fa550601a..cd210e393 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -273,6 +273,12 @@ module UtilityEighWrapper end using .UtilityEighWrapper +module UtilityQRWrapper + include("utility/qr_wrapper.jl") + export utility_qr_wrapper +end +using .UtilityQRWrapper + module UtilitySVDWrapper include("utility/svd_wrapper.jl") export utility_svd_wrapper diff --git a/test/testsuite/utility/qr_wrapper.jl b/test/testsuite/utility/qr_wrapper.jl new file mode 100644 index 000000000..bf41a12a3 --- /dev/null +++ b/test/testsuite/utility/qr_wrapper.jl @@ -0,0 +1,48 @@ +using Test, Random, LinearAlgebra +using TensorKit, MatrixAlgebraKit, PEPSKit, Zygote, Adapt +using ChainRulesCore: rrule, NoTangent, ZeroTangent + +"""Check LQ factors and pullbacks against adjointed QR, including mixed dual fermionic legs.""" +function utility_qr_wrapper(AT) + Random.seed!(73194) + qr_alg, lq_alg = QRAdjoint(), LQAdjoint() + spaces = (ComplexSpace(2), Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionNumber](-1 => 1, 0 => 1, 1 => 1)) + return @testset "QR/LQ wrapper ($AT)" begin + for V in spaces, (nout, nin) in ((2, 3), (3, 2)) + output = prod(ntuple(i -> isodd(i) ? V : V', nout)) + input = prod(ntuple(i -> isodd(i) ? V' : V, nin)) + A = adapt(AT, randn(ComplexF64, output ← input)) + reference = X -> begin + U, R = left_orth(X', qr_alg) + return R', U' + end + factorize = X -> right_orth(X, lq_alg) + L, Q = factorize(A) + Lref, Qref = reference(A) + @test L ≈ Lref + @test Q ≈ Qref + @test L * Q ≈ A + @test right_orth!(copy(A), lq_alg) == (L, Q) + + H = adapt(AT, randn(ComplexF64, output ← output)) + H = H + H' + K = adapt(AT, randn(ComplexF64, input ← input)) + K = K + K' + # Both factors contribute to this gauge-invariant loss. + loss = decomposition -> X -> begin + l, q = decomposition(X) + return real(dot(H, l * l')) + real(dot(K, q' * q)) + end + gradient = only(Zygote.gradient(loss(factorize), A)) + @test gradient ≈ only(Zygote.gradient(loss(reference), A)) + direction = adapt(AT, randn(ComplexF64, space(A))) + direction /= norm(direction) + epsilon = 1.0e-5 + finite = (loss(factorize)(A + epsilon * direction) - loss(factorize)(A - epsilon * direction)) / (2 * epsilon) + @test real(dot(gradient, direction)) ≈ finite atol = 1.0e-7 rtol = 1.0e-7 + + _, back = rrule(right_orth!, copy(A), lq_alg) + @test back((ZeroTangent(), ZeroTangent())) == (NoTangent(), ZeroTangent(), NoTangent()) + end + end +end diff --git a/test/utility/qr_wrapper.jl b/test/utility/qr_wrapper.jl new file mode 100644 index 000000000..f39bea571 --- /dev/null +++ b/test/utility/qr_wrapper.jl @@ -0,0 +1,11 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.utility_qr_wrapper(Vector) +end From 359ef23fe8b7ab06f937f16d414b4ab19227ddc1 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 13 Oct 2026 15:56:35 +0800 Subject: [PATCH 05/14] FP gradient improvements --- .../optimization/implicit_differentiation.jl | 79 +++++++++++++++---- src/environments/split_ctmrg_environments.jl | 53 ++++++++++--- .../gradients/split_ctmrg_gradients.jl | 24 +++++- 3 files changed, 128 insertions(+), 28 deletions(-) diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl index 395de31bc..3c79141fb 100644 --- a/src/algorithms/optimization/implicit_differentiation.jl +++ b/src/algorithms/optimization/implicit_differentiation.jl @@ -216,6 +216,65 @@ function _set_fixed_truncation(alg::SplitCTMRG) return @set alg_fixed.interlayer_trunc = FixedSpaceTruncation() end +"""Prepare the joint pullback of a fixed-space, gauge-aligned CTMRG iteration.""" +function _fixedpoint_vjp(config::RuleConfig, state, env, alg::CTMRGAlgorithm) + alg_fixed = _set_fixed_truncation(alg) + env_conv, info = ctmrg_iteration(InfiniteSquareNetwork(state), env, alg_fixed) + signs, corner_phases, edge_phases = compute_gauge_fix_gauge(env_conv, env, _scrambling_env_gauge(alg)) + function gauge_fixed_iteration(A, x) + return fix_phases( + ctmrg_iteration(InfiniteSquareNetwork(A), x, alg_fixed)[1], + signs, corner_phases, edge_phases, + ) + end + _, env_vjp = rrule_via_ad(config, gauge_fixed_iteration, state, env) + return info, env_vjp +end + +"""Reuse the tracked split sweep for gauge selection and compose its phase pullback.""" +function _fixedpoint_vjp(config::RuleConfig, state, env::SplitCTMRGEnv, alg::SplitCTMRG) + alg_fixed = _set_fixed_truncation(alg) + iteration(A, x) = ctmrg_iteration(InfiniteSquareNetwork(A), x, alg_fixed) + (env_conv, info), sweep_vjp = rrule_via_ad(config, iteration, state, env) + signs, corner_phases, edge_phases = compute_gauge_fix_gauge(env_conv, env, _scrambling_env_gauge(alg)) + align(x) = fix_phases(x, signs, corner_phases, edge_phases) + _, phase_vjp = rrule_via_ad(config, align, env_conv) + env_vjp(x) = sweep_vjp((phase_vjp(x)[2], ZeroTangent())) + return info, env_vjp +end + +"""Separate the environment and state actions of a joint iteration pullback.""" +function _fixedpoint_vjps(env_vjp, state, env) + ∂f∂x(x)::typeof(env) = env_vjp(x)[3] + ∂f∂A(x)::typeof(state) = env_vjp(x)[2] + return ∂f∂x, ∂f∂A +end + +"""Reuse the state cotangent for an exactly matching split seed, owning all cached storage.""" +function _fixedpoint_vjps(env_vjp, state, env::SplitCTMRGEnv) + cache = Ref{Union{Nothing, Tuple{typeof(env), typeof(state)}}}(nothing) + function ∂f∂x(x)::typeof(env) + cache[] = nothing + seed = deepcopy(x) + result = env_vjp(x) + cache[] = (seed, deepcopy(result[2])) + # Krylov solvers mutate the environment result, so retain only the state part. + return _project_split_tangent(result[3], env) + end + function ∂f∂A(x)::typeof(state) + cached = cache[] + if !isnothing(cached) + seed, gradient = cached + if x.ket_first == seed.ket_first && x.corners == seed.corners && + x.edges_ket == seed.edges_ket && x.edges_bra == seed.edges_bra + return deepcopy(gradient) + end + end + return deepcopy(env_vjp(x)[2]) + end + return ∂f∂x, ∂f∂A +end + # compute the CTMRG gradient through fixed-point differentiation function _rrule( gradmode::FixedPointGradient, @@ -231,25 +290,13 @@ function _rrule( env, = leading_boundary(envinit, state, alg) - # prepare iterating function corresponding to a single gauge-fixed CTMRG iteration - alg_fixed = _set_fixed_truncation(alg) # fix spaces during differentiation - alg_gauge = _scrambling_env_gauge(alg) # select appropriate gauge-fixing algorithm - env_conv, info = ctmrg_iteration(InfiniteSquareNetwork(state), env, alg_fixed) - signs, corner_phases, edge_phases = compute_gauge_fix_gauge(env_conv, env, alg_gauge) - function gauge_fixed_iteration(A, x) - return fix_phases( - ctmrg_iteration(InfiniteSquareNetwork(A), x, alg_fixed)[1], - signs, corner_phases, edge_phases, - ) - end - # prepare its pullback - _, env_vjp = rrule_via_ad(config, gauge_fixed_iteration, state, env) - # split off state and environment parts - ∂f∂A(x)::typeof(state) = env_vjp(x)[2] - ∂f∂x(x)::typeof(env) = env_vjp(x)[3] + info, env_vjp = _fixedpoint_vjp(config, state, env, alg) function leading_boundary_fixed_pullback((Δenv′, Δinfo)) Δenv = unthunk(Δenv′) + Δenv isa AbstractZero && return NoTangent(), ZeroTangent(), ZeroTangent(), NoTangent() + Δenv = env isa SplitCTMRGEnv ? _project_split_tangent(Δenv, env) : Δenv + ∂f∂x, ∂f∂A = _fixedpoint_vjps(env_vjp, state, env) # evaluate the geometric sum ∂F∂env = fixedpoint_gradient(Δenv, ∂f∂x, ∂f∂A, Δenv, gradmode.solver_alg) diff --git a/src/environments/split_ctmrg_environments.jl b/src/environments/split_ctmrg_environments.jl index 574098d0b..72521169b 100644 --- a/src/environments/split_ctmrg_environments.jl +++ b/src/environments/split_ctmrg_environments.jl @@ -9,6 +9,7 @@ By default, initialization uses the outer environment space for this bond; conve Use `edge(env, :ket, dir, r, c)` and `edge(env, :bra, dir, r, c)` to access the factors, or `CTMRGEnv(env)` to contract them for existing observable routines. This environment supports [`SplitCTMRG`](@ref) and its fixed-point gradient. +Raw reverse-mode pullbacks may return structural cotangents; `ChainRulesCore.ProjectTo(env)(Δenv)` materializes a concrete environment for vector operations. """ struct SplitCTMRGEnv{C, E} corners::Array{C, 3} @@ -35,7 +36,8 @@ function ChainRulesCore.rrule( function split_env_pullback(Δenv_) Δenv = unthunk(Δenv_) Δenv isa AbstractZero && return NoTangent(), ZeroTangent(), ZeroTangent(), ZeroTangent(), NoTangent() - return NoTangent(), Δenv.corners, Δenv.edges_ket, Δenv.edges_bra, NoTangent() + return NoTangent(), _split_tangent_field(Δenv, :corners), + _split_tangent_field(Δenv, :edges_ket), _split_tangent_field(Δenv, :edges_bra), NoTangent() end return env, split_env_pullback end @@ -51,23 +53,52 @@ function ChainRulesCore.rrule(::typeof(getproperty), env::SplitCTMRGEnv, name::S function split_property_pullback(Δ_) Δ = unthunk(Δ_) Δ isa AbstractZero && return NoTangent(), ZeroTangent(), NoTangent() - Δ = _project_split_array(Δ, result) - corners = name === :corners ? Δ : zerovector.(env.corners) - edges_ket = name === :edges_ket ? Δ : zerovector.(env.edges_ket) - edges_bra = name === :edges_bra ? Δ : zerovector.(env.edges_bra) - return NoTangent(), SplitCTMRGEnv(corners, edges_ket, edges_bra, env.ket_first), NoTangent() + field = NamedTuple{(name,)}((Δ,)) + return NoTangent(), Tangent{typeof(env)}(; field..., ket_first = NoTangent()), NoTangent() end return result, split_property_pullback end -"""Project array cotangents and materialize zero entries in the primal tensor storage.""" -function _project_split_array(Δ, primal) - return map(primal, Δ) do tensor, delta - delta = unthunk(delta) - return delta isa AbstractZero ? zerovector(tensor) : ProjectTo(tensor)(delta) +"""Read a cotangent field, treating omitted fields and Zygote's nothing markers as zero.""" +function _split_tangent_field(Δ, name::Symbol) + Δ = unthunk(Δ) + (Δ === nothing || Δ isa AbstractZero) && return ZeroTangent() + fields = Δ isa Tangent ? ChainRulesCore.backing(Δ) : Δ + return fields isa NamedTuple ? get(fields, name, ZeroTangent()) : getfield(fields, name) +end + +"""Project accumulated array cotangents into fresh tensor storage, materializing zero entries once.""" +function _project_split_array(Δ, primal::Array{<:AbstractTensorMap, 3}) + Δ = unthunk(Δ) + output = similar(primal) + if Δ === nothing || Δ isa AbstractZero + return map!(zerovector, output, primal) end + axes(Δ) == axes(primal) || throw(DimensionMismatch("Split cotangent array axes must match the primal.")) + for I in eachindex(primal) + delta = unthunk(Δ[I]) + output[I] = delta === nothing || delta isa AbstractZero ? + zerovector(primal[I]) : copy(ProjectTo(primal[I])(delta)) + end + return output +end + +"""Materialize an owned split cotangent where differentiation or the fixed-point solver requires a concrete environment.""" +function _project_split_tangent(Δ, primal::SplitCTMRGEnv)::typeof(primal) + return SplitCTMRGEnv( + _project_split_array(_split_tangent_field(Δ, :corners), getfield(primal, :corners)), + _project_split_array(_split_tangent_field(Δ, :edges_ket), getfield(primal, :edges_ket)), + _project_split_array(_split_tangent_field(Δ, :edges_bra), getfield(primal, :edges_bra)), + getfield(primal, :ket_first), + ) end +ChainRulesCore.ProjectTo(env::SplitCTMRGEnv) = ProjectTo{SplitCTMRGEnv}(; primal = env) +(project::ProjectTo{SplitCTMRGEnv})(Δ::Union{NamedTuple, SplitCTMRGEnv}) = _project_split_tangent(Δ, project.primal) +(project::ProjectTo{SplitCTMRGEnv})(Δ::Tangent) = _project_split_tangent(Δ, project.primal) +# Resolve ChainRulesCore's identity projection for a tangent of the same primal type. +(project::ProjectTo{SplitCTMRGEnv})(Δ::Tangent{<:SplitCTMRGEnv}) = _project_split_tangent(Δ, project.primal) + """Factor an edge across its interlayer bond, distributing the singular values symmetrically.""" function _split_edge(E::CTMRGEdgeTensor{T, S, 3}, ket_first::Bool, trunc::TruncationStrategy) where {T, S} p = ket_first ? ((1, 2), (4, 3)) : ((1, 3), (4, 2)) diff --git a/test/testsuite/gradients/split_ctmrg_gradients.jl b/test/testsuite/gradients/split_ctmrg_gradients.jl index 63d82e81d..caadfcea2 100644 --- a/test/testsuite/gradients/split_ctmrg_gradients.jl +++ b/test/testsuite/gradients/split_ctmrg_gradients.jl @@ -1,7 +1,21 @@ using Test, Random, LinearAlgebra using TensorKit, PEPSKit, Zygote, Adapt +using ChainRulesCore: ProjectTo using PEPSKit: ctmrg_iteration, hook_pullback +"""Check public environment gradients against an analytic corner derivative.""" +function _split_environment_gradient(env::SplitCTMRGEnv) + """Read one corner so all other environment derivatives vanish.""" + loss(x) = real(dot(env.corners[1], x.corners[1])) + gradient = only(Zygote.withgradient(loss, env).grad) + @test gradient isa typeof(env) + @test only(Zygote.gradient(loss, env)) isa typeof(env) + @test gradient.corners[1] ≈ env.corners[1] + @test all(C -> iszero(norm(C)), gradient.corners[2:end]) + @test iszero(norm(gradient.edges_ket)) && iszero(norm(gradient.edges_bra)) + return nothing +end + function gradients_split_gauge(AT) return @testset "Split outer and interlayer gauge ($AT)" begin Random.seed!(32591) @@ -18,6 +32,10 @@ function gradients_split_gauge(AT) transformed = PEPSKit.fix_relative_phases(env, signs) restored = gauge_fix(transformed, env, PEPSKit.ScramblingEnvGauge()) @test norm(restored - env) / norm(env) < 1.0e-10 + _split_environment_gradient(env) + S = Vect[FermionParity](0 => 1, 1 => 1) + fermionic = adapt(AT, InfinitePEPS(S, S; unitcell = (2, 3))) + _split_environment_gradient(SplitCTMRGEnv(fermionic, S)) end end @@ -42,6 +60,7 @@ function gradients_split_sweep(AT) loss(A, x) = real(expectation_value(A, H, ctmrg_iteration(InfiniteSquareNetwork(A), x, alg)[1])) _, back = Zygote.pullback(loss, psi, env) dpsi, denv = back(1.0) + denv = ProjectTo(env)(denv) direction = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) direction_env = SplitCTMRGEnv(psi, ComplexSpace(2)) epsilon = 1.0e-6 @@ -61,10 +80,13 @@ function gradients_split_fixedpoint(AT) aligned = gauge_fix(next, env, PEPSKit.ScramblingEnvGauge()) @test norm(aligned - env) / norm(env) < 1.0e-7 gradalg = FixedPointGradient(; solver_alg = (; alg = :GMRES, tol = 1.0e-8, maxiter = 100)) - E, gs = Zygote.withgradient(psi) do A + E, back = Zygote.pullback(psi) do A boundary, = hook_pullback(leading_boundary, env, A, alg; alg_rrule = gradalg) return cost_function(A, boundary, H) end + gs = back(1.0) + @test only(back(2.0)) ≈ 2 * only(gs) atol = 1.0e-6 rtol = 1.0e-5 + @test only(back(1.0)) ≈ only(gs) atol = 1.0e-6 rtol = 1.0e-5 direction = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) epsilon = 1.0e-5 energy(A) = cost_function(A, leading_boundary(env, A, alg)[1], H) From b98c9deec759f2a618e33b63fc50d4029fd4cf83 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 13 Oct 2026 17:16:25 +0800 Subject: [PATCH 06/14] Reduce time to first run --- src/algorithms/ctmrg/projectors/split.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/ctmrg/projectors/split.jl b/src/algorithms/ctmrg/projectors/split.jl index 664ecb573..31279fcc9 100644 --- a/src/algorithms/ctmrg/projectors/split.jl +++ b/src/algorithms/ctmrg/projectors/split.jl @@ -183,11 +183,11 @@ function sequential_projectors( yellow_east, eye = _split_projectors_from_halves(col + 1, yellow_east_halves, env, alg, interlayer_trunc, Val(true), Val(true)) # Teal refines the outer cut using environments with green projectors inserted. - teal_halves = dtmap(r -> _split_teal_environment((r, col), network, env, green, green_east), axes(env, 2)) - teals = dtmap(r -> _split_teal_projector((r, col), teal_halves, env, alg), axes(env, 2)) + teal_halves = Base.inferencebarrier(dtmap)(r -> _split_teal_environment((r, col), network, env, green, green_east), axes(env, 2)) + teals = Base.inferencebarrier(dtmap)(r -> _split_teal_projector((r, col), teal_halves, env, alg), axes(env, 2)) # Red refines the interlayer cut using environments with yellow projectors inserted. - reds = dtmap(axes(env, 2)) do r + reds = Base.inferencebarrier(dtmap)(axes(env, 2)) do r quarters = _split_red_quarters((r, col), network, env, yellow, yellow_east) return _split_red_projector((r, col), quarters, env, alg; trunc = interlayer_trunc) end From 29181c97eec5f1f8b6f60b54b54d96c078a7f46c Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 09:10:02 +0800 Subject: [PATCH 07/14] Use Fishman projectors for split CTMRG --- docs/src/assets/pepskit.bib | 12 +++++++++++ docs/src/man/environments.md | 5 +++-- src/algorithms/ctmrg/projectors/split.jl | 21 +++++++++++++++--- src/algorithms/ctmrg/split.jl | 1 + test/testsuite/ctmrg/split.jl | 27 ++++++++++++++++++++++++ 5 files changed, 61 insertions(+), 5 deletions(-) diff --git a/docs/src/assets/pepskit.bib b/docs/src/assets/pepskit.bib index 0be60cf92..47e1c192e 100644 --- a/docs/src/assets/pepskit.bib +++ b/docs/src/assets/pepskit.bib @@ -24,6 +24,18 @@ @article{naumann_variationally_2025 url = {https://doi.org/10.1103/PhysRevB.111.235116} } +@article{fishman_faster_2018, + title = {Faster methods for contracting infinite two-dimensional tensor networks}, + author = {Fishman, M. T. and Vanderstraeten, L. and Zauner-Stauber, V. and Haegeman, J. and Verstraete, F.}, + journal = {Phys. Rev. B}, + volume = {98}, + pages = {235148}, + year = {2018}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.98.235148}, + url = {https://doi.org/10.1103/PhysRevB.98.235148} +} + @article{evenbly_gauge_2018, title = {Gauge fixing, canonical forms, and optimal truncations in tensor networks with closed loops}, author = {Evenbly, Glen}, diff --git a/docs/src/man/environments.md b/docs/src/man/environments.md index 498d3631f..b447395b1 100644 --- a/docs/src/man/environments.md +++ b/docs/src/man/environments.md @@ -7,7 +7,8 @@ [`SplitCTMRGEnv`](@ref) stores separate ket and bra factors for each CTMRG edge, connected by an interlayer bond. [`SplitCTMRG`](@ref) absorbs and projects the two layers sequentially using the full projectors of [Naumann et al., Phys. Rev. B 111, 235116 (2025)](https://doi.org/10.1103/PhysRevB.111.235116). -Its default [`SplitProjector`](@ref) constructs the four projector types directly from their split environment halves. +[`SplitProjector`](@ref) uses the Fishman construction for all four projector types, factorizing each environment half by SVD and keeping its square-root-weighted factor facing the cut before computing the projectors. +This improves the conditioning of the projector calculation at the cost of two additional SVDs per projector pair, following [Fishman et al.](@cite fishman_faster_2018), Appendix A. The outer and interlayer spaces can be chosen independently: ```julia @@ -23,7 +24,7 @@ When omitted, `interlayer_trunc` inherits `trunc`; the default `FixedSpaceTrunca `SplitCTMRGEnv(env::CTMRGEnv)` factorizes an existing environment exactly, while its `trunc` keyword can reduce the interlayer spaces. `CTMRGEnv(env::SplitCTMRGEnv)` contracts the interlayer bonds so the existing expectation-value routines can be reused. This version supports PEPS overlap networks, full projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref). -For comparisons with variPEPS, select its `Projector_Method.FULL` option; its default `FISHMAN` construction reweights the environment halves and uses different projectors. +The half-environment factors use direct square roots of their singular values; the projector pseudoinverse masks singular values below the relative cutoff `eps(real(scalartype(ψ)))^(3 / 4)`. Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref PEPSKit.ScramblingEnvGauge) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. The ordinary and split environments use the same boundary-MPS gauge method, with each split edge contributing two MPS sites instead of one. diff --git a/src/algorithms/ctmrg/projectors/split.jl b/src/algorithms/ctmrg/projectors/split.jl index 31279fcc9..01f2307b1 100644 --- a/src/algorithms/ctmrg/projectors/split.jl +++ b/src/algorithms/ctmrg/projectors/split.jl @@ -2,8 +2,9 @@ $(TYPEDEF) Projector algorithm for [`SplitCTMRG`](@ref), constructing the green (Fig. 3), teal (Fig. 7), yellow (Fig. 8), and red (Fig. 9) projector pairs of [Naumann et al.](@cite naumann_variationally_2025). -Each projector pair is constructed from a truncated SVD of the product of its environment halves. -QR compression is applied to the environment halves for the yellow projectors and, when the truncation permits, the quarter-environments for the red projectors before multiplication. +Each projector pair follows [Fishman et al.](@cite fishman_faster_2018), Appendix A: SVD-factorize each environment half, keep its square-root-weighted factor facing the cut, and truncate the product of these factors. +The half-environment singular values are weighted by direct square roots; the final projector pseudoinverse uses the relative cutoff of [`sdiag_pow`](@ref). +QR compression is applied to the environment halves for the yellow projectors and, when the truncation permits, the quarter-environments for the red projectors before factorization. ## Fields @@ -15,6 +16,7 @@ $(TYPEDFIELDS) Construct a split projector algorithm with `decomposition_alg`, `trunc`, and `verbosity`, following [`FullInfiniteProjector`](@ref). `trunc` controls the green and teal projectors; [`SplitCTMRG`](@ref) supplies `interlayer_trunc` for the yellow and red projectors. +The half-environment SVDs are untruncated; when `decomposition_alg` uses [`IterSVD`](@ref), these preliminary SVDs use the default dense algorithm instead. """ struct SplitProjector{S <: SVDAdjoint, T} <: ProjectorAlgorithm decomposition_alg::S @@ -27,6 +29,18 @@ end PROJECTOR_SYMBOLS[:SplitProjector] = SplitProjector +"""Factor split environment halves into square-root-weighted Fishman factors facing their shared cut.""" +function _split_fishman_halves(Q::AbstractTensorMap, Q_next::AbstractTensorMap, alg::SplitProjector) + half_alg = decomposition_algorithm(alg) + half_alg.fwd_alg isa IterSVD && (half_alg = SVDAdjoint()) + half_alg = @set half_alg.fwd_alg = TruncatedAlgorithm(half_alg.fwd_alg, notrunc()) + _, S, V = svd_trunc_no_error(Q, half_alg) + U_next, S_next, _ = svd_trunc_no_error(Q_next, half_alg) + Q = sdiag_pow(S, 0.5) * V + Q_next = U_next * sdiag_pow(S_next, 0.5) + return Q / norm(Q), Q_next / norm(Q_next) +end + """Compute a split projector pair, optionally QR-compressing the environment halves for the yellow or red projectors.""" function compute_projector( halves, alg::SplitProjector; @@ -41,7 +55,8 @@ function compute_projector( Q_next, _ = right_orth(Q_next, LQAdjoint()) end - # All four colors use an SVD of the product of their prepared halves. + # All four colors use Fishman factors before the final truncated SVD. + Q, Q_next = _split_fishman_halves(Q, Q_next, alg) environment = Q * Q_next svd_alg = decomposition_algorithm(_set_decomposition_truncation(alg, trunc)) U, S, V, truncation_error = svd_trunc!(environment / norm(environment), svd_alg) diff --git a/src/algorithms/ctmrg/split.jl b/src/algorithms/ctmrg/split.jl index 582237243..e8165a048 100644 --- a/src/algorithms/ctmrg/split.jl +++ b/src/algorithms/ctmrg/split.jl @@ -6,6 +6,7 @@ Uses [`SplitCTMRGEnv`](@ref), with independently truncated outer and interlayer The common iteration, `trunc`, and `decomposition_alg` keywords follow [`SequentialCTMRG`](@ref). `trunc` controls the green and teal projectors; `interlayer_trunc` controls the yellow and red projectors and defaults to `trunc`. Uses [`SplitProjector`](@ref) to construct the split projector pairs. +All four projector types use the square-root factorization of [Fishman et al.](@cite fishman_faster_2018), Appendix A, for numerical stability. The fixed-point gradient aligns the outer and interlayer gauges and holds both environment spaces fixed. """ struct SplitCTMRG{P <: SplitProjector, I <: TruncationStrategy} <: CTMRGAlgorithm diff --git a/test/testsuite/ctmrg/split.jl b/test/testsuite/ctmrg/split.jl index 4f9948448..01f2c07df 100644 --- a/test/testsuite/ctmrg/split.jl +++ b/test/testsuite/ctmrg/split.jl @@ -2,16 +2,29 @@ using Test, Random, LinearAlgebra using TensorKit, MatrixAlgebraKit, PEPSKit, Adapt using PEPSKit: ctmrg_iteration, WEST +"""Construct a dense Fishman cut projector independently from fourth roots of the Gram matrices.""" +function _split_fishman_reference(Q::AbstractMatrix, Q_next::AbstractMatrix, rank::Int) + A = normalize(Matrix(Hermitian(Q' * Q)^0.25)) + B = normalize(Matrix(Hermitian(Q_next * Q_next')^0.25)) + U, S, V = svd(normalize(A * B)) + return B * V[:, 1:rank] * Diagonal(1 ./ S[1:rank]) * U[:, 1:rank]' * A +end + function ctmrg_split_projectors(AT) Random.seed!(37913) return @testset "Split projector stability ($AT)" begin alg = SplitCTMRG(; verbosity = 0).projector_alg + iterative = SplitProjector(; decomposition_alg = (; fwd_alg = (; alg = :GKL, krylovdim = 10)), verbosity = 0) trunc = truncspace(ComplexSpace(3)) # Cover the two-index yellow and three-index red cuts. for cut in (ComplexSpace(3) ⊗ ComplexSpace(2), ComplexSpace(3) ⊗ ComplexSpace(2) ⊗ ComplexSpace(2)) Q = adapt(AT, randn(ComplexF64, ComplexSpace(12) ← cut)) Q_next = adapt(AT, randn(ComplexF64, cut ← ComplexSpace(12))) reference, refinfo = PEPSKit.compute_projector((Q, Q_next), alg; trunc) + dense = _split_fishman_reference( + reshape(convert(Array, Q), 12, :), reshape(convert(Array, Q_next), :, 12), 3, + ) + @test reshape(convert(Array, reference[1] * reference[2]), size(dense)) ≈ dense compressed, info = PEPSKit.compute_projector((Q, Q_next), alg; trunc, compress = true) # QR compression preserves the cut projector and the nonzero singular spectrum. @test compressed[1] * compressed[2] ≈ reference[1] * reference[2] @@ -21,6 +34,20 @@ function ctmrg_split_projectors(AT) @test scaled[1] * scaled[2] ≈ reference[1] * reference[2] end end + + # Check conditioning, projector pseudoinverse masking, and fixed output spaces with both final solvers. + V = ComplexSpace(4) + spectra = (([1.0, 1.0e-3, 1.0e-6, 1.0e-9], 4), ([1.0, 1.0e-6, 1.0e-15, 0.0], 2)) + for (values, retained) in spectra, (solver, rank) in ((alg, 4), (iterative, 1)) + Q = adapt(AT, TensorMap(Matrix(Diagonal(values)), V ← V)) + projectors, info = PEPSKit.compute_projector((Q, Q), solver; trunc = truncspace(ComplexSpace(rank))) + spectrum = values + support = [ones(min(retained, rank)); zeros(4 - min(retained, rank))] + @test Array(info.S.data) ≈ normalize(spectrum)[1:rank] + @test count(!iszero, info.S.data) == min(count(!iszero, values), rank) + @test convert(Array, projectors[1] * projectors[2]) ≈ (norm(spectrum) / sum(spectrum)) * Diagonal(support) + @test all(P -> isfinite(norm(P)), projectors) + end end end From f3dfa0205a102a93a340b65390879db383cbe3fe Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 13:56:29 +0800 Subject: [PATCH 08/14] Norm-preserving preconditioning --- .../optimization/peps_optimization.jl | 36 +++++++---- src/utility/retractions.jl | 32 +++++----- test/optimization/norm_preservation.jl | 12 ++++ test/testsuite/TestSuite.jl | 8 +++ .../optimization/norm_preservation.jl | 62 +++++++++++++++++++ test/testsuite/utility/retractions.jl | 20 +++++- 6 files changed, 140 insertions(+), 30 deletions(-) create mode 100644 test/optimization/norm_preservation.jl create mode 100644 test/testsuite/optimization/norm_preservation.jl diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index ac950fe56..645193581 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -115,6 +115,7 @@ end Find the fixed point of `operator` (i.e. the ground state) starting from `peps₀` according to the supplied optimization parameters. The initial environment `env₀` serves as an initial guess for the first boundary contraction run. By default, a random initial environment is used. +Each local tensor is normalized initially, and optimization preserves its Euclidean norm by projecting gradients and preconditioned directions onto its real tangent space. The optimization parameters can be supplied via the keyword arguments or directly as a `PEPSOptimize` struct. The following keyword arguments are supported: @@ -271,7 +272,8 @@ function fixedpoint( alg.precondition_alg; iter = tracked_finalizer.tol_state.iter, g_global = tracked_finalizer.tol_state.gradnorm ) - return peps_precondition(x, g, tracked_finalizer.tol_state, precondition_alg) + g_prec = peps_precondition(x, g, tracked_finalizer.tol_state, precondition_alg) + return peps_project_tangent(g_prec, x[1]) end # optimize operator cost function @@ -302,7 +304,7 @@ function fixedpoint( end return cost_function(ψ, env′, operator) end - g = only(gs) # `withgradient` returns tuple of gradients `gs` + g = peps_project_tangent(only(gs), peps) tracked_finalizer.gradnorms_unitcell[end] = norm.(unitcell(g)) tracked_finalizer.times[end] = time() - start_time return E, g @@ -352,11 +354,20 @@ function peps_normalize(A::InfinitePEPS) return InfinitePEPS(normalized_tensors) end +""" +Project each direction onto the real tangent space of its fixed-norm PEPS tensor without mutating the inputs. +""" +function peps_project_tangent(η::InfinitePEPS, peps::InfinitePEPS) + tangents = map(unitcell(η), unitcell(peps)) do direction, A + return add(direction, A, -real(inner(A, direction)) / norm(A)^2) + end + return InfinitePEPS(tangents) +end + """ $(SIGNATURES) -Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step -size `α`, giving a new PEPS `A´`, +Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step size `α`, giving a new PEPS `A´`, ```math A' ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, ``` @@ -364,7 +375,8 @@ and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +Each tensor and its direction must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = 0``. """ function peps_retract(x, η, α) peps = x[1] @@ -380,17 +392,17 @@ end """ $(SIGNATURES) -Transports a direction at `A = x[1]` to a valid direction at `A´ = x´[1]` corresponding to -the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting -from a direction `η` of the form +Transports a direction at `A = x[1]` to a valid direction at `A´ = x´[1]` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. +In particular, starting from a direction `ξ` of the form ```math -ξ = ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ +ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ ``` -where ``⟨ Δξ, A ⟩ = ⟨ Δξ, η ⟩ = 0``, it returns +where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ = 0``, it returns ```math -ξ(α) = ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ +ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. +Each tensor and its directions `η` and `ξ` must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = \\operatorname{Re} ⟨ A, ξ ⟩ = 0``. """ function peps_transport!(ξ, x, η, α, x´) peps = x[1] diff --git a/src/utility/retractions.jl b/src/utility/retractions.jl index fe4aee991..75a6e532b 100644 --- a/src/utility/retractions.jl +++ b/src/utility/retractions.jl @@ -5,8 +5,7 @@ Utilities for preserving the norm of (VectorInterface-compliant) vectors during """ $(SIGNATURES) -Performs a norm-preserving retraction of vector `A` along the direction `η` with step size -`α`, giving a new vector `A´`, +Performs a norm-preserving retraction of vector `A` along the direction `η` with step size `α`, giving a new vector `A´`, ```math A' ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, ``` @@ -14,16 +13,18 @@ and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +The inputs must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = 0``. +For `η = 0`, return copies of `A` and `η`. !!! note - The vectors `A` and `η` should satisfy the interface specified by - [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl) + The vectors `A` and `η` should satisfy the interface specified by [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl). """ function norm_preserving_retract(A, η, α) n_A = norm(A) n_η = norm(η) + iszero(n_η) && return scale(A, true), zerovector(η) sn, cs = sincos(α * n_η / n_A) A´ = add(A, η, sn * n_A / n_η, cs) @@ -35,29 +36,30 @@ end """ $(SIGNATURES) -Transports a direction `ξ` at `A` to a valid direction at `A´` corresponding to -the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting -from a direction `η` of the form +Transports a direction `ξ` at `A` to a valid direction at `A´` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. +In particular, starting from a direction `ξ` of the form ```math -ξ = ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ +ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ ``` -where ``⟨ Δξ, A ⟩ = ⟨ Δξ, η ⟩ = 0``, it returns +where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ = 0``, it returns ```math -ξ(α) = ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin ( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ +ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin ( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. +Both `η` and `ξ` must be real-orthogonal to `A`. +For `η = 0`, return `ξ` unchanged. !!! note - The vectors `A` and `η` should satisfy the interface specified by - [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl) + The vectors `A` and `η` should satisfy the interface specified by [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl). """ function norm_preserving_transport!(ξ, A, η, α, A´) n_A = norm(A) n_η = norm(η) + iszero(n_η) && return ξ sn, cs = sincos(α * n_η / n_A) - overlaps_η_ξ = inner(η, ξ) / n_η + overlaps_η_ξ = real(inner(η, ξ)) / n_η add!(ξ, η, (cs - 1) * overlaps_η_ξ / n_η) add!(ξ, A, -sn * overlaps_η_ξ / n_A) diff --git a/test/optimization/norm_preservation.jl b/test/optimization/norm_preservation.jl new file mode 100644 index 000000000..a145faf9c --- /dev/null +++ b/test/optimization/norm_preservation.jl @@ -0,0 +1,12 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.optimization_tangent_projection(Vector) + TestSuite.optimization_norm_preservation(Vector) +end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..44eeebf4a 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -197,6 +197,14 @@ module GradientsCTMRG end using .GradientsCTMRG +# Optimization +# ------------ +module OptimizationNormPreservation + include("optimization/norm_preservation.jl") + export optimization_tangent_projection, optimization_norm_preservation +end +using .OptimizationNormPreservation + # Time evolution # -------------- module TimeEvolClusterProjectors diff --git a/test/testsuite/optimization/norm_preservation.jl b/test/testsuite/optimization/norm_preservation.jl new file mode 100644 index 000000000..41ec3d4f0 --- /dev/null +++ b/test/testsuite/optimization/norm_preservation.jl @@ -0,0 +1,62 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using PEPSKit +using Adapt + +"""Check local tangent projection and norm-preserving retraction of PEPS tensors.""" +function optimization_tangent_projection(AT) + return @testset "PEPS tangent projection ($T, $AT)" for T in (Float64, ComplexF64) + Random.seed!(1234) + peps = adapt(AT, InfinitePEPS(randn, T, ℂ^2, ℂ^2; unitcell = (1, 2))) + rmul!(peps[1, 2], 3) + direction = adapt(AT, InfinitePEPS(randn, T, ℂ^2, ℂ^2; unitcell = (1, 2))) + original_peps = deepcopy(peps) + original_direction = deepcopy(direction) + + tangent = PEPSKit.peps_project_tangent(direction, peps) + (retracted, _), _ = PEPSKit.peps_retract((peps, nothing), tangent, 0.3) + for i in eachindex(peps) + @test real(dot(peps[i], tangent[i])) ≈ 0 atol = 1.0e-12 + @test norm(retracted[i]) ≈ norm(peps[i]) rtol = 1.0e-12 + @test peps[i] == original_peps[i] + @test direction[i] == original_direction[i] + end + if T <: Complex + phase_direction = im * peps + @test norm(PEPSKit.peps_project_tangent(phase_direction, peps) - phase_direction) < 1.0e-12 + end + end +end + +"""Check local tensor norms during CTMRG optimization with and without preconditioning.""" +function optimization_norm_preservation(AT) + return @testset "Traditional CTMRG optimization preserves local tensor norms ($AT)" begin + Random.seed!(123) + H = adapt(AT, heisenberg_XYZ(InfiniteSquare(); Jx = -1, Jy = 1, Jz = -1)) + peps₀ = adapt(AT, InfinitePEPS(randn, ComplexF64, ℂ^2, ℂ^2)) + env₀ = adapt(AT, CTMRGEnv(peps₀, ℂ^4)) + + @testset "preconditioner = $precondition_alg" for precondition_alg in (nothing, (;)) + norms = Float64[] + radial_gradients = Float64[] + """Record local tensor norms and radial gradient components after each accepted step.""" + function record_norms((peps, env), E, g, iter) + append!(norms, norm.(PEPSKit.unitcell(peps))) + append!(radial_gradients, abs.(real.(dot.(PEPSKit.unitcell(peps), PEPSKit.unitcell(g))))) + return (peps, env), E, g + end + fixedpoint( + H, peps₀, deepcopy(env₀); + boundary_alg = (; alg = :SequentialCTMRG, maxiter = 100), + gradient_alg = (; solver_alg = (; alg = :GMRES)), + optimizer_alg = (; maxiter = 4, tol = 1.0e-8), + precondition_alg, verbosity = 0, (finalize!) = record_norms, + ) + @test !isempty(norms) + @test all(n -> isapprox(n, 1; atol = 1.0e-12), norms) + @test all(r -> r < 1.0e-12, radial_gradients) + end + end +end diff --git a/test/testsuite/utility/retractions.jl b/test/testsuite/utility/retractions.jl index c4e0311c0..400b73622 100644 --- a/test/testsuite/utility/retractions.jl +++ b/test/testsuite/utility/retractions.jl @@ -10,6 +10,7 @@ dtype = ComplexF64 Vphyss = [ℂ^2, U1Space(0 => 1, -1 => 1, 1 => 1)] Vpepss = [ℂ^4, U1Space(0 => 2, -1 => 1, 1 => 1)] +"""Check real-sphere retractions and isometric transport of tensor directions.""" function utility_retractions(AT) return @testset "Norm-preserving tensor retractions for sectortype $(sectortype(Vphyss[i])) ($AT)" for i in eachindex( @@ -24,13 +25,26 @@ function utility_retractions(AT) normalized_A = scale(A, inv(norm(A))) η = adapt(AT, randn(dtype, peps_space)) ζ = adapt(AT, randn(dtype, peps_space)) - add!(η, normalized_A, -inner(normalized_A, η)) - add!(ζ, normalized_A, -inner(normalized_A, ζ)) + add!(η, normalized_A, -real(inner(normalized_A, η))) + add!(ζ, normalized_A, -real(inner(normalized_A, ζ))) A´, ξ = PEPSKit.norm_preserving_retract(A, η, α) @test norm(A´) ≈ norm(A) rtol = 1.0e-12 + @test real(inner(A´, ξ)) ≈ 0 atol = 1.0e-12 + norm_ζ = norm(ζ) PEPSKit.norm_preserving_transport!(ζ, A, η, α, A´) - @test inner(ζ, A´) ≈ 0 atol = 1.0e-12 + @test real(inner(ζ, A´)) ≈ 0 atol = 1.0e-12 + @test norm(ζ) ≈ norm_ζ rtol = 1.0e-12 + + for point in (A, [A], VectorInterface.MinimalMVec([A])) + norm_point = norm(point) + unchanged_A, zero_derivative = PEPSKit.norm_preserving_retract(point, zerovector(point), α) + @test norm(unchanged_A) ≈ norm_point + @test norm(zero_derivative) == 0 + scale!(unchanged_A, 0) + @test norm(point) ≈ norm_point + end + @test PEPSKit.norm_preserving_transport!(ζ, A, zerovector(A), α, A) === ζ end end From 813500dc56ee3d330aa73a33b2173916ca10d486 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 15:33:07 +0800 Subject: [PATCH 09/14] Use `\prime` instead of backticks or single-quote --- .../optimization/peps_optimization.jl | 26 +++++++++---------- src/utility/retractions.jl | 16 ++++++------ 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 645193581..6649b15ff 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -367,15 +367,15 @@ end """ $(SIGNATURES) -Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step size `α`, giving a new PEPS `A´`, +Performs a norm-preserving retraction of an infinite PEPS `A = x[1]` along `η` with step size `α`, giving a new PEPS ``A′``, ```math -A' ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, +A′ ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, ``` and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A′, ξ ⟩ = 0`` and ``‖A′‖ = ‖A‖``. Each tensor and its direction must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = 0``. """ function peps_retract(x, η, α) @@ -383,16 +383,16 @@ function peps_retract(x, η, α) env = deepcopy(x[2]) retractions = norm_preserving_retract.(unitcell(peps), unitcell(η), α) - peps´ = InfinitePEPS(map(first, retractions)) + peps′ = InfinitePEPS(map(first, retractions)) ξ = InfinitePEPS(map(last, retractions)) - return (peps´, env), ξ + return (peps′, env), ξ end """ $(SIGNATURES) -Transports a direction at `A = x[1]` to a valid direction at `A´ = x´[1]` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. +Transports a direction at `A = x[1]` to a valid direction at ``A′ = x′[1]`` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting from a direction `ξ` of the form ```math ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ @@ -401,15 +401,15 @@ where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ ```math ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A′, ξ(α) ⟩ = 0``. Each tensor and its directions `η` and `ξ` must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = \\operatorname{Re} ⟨ A, ξ ⟩ = 0``. """ -function peps_transport!(ξ, x, η, α, x´) +function peps_transport!(ξ, x, η, α, x′) peps = x[1] - peps´ = x´[1] + peps′ = x′[1] norm_preserving_transport!.( - unitcell(ξ), unitcell(peps), unitcell(η), α, unitcell(peps´) + unitcell(ξ), unitcell(peps), unitcell(η), α, unitcell(peps′) ) return ξ @@ -434,11 +434,11 @@ function symmetrize_retract_and_finalize!( end function retract_then_symmetrize((peps, env), η, α) # retract - (peps´, env´), ξ = retract((peps, env), η, α) + (peps′, env′), ξ = retract((peps, env), η, α) # symmetrize retracted point and directional derivative - peps´_symm = symmetrize!(peps´, symm) + peps′_symm = symmetrize!(peps′, symm) ξ_symm = symmetrize!(ξ, symm) - return (peps´_symm, env´), ξ_symm + return (peps′_symm, env′), ξ_symm end return retract_then_symmetrize, symmetrize_then_finalize! end diff --git a/src/utility/retractions.jl b/src/utility/retractions.jl index 75a6e532b..9de4d5f6d 100644 --- a/src/utility/retractions.jl +++ b/src/utility/retractions.jl @@ -5,15 +5,15 @@ Utilities for preserving the norm of (VectorInterface-compliant) vectors during """ $(SIGNATURES) -Performs a norm-preserving retraction of vector `A` along the direction `η` with step size `α`, giving a new vector `A´`, +Performs a norm-preserving retraction of vector `A` along the direction `η` with step size `α`, giving a new vector ``A′``, ```math -A' ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, +A′ ← \\cos ( α ‖η‖ / ‖A‖ ) A + \\sin ( α ‖η‖ / ‖A‖ ) ‖A‖ η / ‖η‖, ``` and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A′, ξ ⟩ = 0`` and ``‖A′‖ = ‖A‖``. The inputs must satisfy ``\\operatorname{Re} ⟨ A, η ⟩ = 0``. For `η = 0`, return copies of `A` and `η`. @@ -27,16 +27,16 @@ function norm_preserving_retract(A, η, α) iszero(n_η) && return scale(A, true), zerovector(η) sn, cs = sincos(α * n_η / n_A) - A´ = add(A, η, sn * n_A / n_η, cs) + A′ = add(A, η, sn * n_A / n_η, cs) ξ = add(A, η, cs, -sn * n_η / n_A) - return A´, ξ + return A′, ξ end """ $(SIGNATURES) -Transports a direction `ξ` at `A` to a valid direction at `A´` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. +Transports a direction `ξ` at `A` to a valid direction at ``A′`` corresponding to the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting from a direction `ξ` of the form ```math ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ @@ -45,7 +45,7 @@ where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ ```math ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin ( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A′, ξ(α) ⟩ = 0``. Both `η` and `ξ` must be real-orthogonal to `A`. For `η = 0`, return `ξ` unchanged. @@ -53,7 +53,7 @@ For `η = 0`, return `ξ` unchanged. The vectors `A` and `η` should satisfy the interface specified by [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl). """ -function norm_preserving_transport!(ξ, A, η, α, A´) +function norm_preserving_transport!(ξ, A, η, α, A′) n_A = norm(A) n_η = norm(η) iszero(n_η) && return ξ From 8c60dc9b4493004eb0d23166672957a3b3109e82 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 15:34:43 +0800 Subject: [PATCH 10/14] Fix regularization formula in docstring --- src/algorithms/optimization/preconditioning.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/optimization/preconditioning.jl b/src/algorithms/optimization/preconditioning.jl index 31494d7d1..c284b6768 100644 --- a/src/algorithms/optimization/preconditioning.jl +++ b/src/algorithms/optimization/preconditioning.jl @@ -40,7 +40,7 @@ which is solved separately for every tensor in the unit cell using an iterative The regularization ``δ`` prevents the (generally singular) metric from being inverted exactly, and decays during the optimization according to ```math -δ = \\texttt{regularization} × ‖∇f‖^2 / \\texttt{iter}^2 +δ = \\texttt{regularization} × ‖∇f‖^2 / \\max(\\texttt{iter}, 1) ``` so that the preconditioner acts conservatively far from the minimum and becomes increasingly aggressive as the optimization converges. From 54cd1bb8751418d220e2ca355c032c4398ee3ac6 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 15:41:33 +0800 Subject: [PATCH 11/14] More backtick to prime renaming --- test/testsuite/utility/retractions.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/testsuite/utility/retractions.jl b/test/testsuite/utility/retractions.jl index 400b73622..1c44d5bd2 100644 --- a/test/testsuite/utility/retractions.jl +++ b/test/testsuite/utility/retractions.jl @@ -28,13 +28,13 @@ function utility_retractions(AT) add!(η, normalized_A, -real(inner(normalized_A, η))) add!(ζ, normalized_A, -real(inner(normalized_A, ζ))) - A´, ξ = PEPSKit.norm_preserving_retract(A, η, α) - @test norm(A´) ≈ norm(A) rtol = 1.0e-12 + A′, ξ = PEPSKit.norm_preserving_retract(A, η, α) + @test norm(A′) ≈ norm(A) rtol = 1.0e-12 - @test real(inner(A´, ξ)) ≈ 0 atol = 1.0e-12 + @test real(inner(A′, ξ)) ≈ 0 atol = 1.0e-12 norm_ζ = norm(ζ) - PEPSKit.norm_preserving_transport!(ζ, A, η, α, A´) - @test real(inner(ζ, A´)) ≈ 0 atol = 1.0e-12 + PEPSKit.norm_preserving_transport!(ζ, A, η, α, A′) + @test real(inner(ζ, A′)) ≈ 0 atol = 1.0e-12 @test norm(ζ) ≈ norm_ζ rtol = 1.0e-12 for point in (A, [A], VectorInterface.MinimalMVec([A])) From 5e206afbe70c3258a40617bca299379a9b88a722 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 17:31:33 +0800 Subject: [PATCH 12/14] Avoid redundant second projection when preconditioning is off --- src/algorithms/optimization/peps_optimization.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 6649b15ff..a72bb51e6 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -268,6 +268,7 @@ function fixedpoint( # initialize the preconditioner function precondition(x, g) + isnothing(alg.precondition_alg) && return g precondition_alg = adapt_solver( alg.precondition_alg; iter = tracked_finalizer.tol_state.iter, g_global = tracked_finalizer.tol_state.gradnorm @@ -359,7 +360,7 @@ Project each direction onto the real tangent space of its fixed-norm PEPS tensor """ function peps_project_tangent(η::InfinitePEPS, peps::InfinitePEPS) tangents = map(unitcell(η), unitcell(peps)) do direction, A - return add(direction, A, -real(inner(A, direction)) / norm(A)^2) + return add(direction, A, -real(inner(A, direction)) / real(inner(A, A))) end return InfinitePEPS(tangents) end From ec998c18a68e22110b338fb8c08b2cd0b25a15b7 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 15 Oct 2026 17:49:40 +0800 Subject: [PATCH 13/14] Preserve local tensor norms during PEPS optimization --- .../optimization/peps_optimization.jl | 26 ++++++++++++++----- src/utility/retractions.jl | 14 +++++----- test/optimization/norm_preservation.jl | 19 ++++++++++++++ test/testsuite/utility/retractions.jl | 26 ++++++++++++++----- 4 files changed, 66 insertions(+), 19 deletions(-) create mode 100644 test/optimization/norm_preservation.jl diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 1ad9bad81..84079c574 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -115,6 +115,7 @@ end Find the fixed point of `operator` (i.e. the ground state) starting from `peps₀` according to the supplied optimization parameters. The initial environment `env₀` serves as an initial guess for the first boundary contraction run. By default, a random initial environment is used. +Each local tensor is normalized initially, and optimization preserves its Euclidean norm by projecting gradients and preconditioned directions onto its real tangent space. The optimization parameters can be supplied via the keyword arguments or directly as a `PEPSOptimize` struct. The following keyword arguments are supported: @@ -271,7 +272,8 @@ function fixedpoint( alg.precondition_alg; iter = tracked_finalizer.tol_state.iter, g_global = tracked_finalizer.tol_state.gradnorm ) - return peps_precondition(x, g, tracked_finalizer.tol_state, precondition_alg) + g_prec = peps_precondition(x, g, tracked_finalizer.tol_state, precondition_alg) + return peps_project_tangent(g_prec, x[1]) end # optimize operator cost function @@ -302,7 +304,7 @@ function fixedpoint( end return cost_function(ψ, env′, operator) end - g = only(gs) # `withgradient` returns tuple of gradients `gs` + g = peps_project_tangent(only(gs), peps) tracked_finalizer.gradnorms_unitcell[end] = norm.(unitcell(g)) tracked_finalizer.times[end] = time() - start_time return E, g @@ -352,6 +354,16 @@ function peps_normalize(A::InfinitePEPS) return InfinitePEPS(normalized_tensors) end +""" +Project each direction onto the real tangent space of its fixed-norm PEPS tensor without mutating the inputs. +""" +function peps_project_tangent(η::InfinitePEPS, peps::InfinitePEPS) + tangents = map(unitcell(η), unitcell(peps)) do direction, A + return add(direction, A, -real(inner(A, direction)) / norm(A)^2) + end + return InfinitePEPS(tangents) +end + """ $(SIGNATURES) @@ -364,7 +376,7 @@ and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. """ function peps_retract(x, η, α) peps = x[1] @@ -384,13 +396,13 @@ Transports a direction at `A = x[1]` to a valid direction at `A´ = x´[1]` corr the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting from a direction `η` of the form ```math -ξ = ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ +ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ ``` -where ``⟨ Δξ, A ⟩ = ⟨ Δξ, η ⟩ = 0``, it returns +where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ = 0``, it returns ```math -ξ(α) = ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ +ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. """ function peps_transport!(ξ, x, η, α, x´) peps = x[1] diff --git a/src/utility/retractions.jl b/src/utility/retractions.jl index fe4aee991..86a7078c7 100644 --- a/src/utility/retractions.jl +++ b/src/utility/retractions.jl @@ -14,7 +14,7 @@ and corresponding directional derivative `ξ`, ```math ξ = \\cos ( α ‖η‖ / ‖A‖ ) η - \\sin ( α ‖η‖ / ‖A‖ ) ‖η‖ A / ‖A‖, ``` -such that ``⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. +such that ``\\operatorname{Re} ⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. !!! note The vectors `A` and `η` should satisfy the interface specified by @@ -24,6 +24,7 @@ such that ``⟨ A', ξ ⟩ = 0`` and ``‖A'‖ = ‖A‖``. function norm_preserving_retract(A, η, α) n_A = norm(A) n_η = norm(η) + iszero(n_η) && return scale(A, true), zerovector(η) sn, cs = sincos(α * n_η / n_A) A´ = add(A, η, sn * n_A / n_η, cs) @@ -39,13 +40,13 @@ Transports a direction `ξ` at `A` to a valid direction at `A´` corresponding t the norm-preserving retraction of `A` along `η` with step size `α`. In particular, starting from a direction `η` of the form ```math -ξ = ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ +ξ = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ η / ‖η‖ + Δξ ``` -where ``⟨ Δξ, A ⟩ = ⟨ Δξ, η ⟩ = 0``, it returns +where ``\\operatorname{Re} ⟨ Δξ, A ⟩ = \\operatorname{Re} ⟨ Δξ, η ⟩ = 0``, it returns ```math -ξ(α) = ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin ( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ +ξ(α) = \\operatorname{Re} ⟨ η / ‖η‖, ξ ⟩ ( \\cos ( α ‖η‖ / ‖A‖ ) η / ‖η‖ - \\sin ( α ‖η‖ / ‖A‖ ) A / ‖A‖ ) + Δξ ``` -such that ``‖ξ(α)‖ = ‖ξ‖, ⟨ A', ξ(α) ⟩ = 0``. +such that ``‖ξ(α)‖ = ‖ξ‖, \\operatorname{Re} ⟨ A', ξ(α) ⟩ = 0``. !!! note The vectors `A` and `η` should satisfy the interface specified by @@ -55,9 +56,10 @@ such that ``‖ξ(α)‖ = ‖ξ‖, ⟨ A', ξ(α) ⟩ = 0``. function norm_preserving_transport!(ξ, A, η, α, A´) n_A = norm(A) n_η = norm(η) + iszero(n_η) && return ξ sn, cs = sincos(α * n_η / n_A) - overlaps_η_ξ = inner(η, ξ) / n_η + overlaps_η_ξ = real(inner(η, ξ)) / n_η add!(ξ, η, (cs - 1) * overlaps_η_ξ / n_η) add!(ξ, A, -sn * overlaps_η_ξ / n_A) diff --git a/test/optimization/norm_preservation.jl b/test/optimization/norm_preservation.jl new file mode 100644 index 000000000..3b4c709b4 --- /dev/null +++ b/test/optimization/norm_preservation.jl @@ -0,0 +1,19 @@ +using Test, Random, LinearAlgebra, TensorKit, PEPSKit + +@testset "Norm-preserving optimization geometry" begin + Random.seed!(1234) + peps = PEPSKit.peps_normalize(InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell = (1, 2))) + direction = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell = (1, 2)) + tangent = PEPSKit.peps_project_tangent(direction, peps) + (retracted, _), velocity = PEPSKit.peps_retract((peps, nothing), tangent, 0.3) + transported = deepcopy(tangent) + PEPSKit.peps_transport!(transported, (peps, nothing), tangent, 0.3, (retracted, nothing)) + for i in eachindex(peps) + @test abs(real(dot(peps[i], tangent[i]))) < 1e-12 + @test norm(retracted[i]) ≈ 1 atol = 1e-12 + @test norm(transported[i]) ≈ norm(tangent[i]) rtol = 1e-12 + @test abs(real(dot(retracted[i], transported[i]))) < 1e-12 + @test transported[i] ≈ velocity[i] + end + @test PEPSKit.peps_project_tangent(im * peps, peps) ≈ im * peps +end diff --git a/test/testsuite/utility/retractions.jl b/test/testsuite/utility/retractions.jl index c4e0311c0..1c44d5bd2 100644 --- a/test/testsuite/utility/retractions.jl +++ b/test/testsuite/utility/retractions.jl @@ -10,6 +10,7 @@ dtype = ComplexF64 Vphyss = [ℂ^2, U1Space(0 => 1, -1 => 1, 1 => 1)] Vpepss = [ℂ^4, U1Space(0 => 2, -1 => 1, 1 => 1)] +"""Check real-sphere retractions and isometric transport of tensor directions.""" function utility_retractions(AT) return @testset "Norm-preserving tensor retractions for sectortype $(sectortype(Vphyss[i])) ($AT)" for i in eachindex( @@ -24,13 +25,26 @@ function utility_retractions(AT) normalized_A = scale(A, inv(norm(A))) η = adapt(AT, randn(dtype, peps_space)) ζ = adapt(AT, randn(dtype, peps_space)) - add!(η, normalized_A, -inner(normalized_A, η)) - add!(ζ, normalized_A, -inner(normalized_A, ζ)) + add!(η, normalized_A, -real(inner(normalized_A, η))) + add!(ζ, normalized_A, -real(inner(normalized_A, ζ))) - A´, ξ = PEPSKit.norm_preserving_retract(A, η, α) - @test norm(A´) ≈ norm(A) rtol = 1.0e-12 + A′, ξ = PEPSKit.norm_preserving_retract(A, η, α) + @test norm(A′) ≈ norm(A) rtol = 1.0e-12 - PEPSKit.norm_preserving_transport!(ζ, A, η, α, A´) - @test inner(ζ, A´) ≈ 0 atol = 1.0e-12 + @test real(inner(A′, ξ)) ≈ 0 atol = 1.0e-12 + norm_ζ = norm(ζ) + PEPSKit.norm_preserving_transport!(ζ, A, η, α, A′) + @test real(inner(ζ, A′)) ≈ 0 atol = 1.0e-12 + @test norm(ζ) ≈ norm_ζ rtol = 1.0e-12 + + for point in (A, [A], VectorInterface.MinimalMVec([A])) + norm_point = norm(point) + unchanged_A, zero_derivative = PEPSKit.norm_preserving_retract(point, zerovector(point), α) + @test norm(unchanged_A) ≈ norm_point + @test norm(zero_derivative) == 0 + scale!(unchanged_A, 0) + @test norm(point) ≈ norm_point + end + @test PEPSKit.norm_preserving_transport!(ζ, A, zerovector(A), α, A) === ζ end end From 557661cd90d63750160ac04e71ee358066c2e4fa Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 16 Oct 2026 12:25:47 +0800 Subject: [PATCH 14/14] Stabilize split CTMRG near product states and default to GMRES --- docs/src/man/environments.md | 6 +- src/algorithms/ctmrg/projectors/split.jl | 51 ++++++++--- .../optimization/implicit_differentiation.jl | 4 +- .../optimization/peps_optimization.jl | 2 +- src/algorithms/select_algorithm.jl | 15 +++- test/ctmrg/convergence.jl | 9 ++ test/ctmrg/split.jl | 14 ++- test/gradients/split_ctmrg_gradients.jl | 10 ++- test/testsuite/TestSuite.jl | 8 +- test/testsuite/ctmrg/convergence.jl | 61 +++++++++++++ test/testsuite/ctmrg/split.jl | 87 ++++++++++++++++++- 11 files changed, 241 insertions(+), 26 deletions(-) create mode 100644 test/ctmrg/convergence.jl create mode 100644 test/testsuite/ctmrg/convergence.jl diff --git a/docs/src/man/environments.md b/docs/src/man/environments.md index b447395b1..e75b7649f 100644 --- a/docs/src/man/environments.md +++ b/docs/src/man/environments.md @@ -23,8 +23,10 @@ Use `trunc` for the outer projectors and `interlayer_trunc` for the interlayer p When omitted, `interlayer_trunc` inherits `trunc`; the default `FixedSpaceTruncation()` preserves the existing spaces of each bond independently. `SplitCTMRGEnv(env::CTMRGEnv)` factorizes an existing environment exactly, while its `trunc` keyword can reduce the interlayer spaces. `CTMRGEnv(env::SplitCTMRGEnv)` contracts the interlayer bonds so the existing expectation-value routines can be reused. -This version supports PEPS overlap networks, full projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref). -The half-environment factors use direct square roots of their singular values; the projector pseudoinverse masks singular values below the relative cutoff `eps(real(scalartype(ψ)))^(3 / 4)`. +This version supports PEPS overlap networks, Fishman projectors, and automatic differentiation through finite sweeps or [`FixedPointGradient`](@ref). +The half-environment factors use direct square roots of their singular values, masking weights whose original singular values fall below the relative cutoff `eps(real(scalartype(ψ)))^(3 / 4)`; the projector pseudoinverse uses the same cutoff. + +GMRES is the default solver for `FixedPointGradient` when optimization uses `SplitCTMRG`; select another solver with `gradient_alg = (; solver_alg = (; alg = :Arnoldi))` or supply a gradient-algorithm instance. Fixed-point differentiation uses [`ScramblingEnvGauge`](@ref PEPSKit.ScramblingEnvGauge) on the alternating ket and bra boundary factors to align both outer and interlayer gauges. The ordinary and split environments use the same boundary-MPS gauge method, with each split edge contributing two MPS sites instead of one. diff --git a/src/algorithms/ctmrg/projectors/split.jl b/src/algorithms/ctmrg/projectors/split.jl index 01f2307b1..194c77e62 100644 --- a/src/algorithms/ctmrg/projectors/split.jl +++ b/src/algorithms/ctmrg/projectors/split.jl @@ -3,8 +3,9 @@ $(TYPEDEF) Projector algorithm for [`SplitCTMRG`](@ref), constructing the green (Fig. 3), teal (Fig. 7), yellow (Fig. 8), and red (Fig. 9) projector pairs of [Naumann et al.](@cite naumann_variationally_2025). Each projector pair follows [Fishman et al.](@cite fishman_faster_2018), Appendix A: SVD-factorize each environment half, keep its square-root-weighted factor facing the cut, and truncate the product of these factors. -The half-environment singular values are weighted by direct square roots; the final projector pseudoinverse uses the relative cutoff of [`sdiag_pow`](@ref). +The square-root weights of the half-environment singular values are masked using the relative cutoff of [`sdiag_pow`](@ref); the final projector pseudoinverse uses the same cutoff. QR compression is applied to the environment halves for the yellow projectors and, when the truncation permits, the quarter-environments for the red projectors before factorization. +SVD compression is used instead when a QR/LQ triangular factor is numerically rank deficient. ## Fields @@ -29,15 +30,45 @@ end PROJECTOR_SYMBOLS[:SplitProjector] = SplitProjector -"""Factor split environment halves into square-root-weighted Fishman factors facing their shared cut.""" -function _split_fishman_halves(Q::AbstractTensorMap, Q_next::AbstractTensorMap, alg::SplitProjector) +"""Use an untruncated dense decomposition for split environment factors.""" +function _split_half_svd_algorithm(alg::SplitProjector) half_alg = decomposition_algorithm(alg) half_alg.fwd_alg isa IterSVD && (half_alg = SVDAdjoint()) - half_alg = @set half_alg.fwd_alg = TruncatedAlgorithm(half_alg.fwd_alg, notrunc()) + return @set half_alg.fwd_alg = TruncatedAlgorithm(half_alg.fwd_alg, notrunc()) +end + +"""Compress a split environment with QR/LQ, falling back to SVD for a rank-deficient triangular factor.""" +function _split_compress_half(t::AbstractTensorMap, alg::SplitProjector, ::Val{left}) where {left} + factor = left ? last(left_orth(t, QRAdjoint())) : first(right_orth(t, LQAdjoint())) + fullrank = ignore_derivatives() do + all(blocks(factor)) do (_, b) + tolerance = MatrixAlgebraKit.default_pullback_rank_atol(b) + return all(d -> abs(d) > tolerance, diagview(b)) + end + end + fullrank && return factor + # Unpivoted QR/LQ can have internal zero pivots followed by nonzero pivots. + # Its triangular pullback solve is then singular even when the last pivot is nonzero. + U, S, V = svd_trunc_no_error(t, _split_half_svd_algorithm(alg)) + return left ? S * V : U * S +end + +"""Mask negligible half-environment square-root weights to suppress amplified roundoff.""" +function _split_fishman_sqrt(S::DiagonalTensorMap) + mask = ignore_derivatives() do + cutoff = eps(real(scalartype(S)))^(3 / 4) * norm(S, Inf) + DiagonalTensorMap(ifelse.(S.data .>= cutoff, one(scalartype(S)), zero(scalartype(S))), space(S, 1)) + end + return sdiag_pow(S, 0.5) * mask +end + +"""Factor split environment halves into square-root-weighted Fishman factors facing their shared cut.""" +function _split_fishman_halves(Q::AbstractTensorMap, Q_next::AbstractTensorMap, alg::SplitProjector) + half_alg = _split_half_svd_algorithm(alg) _, S, V = svd_trunc_no_error(Q, half_alg) U_next, S_next, _ = svd_trunc_no_error(Q_next, half_alg) - Q = sdiag_pow(S, 0.5) * V - Q_next = U_next * sdiag_pow(S_next, 0.5) + Q = _split_fishman_sqrt(S) * V + Q_next = U_next * _split_fishman_sqrt(S_next) return Q / norm(Q), Q_next / norm(Q_next) end @@ -51,8 +82,8 @@ function compute_projector( Q = Q / norm(Q) Q_next = Q_next / norm(Q_next) if compress - _, Q = left_orth(Q, QRAdjoint()) - Q_next, _ = right_orth(Q_next, LQAdjoint()) + Q = _split_compress_half(Q, alg, Val(true)) + Q_next = _split_compress_half(Q_next, alg, Val(false)) end # All four colors use Fishman factors before the final truncated SVD. @@ -174,8 +205,8 @@ function _split_red_projector( if reduce # SE = Qb * Rb and NE = Lt * Qt with isometric Qb and Qt. # Their factors retain the norms and nonzero spectrum of (SE * SW) * (NW * NE). - _, Rb = left_orth(SE, QRAdjoint()) - Lt, _ = right_orth(NE, LQAdjoint()) + Rb = _split_compress_half(SE, alg, Val(true)) + Lt = _split_compress_half(NE, alg, Val(false)) return compute_projector((Rb * SW, NW * Lt), alg; trunc) end # Preserve zero-spectrum dimensions and custom truncation ordering on the full halves. diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl index 3c79141fb..d12f82c6f 100644 --- a/src/algorithms/optimization/implicit_differentiation.jl +++ b/src/algorithms/optimization/implicit_differentiation.jl @@ -288,9 +288,9 @@ function _rrule( envinit isa SplitCTMRGEnv && scalartype(envinit) <: Real && throw(ArgumentError("Split fixed-point gradients require a complex environment; use complex(envinit).")) - env, = leading_boundary(envinit, state, alg) + env, info = leading_boundary(envinit, state, alg) - info, env_vjp = _fixedpoint_vjp(config, state, env, alg) + _, env_vjp = _fixedpoint_vjp(config, state, env, alg) function leading_boundary_fixed_pullback((Δenv′, Δinfo)) Δenv = unthunk(Δenv′) diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 84079c574..021dd2741 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -63,7 +63,7 @@ function PEPSOptimize(; reuse_env = Defaults.reuse_env, symmetrization = nothing, ) boundary_algorithm = _alg_or_nt(CTMRGAlgorithm, boundary_alg) - gradient_algorithm = _alg_or_nt(GradientAlgorithm, gradient_alg) + gradient_algorithm = _gradient_algorithm(gradient_alg, parent_alg(boundary_algorithm)) optimizer_algorithm = _alg_or_nt(OptimKit.OptimizationAlgorithm, optimizer_alg) precondition_algorithm = _alg_or_nt(PreconditionAlgorithm, precondition_alg) diff --git a/src/algorithms/select_algorithm.jl b/src/algorithms/select_algorithm.jl index 6324b71a6..15d6be7d5 100644 --- a/src/algorithms/select_algorithm.jl +++ b/src/algorithms/select_algorithm.jl @@ -14,6 +14,19 @@ other input, including `nothing`. parent_alg(alg) = alg parent_alg(alg::DynamicTol) = parent_alg(alg.alg) +"""Construct a gradient algorithm with split-CTMRG solver defaults, preserving explicit choices.""" +function _gradient_algorithm(gradient_alg, boundary_alg) + if boundary_alg isa SplitCTMRG && gradient_alg isa NamedTuple && + get(gradient_alg, :alg, Defaults.gradient_alg) === :FixedPointGradient + solver_alg = get(gradient_alg, :solver_alg, (;)) + if solver_alg isa NamedTuple + solver_alg = merge((; alg = :GMRES), solver_alg) + gradient_alg = merge(gradient_alg, (; solver_alg)) + end + end + return _alg_or_nt(GradientAlgorithm, gradient_alg) +end + """ _dynamic_tol_or_alg(alg; dynamic_tols::Bool, tol_min::Real, tol_max::Real, tol_factor::Real) @@ -101,7 +114,7 @@ function select_algorithm( ) gradient_kwargs = merge(defaults, gradient_alg) dynamic_tol_kwargs, gradient_kwargs = _pop_dynamic_tol_kwargs(gradient_kwargs) - gradient_alg = GradientAlgorithm(; gradient_kwargs...) + gradient_alg = _gradient_algorithm(gradient_kwargs, parent_alg(boundary_alg)) gradient_alg = _dynamic_tol_or_alg(gradient_alg; dynamic_tol_kwargs...) end diff --git a/test/ctmrg/convergence.jl b/test/ctmrg/convergence.jl new file mode 100644 index 000000000..9c81c83c3 --- /dev/null +++ b/test/ctmrg/convergence.jl @@ -0,0 +1,9 @@ +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_split_gradient_defaults(Vector) + TestSuite.ctmrg_fixedpoint_forward_info(Vector) +end diff --git a/test/ctmrg/split.jl b/test/ctmrg/split.jl index bc80b15e3..9bb5bb9de 100644 --- a/test/ctmrg/split.jl +++ b/test/ctmrg/split.jl @@ -1,7 +1,13 @@ @isdefined(TestSuite) || include("../testsuite/TestSuite.jl") using .TestSuite -TestSuite.ctmrg_split_projectors(Vector) -TestSuite.ctmrg_split_environment(Vector) -TestSuite.ctmrg_split_unitcell(Vector) -TestSuite.ctmrg_split_expectation_values(Vector) +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_split_projectors(Vector) + TestSuite.ctmrg_split_compression(Vector) + TestSuite.ctmrg_split_fishman_cutoff(Vector) + TestSuite.ctmrg_split_environment(Vector) + TestSuite.ctmrg_split_unitcell(Vector) + TestSuite.ctmrg_split_expectation_values(Vector) +end diff --git a/test/gradients/split_ctmrg_gradients.jl b/test/gradients/split_ctmrg_gradients.jl index 5f0bb2611..7a9605a7d 100644 --- a/test/gradients/split_ctmrg_gradients.jl +++ b/test/gradients/split_ctmrg_gradients.jl @@ -1,6 +1,10 @@ @isdefined(TestSuite) || include("../testsuite/TestSuite.jl") using .TestSuite -TestSuite.gradients_split_gauge(Vector) -TestSuite.gradients_split_sweep(Vector) -TestSuite.gradients_split_fixedpoint(Vector) +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.gradients_split_gauge(Vector) + TestSuite.gradients_split_sweep(Vector) + TestSuite.gradients_split_fixedpoint(Vector) +end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index cd210e393..5f82aa434 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -85,11 +85,17 @@ using .CompressLocal # ----- module CTMRGSplit include("ctmrg/split.jl") - export ctmrg_split_projectors + export ctmrg_split_projectors, ctmrg_split_compression, ctmrg_split_fishman_cutoff export ctmrg_split_environment, ctmrg_split_unitcell, ctmrg_split_expectation_values end using .CTMRGSplit +module CTMRGConvergence + include("ctmrg/convergence.jl") + export ctmrg_split_gradient_defaults, ctmrg_fixedpoint_forward_info +end +using .CTMRGConvergence + module CTMRGContractions include("ctmrg/contractions.jl") export ctmrg_contractions_random_cartesian_spaces, ctmrg_contractions_specific_u1_spaces diff --git a/test/testsuite/ctmrg/convergence.jl b/test/testsuite/ctmrg/convergence.jl new file mode 100644 index 000000000..503861a8d --- /dev/null +++ b/test/testsuite/ctmrg/convergence.jl @@ -0,0 +1,61 @@ +using Test, LinearAlgebra, KrylovKit, PEPSKit +using Random, TensorKit, Zygote, Adapt + +"""Check split-specific solver defaults across both optimizer entry points and explicit overrides.""" +function ctmrg_split_gradient_defaults(AT) + return @testset "Split fixed-point gradient defaults ($AT)" begin + Random.seed!(73418) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) + ordinary_env = CTMRGEnv(psi, ComplexSpace(2)) + split_env = SplitCTMRGEnv(ordinary_env) + constructors = ( + kwargs -> PEPSOptimize(; boundary_alg = SplitCTMRG(), kwargs...), + kwargs -> PEPSKit.select_algorithm(PEPSKit.fixedpoint, split_env; kwargs...), + ) + for construct in constructors + solver = PEPSKit.parent_alg(construct((;)).gradient_alg).solver_alg + @test solver isa GMRES + partial = (; solver_alg = (; tol = 1.0e-7, krylovdim = 17)) + solver = PEPSKit.parent_alg(construct((; gradient_alg = partial)).gradient_alg).solver_alg + @test solver isa GMRES + @test solver.tol == 1.0e-7 + @test solver.krylovdim == 17 + explicit = (; solver_alg = (; alg = :Arnoldi)) + solver = PEPSKit.parent_alg(construct((; gradient_alg = explicit)).gradient_alg).solver_alg + @test solver isa Arnoldi + requested = Arnoldi(; tol = 1.0e-7) + explicit = (; solver_alg = requested) + @test PEPSKit.parent_alg(construct((; gradient_alg = explicit)).gradient_alg).solver_alg === requested + explicit = FixedPointGradient(; solver_alg = Arnoldi(; tol = 1.0e-7)) + @test PEPSKit.parent_alg(construct((; gradient_alg = explicit)).gradient_alg) === explicit + @test construct((; gradient_alg = nothing)).gradient_alg === nothing + end + wrapped = PEPSKit.DynamicTol(SplitCTMRG(), 1.0e-12, 1.0e-4, 1.0e-3) + @test PEPSKit.parent_alg(PEPSOptimize(; boundary_alg = wrapped).gradient_alg).solver_alg isa GMRES + @test PEPSOptimize().gradient_alg.solver_alg isa Arnoldi + ordinary = PEPSKit.select_algorithm(PEPSKit.fixedpoint, ordinary_env) + @test PEPSKit.parent_alg(ordinary.gradient_alg).solver_alg isa Arnoldi + @test FixedPointGradient().solver_alg isa Arnoldi + end +end + +"""Preserve forward convergence diagnostics when attaching a fixed-point pullback.""" +function ctmrg_fixedpoint_forward_info(AT) + return @testset "Fixed-point forward convergence information ($AT)" begin + Random.seed!(73419) + psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(2))) + for alg in (SequentialCTMRG(; tol = 0.0, maxiter = 1, verbosity = 0), SplitCTMRG(; tol = 0.0, maxiter = 1, verbosity = 0)) + initial = CTMRGEnv(psi, ComplexSpace(2)) + env = alg isa SplitCTMRG ? SplitCTMRGEnv(initial) : initial + expected, expected_info = leading_boundary(env, psi, alg) + gradalg = FixedPointGradient(; solver_alg = (; alg = :GMRES)) + (actual, info), _ = Zygote.pullback(psi) do A + return PEPSKit.hook_pullback(leading_boundary, env, A, alg; alg_rrule = gradalg) + end + @test norm(actual - expected) ≤ 1.0e-12 * max(norm(expected), 1) + @test !info.converged + @test info.convergence_error ≈ expected_info.convergence_error + @test info.contraction_metrics == expected_info.contraction_metrics + end + end +end diff --git a/test/testsuite/ctmrg/split.jl b/test/testsuite/ctmrg/split.jl index 01f2c07df..df5cf2ab6 100644 --- a/test/testsuite/ctmrg/split.jl +++ b/test/testsuite/ctmrg/split.jl @@ -1,5 +1,5 @@ using Test, Random, LinearAlgebra -using TensorKit, MatrixAlgebraKit, PEPSKit, Adapt +using TensorKit, MatrixAlgebraKit, PEPSKit, Adapt, Zygote using PEPSKit: ctmrg_iteration, WEST """Construct a dense Fishman cut projector independently from fourth roots of the Gram matrices.""" @@ -44,13 +44,96 @@ function ctmrg_split_projectors(AT) spectrum = values support = [ones(min(retained, rank)); zeros(4 - min(retained, rank))] @test Array(info.S.data) ≈ normalize(spectrum)[1:rank] - @test count(!iszero, info.S.data) == min(count(!iszero, values), rank) + @test count(!iszero, info.S.data) == min(retained, rank) @test convert(Array, projectors[1] * projectors[2]) ≈ (norm(spectrum) / sum(spectrum)) * Diagonal(support) @test all(P -> isfinite(norm(P)), projectors) end end end +"""Check that Fishman square roots suppress subthreshold modes in both values and derivatives.""" +function ctmrg_split_fishman_cutoff(AT) + Random.seed!(37914) + return @testset "Fishman half-spectrum cutoff ($AT)" begin + V = ComplexSpace(4) + values = [1.0, 1.0e-6, 1.0e-15, 0.0] + expected = [1.0, 1.0e-3, 0.0, 0.0] + for scale in (1.0, 1.0e200, 1.0e-200) + S = adapt(AT, DiagonalTensorMap(scale * values, V)) + root = PEPSKit._split_fishman_sqrt(S) + @test Array(root.data) / sqrt(scale) ≈ expected + @test space(root) == space(S) + end + input = adapt(AT, values) + loss(x) = sum(abs2, PEPSKit._split_fishman_sqrt(DiagonalTensorMap(x, V)).data) + gradient = only(Zygote.gradient(loss, input)) + @test Array(gradient) ≈ [1.0, 1.0, 0.0, 0.0] + @test all(isfinite, gradient) + + # Tensor contractions send a dense cotangent into the diagonal square root. + S = adapt(AT, DiagonalTensorMap(values, V)) + weight = adapt(AT, randn(ComplexF64, V ← V)) + weight += weight' + gradient = only(Zygote.gradient(X -> real(dot(weight, PEPSKit._split_fishman_sqrt(X))), S)) + expected = real.(diag(convert(Array, weight))) .* [0.5, 500.0, 0.0, 0.0] + @test gradient ≈ adapt(AT, DiagonalTensorMap(expected, V)) + @test all(isfinite, gradient.data) + end +end + +"""Check split compression and its pullback when unpivoted QR/LQ has internal zero pivots.""" +function ctmrg_split_compression(AT) + Random.seed!(37915) + return @testset "Rank-deficient split compression ($AT)" begin + alg = SplitCTMRG(; verbosity = 0).projector_alg + for T in (Float64, ComplexF64), left in (true, false), deficient in (false, true) + dense = randn(T, 8, 6) + if deficient + dense[:, 2] .= 0 + dense[:, 4] .= 0 + end + left || (dense = Matrix(dense')) + A = adapt(AT, TensorMap(dense, ComplexSpace(size(dense, 1)) ← ComplexSpace(size(dense, 2)))) + factor = PEPSKit._split_compress_half(A, alg, Val(left)) + gram(X) = left ? X' * X : X * X' + @test gram(factor) ≈ gram(A) + if !deficient + reference = left ? last(left_orth(A, QRAdjoint())) : first(right_orth(A, LQAdjoint())) + @test factor == reference + end + V = left ? domain(A)[1] : codomain(A)[1] + H = adapt(AT, randn(T, V ← V)) + H += H' + loss(X) = real(dot(H, gram(PEPSKit._split_compress_half(X, alg, Val(left))))) + gradient = only(Zygote.gradient(loss, A)) + exact = left ? 2 * A * H : 2 * H * A + @test isfinite(norm(gradient)) + @test gradient ≈ exact atol = 1.0e-10 rtol = 1.0e-10 + direction = normalize(adapt(AT, randn(T, space(A)))) + epsilon = 1.0e-5 + finite = (loss(A + epsilon * direction) - loss(A - epsilon * direction)) / (2 * epsilon) + @test real(dot(gradient, direction)) ≈ finite atol = 1.0e-7 rtol = 1.0e-7 + if deficient + trunc = truncspace(ComplexSpace(3)) + direct, = PEPSKit.compute_projector((A, A'), alg; trunc) + compressed, = PEPSKit.compute_projector((A, A'), alg; trunc, compress = true) + @test compressed[1] * compressed[2] ≈ direct[1] * direct[2] + W = domain(A)[1] + weight = adapt(AT, randn(T, W ← W)) + weight += weight' + cut_loss = X -> begin + pair, = PEPSKit.compute_projector((X, X'), alg; trunc, compress = true) + return real(dot(weight, pair[1] * pair[2])) + end + gradient = only(Zygote.gradient(cut_loss, A)) + finite = (cut_loss(A + epsilon * direction) - cut_loss(A - epsilon * direction)) / (2 * epsilon) + @test isfinite(norm(gradient)) + @test real(dot(gradient, direction)) ≈ finite atol = 1.0e-6 rtol = 1.0e-5 + end + end + end +end + function ctmrg_split_environment(AT) return @testset "Split environment ($AT)" begin Random.seed!(47119)