From 13827ecc2cb800cb119e4adf4fa38cf04cd8a307 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 28 Nov 2024 18:31:31 +0100 Subject: [PATCH 01/19] Add HalfInfiniteProjector and FullInfiniteProjector, add beginnings of implementation --- src/PEPSKit.jl | 4 +- src/algorithms/ctmrg/ctmrg.jl | 102 +++++++++++++++++++++++++--------- 2 files changed, 79 insertions(+), 27 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index c8a6b032..388b6efe 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -59,9 +59,10 @@ include("utility/symmetrization.jl") const ctmrgscheme = :simultaneous const reuse_env = true const trscheme = FixedSpaceTruncation() - const fwd_alg = TensorKit.SVD() + const fwd_alg = TensorKit.SDD() const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1) const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg) + const projector_alg = HalfInfiniteProjector const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -84,6 +85,7 @@ Module containing default values that represent typical algorithm parameters. - `fwd_alg`: SVD algorithm that is used in the forward pass - `rrule_alg`: Reverse-rule for differentiating that SVD - `svd_alg`: Combination of `fwd_alg` and `rrule_alg` +- `projector_alg`: Algorithm to compute CTMRG projectors - `optimizer`: Optimization algorithm for PEPS ground-state optimization - `gradient_linsolver`: Default linear solver for the `LinSolver` gradient algorithm - `iterscheme`: Scheme for differentiating one CTMRG iteration diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 000e91e2..7ee2bc7c 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -8,31 +8,32 @@ have different spaces, this truncation style is different from `TruncationSpace` struct FixedSpaceTruncation <: TensorKit.TruncationScheme end """ - struct ProjectorAlg{S}(; svd_alg=TensorKit.SVD(), trscheme=TensorKit.notrunc(), - fixedspace=false, verbosity=0) - -Algorithm struct collecting all projector related parameters. The truncation scheme has to be -a `TensorKit.TruncationScheme`, and some SVD algorithms might have further restrictions on what -kind of truncation scheme can be used. If `fixedspace` is true, the truncation scheme is set to -`truncspace(V)` where `V` is the environment bond space, adjusted to the corresponding -environment direction/unit cell entry. + struct HalfInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the half-infinite CTMRG environment. """ -@kwdef struct ProjectorAlg{S<:SVDAdjoint,T} +@kwdef struct HalfInfiniteProjector{S<:SVDAdjoint,T} svd_alg::S = Defaults.svd_alg trscheme::T = Defaults.trscheme verbosity::Int = 0 end -# TODO: add option for different projector styles (half-infinite, full-infinite, etc.) -function truncation_scheme(alg::ProjectorAlg, edge) - if alg.trscheme isa FixedSpaceTruncation - return truncspace(space(edge, 1)) - else - return alg.trscheme - end +""" + struct FullInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG environment. +""" +@kwdef struct FullInfiniteProjector{S<:SVDAdjoint,T} + svd_alg::S = Defaults.svd_alg + trscheme::T = Defaults.trscheme + verbosity::Int = 0 end -function svd_algorithm(alg::ProjectorAlg, (dir, r, c)) +const ProjectorAlgs = Union{HalfInfiniteProjector,FullInfiniteProjector} + +function svd_algorithm(alg::ProjectorAlgs, (dir, r, c)) if alg.svd_alg isa SVDAdjoint{<:FixedSVD} fwd_alg = alg.svd_alg.fwd_alg fix_svd = FixedSVD(fwd_alg.U[dir, r, c], fwd_alg.S[dir, r, c], fwd_alg.V[dir, r, c]) @@ -42,6 +43,14 @@ function svd_algorithm(alg::ProjectorAlg, (dir, r, c)) end end +function truncation_scheme(alg::ProjectorAlgs, edge) + if alg.trscheme isa FixedSpaceTruncation + return truncspace(space(edge, 1)) + else + return alg.trscheme + end +end + """ CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, miniter=Defaults.ctmrg_miniter, verbosity=0, @@ -64,12 +73,12 @@ computed on the western side, and then applied and rotated. Or with `:simultaneo are computed and applied simultaneously on all sides, where in particular the corners get contracted with two projectors at the same time. """ -struct CTMRG{S} +struct CTMRG{S,P<:ProjectorAlgs} tol::Float64 maxiter::Int miniter::Int verbosity::Int - projector_alg::ProjectorAlg + projector_alg::P end function CTMRG(; tol=Defaults.ctmrg_tol, @@ -81,15 +90,19 @@ function CTMRG(; ctmrgscheme::Symbol=Defaults.ctmrgscheme, ) return CTMRG{ctmrgscheme}( - tol, maxiter, miniter, verbosity, ProjectorAlg(; svd_alg, trscheme, verbosity) + tol, + maxiter, + miniter, + verbosity, + Defaults.projector_alg(; svd_alg, trscheme, verbosity), ) end ctmrgscheme(::CTMRG{S}) where {S} = S # aliases for the different CTMRG schemes -const SequentialCTMRG = CTMRG{:sequential} -const SimultaneousCTMRG = CTMRG{:simultaneous} +const SequentialCTMRG{P} = CTMRG{:sequential,P} +const SimultaneousCTMRG{P} = CTMRG{:simultaneous,P} # supply correct constructor for Accessors.@set Accessors.constructorof(::Type{CTMRG{S}}) where {S} = CTMRG{S} @@ -139,7 +152,7 @@ end Perform one iteration of CTMRG that maps the `state` and `envs` to a new environment, and also returns the truncation error. """ -function ctmrg_iter(state, envs::CTMRGEnv, alg::SequentialCTMRG) +function ctmrg_iter(state, envs::CTMRGEnv, alg::SequentialCTMRG{}) ϵ = zero(real(scalartype(state))) for _ in 1:4 # left move @@ -155,7 +168,7 @@ function ctmrg_iter(state, envs::CTMRGEnv, alg::SequentialCTMRG) return envs, (; err=ϵ) end -function ctmrg_iter(state, envs::CTMRGEnv, alg::SimultaneousCTMRG) +function ctmrg_iter(state, envs::CTMRGEnv, alg::SimultaneousCTMRG{}) enlarged_envs = ctmrg_expand(state, envs, alg) projectors, info = ctmrg_projectors(enlarged_envs, envs, alg) envs′ = ctmrg_renormalize(enlarged_envs, projectors, state, envs, alg) @@ -206,7 +219,7 @@ Compute the CTMRG projectors based from enlarged environments. In the `:simultaneous` mode, the environment SVD is run in parallel. """ function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG + enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:HalfInfiniteProjector} ) where {C,E} projector_alg = alg.projector_alg ϵ = zero(real(scalartype(envs))) @@ -237,7 +250,44 @@ function ctmrg_projectors( return (map(first, projectors), map(last, projectors)), (; err=ϵ) end function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG + enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:FullInfiniteProjector} +) where {C,E} + projector_alg = alg.projector_alg + ϵ = zero(real(scalartype(envs))) + + coordinates = eachcoordinate(envs) + projectors = dtmap(coordinates) do (r, c) + # SVD half-infinite environment + r′ = _prev(r, size(envs.corners, 2)) + c′ = _next(c, size(envs.corners, 3)) + QQ = fullinfinite_environment( + enlarged_envs[1, r, c], + enlarged_envs[2, r′, c], + enlarged_envs[3, r′, c′], + enlarged_envs[4, r, c′], # TODO: break up the CTMRG expand/project/renormalize scheme! + ) + + trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) + svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) + U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) + ϵ = max(ϵ, ϵ_local / norm(S)) + + # Compute SVD truncation error and check for degenerate singular values + 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 + + # Compute projectors + return build_projectors(U, S, V, enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) + end + + return (map(first, projectors), map(last, projectors)), (; err=ϵ) +end +function ctmrg_projectors( + enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG{<:HalfInfiniteProjector} ) where {C,E} projector_alg = alg.projector_alg # pre-allocation From ccb366f70daeb45a17fcdbe221434c147dfde680 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 3 Dec 2024 17:50:55 +0100 Subject: [PATCH 02/19] Add half-infinite contractions methods and diagrams, refactor names slightly --- .../contractions/ctmrg_contractions.jl | 143 ++++++++++++++++-- src/algorithms/ctmrg/ctmrg.jl | 5 +- src/algorithms/ctmrg/sparse_environments.jl | 10 +- 3 files changed, 142 insertions(+), 16 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index e53bf2ac..1ef726a6 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -210,12 +210,12 @@ function right_projector(E_1, C, E_2, U, isqS, ket::PEPSTensor, bra::PEPSTensor= end """ - halfinfinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) - halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, + half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) + half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} - halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, + half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} - halfinfinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, + half_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} Contract two quadrants (enlarged corners) to form a half-infinite environment. @@ -236,7 +236,7 @@ The environment can also be contracted directly from all its constituent tensors | || || | ``` -Alternatively, contract environment with a vector `x` acting on it +Alternatively, contract the environment with a vector `x` acting on it ``` C_1 -- E_2 -- E_3 -- C_2 @@ -249,14 +249,14 @@ Alternatively, contract environment with a vector `x` acting on it or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. """ -function halfinfinite_environment( +function half_infinite_environment( quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3} ) where {S} return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := quadrant1[χ_in D_inabove D_inbelow; χ D1 D2] * quadrant2[χ D1 D2; χ_out D_outabove D_outbelow] end -function halfinfinite_environment( +function half_infinite_environment( C_1, C_2, E_1, E_2, E_3, E_4, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2 ) where {P<:PEPSTensor} return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := @@ -271,7 +271,7 @@ function halfinfinite_environment( C_2[χ4; χ5] * E_4[χ5 D7 D8; χ_out] end -function halfinfinite_environment( +function half_infinite_environment( C_1, C_2, E_1, @@ -297,7 +297,7 @@ function halfinfinite_environment( E_4[χ5 D7 D8; χ6] * x[χ6 D11 D12] end -function halfinfinite_environment( +function half_infinite_environment( x::AbstractTensor{S,3}, C_1, C_2, @@ -324,6 +324,131 @@ function halfinfinite_environment( conj(E_4[χ6 D9 D10; χ_in]) end +""" + +Contract four quadrants (enlarged corners) to form a full-infinite environment where the +open bonds can be on any side (diagrams show open bonds in the west). + +``` + |~~~~~~~~~| -- |~~~~~~~~~| + |quadrant1| |quadrant2| + |~~~~~~~~~| == |~~~~~~~~~| + | || || | + || | + | || || | + |~~~~~~~~~| -- |~~~~~~~~~| + |quadrant4| |quadrant3| + |~~~~~~~~~| == |~~~~~~~~~| +``` + +The environment can also be contracted directly from all its constituent tensors. + +``` + C_1 -- E_2 -- E_3 -- C_2 + | || || | + E_1 == ket_bra_1 == ket_bra_2 == E_4 + | || || | + || | + | || || | + E_8 == ket_bra_4 == ket_bra_3 == E_5 + | || || | + C_4 -- E_7 -- E_6 -- C_3 +``` + +Alternatively, contract the environment with a vector `x` acting on it + +``` + C_1 -- E_2 -- E_3 -- C_2 + | || || | + E_1 == ket_bra_1 == ket_bra_2 == E_4 + | || || | + || | + || | + | || || | + [~~~~x~~~~~] || | + | || || | + E_8 == ket_bra_4 == ket_bra_3 == E_5 + | || || | + C_4 -- E_7 -- E_6 -- C_3 + +``` +or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. +""" +function full_infinite_environment( + quadrant1::AbstractTensorMap{S,3,3}, + quadrant2::AbstractTensorMap{S,3,3}, + quadrant3::AbstractTensorMap{S,3,3}, + quadrant4::AbstractTensorMap{S,3,3}, +) where {S} end +function full_infinite_environment( + C_1, + C_2, + C_3, + C_4, + E_1, + E_2, + E_3, + E_4, + E_5, + E_6, + E_7, + E_8, + ket_1::P, + ket_2::P, + ket_3::P, + ket_4::P, + bra_1::P=ket_1, + bra_2::P=ket_2, + bra_3::P=ket_3, + bra_4::P=ket_4, +) where {P<:PEPSTensor} end +function full_infinite_environment( + C_1, + C_2, + C_3, + C_4, + E_1, + E_2, + E_3, + E_4, + E_5, + E_6, + E_7, + E_8, + x::AbstractTensor{S,3}, + ket_1::P, + ket_2::P, + ket_3::P, + ket_4::P, + bra_1::P=ket_1, + bra_2::P=ket_2, + bra_3::P=ket_3, + bra_4::P=ket_4, +) where {S,P<:PEPSTensor} end +function full_infinite_environment( + x::AbstractTensor{S,3}, + C_1, + C_2, + C_3, + C_4, + E_1, + E_2, + E_3, + E_4, + E_5, + E_6, + E_7, + E_8, + ket_1::P, + ket_2::P, + ket_3::P, + ket_4::P, + bra_1::P=ket_1, + bra_2::P=ket_2, + bra_3::P=ket_3, + bra_4::P=ket_4, +) where {S,P<:PEPSTensor} end + # Renormalization contractions # ---------------------------- diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 7ee2bc7c..ea4a0805 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -31,6 +31,7 @@ Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG envir verbosity::Int = 0 end +# TODO: do AbstractProjectorAlg type instead? -> would make it easier for users to implement custom projector alg const ProjectorAlgs = Union{HalfInfiniteProjector,FullInfiniteProjector} function svd_algorithm(alg::ProjectorAlgs, (dir, r, c)) @@ -228,7 +229,7 @@ function ctmrg_projectors( projectors = dtmap(coordinates) do (r, c) # SVD half-infinite environment r′ = _prev(r, size(envs.corners, 2)) - QQ = halfinfinite_environment(enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) + QQ = half_infinite_environment(enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) @@ -310,7 +311,7 @@ function ctmrg_projectors( end # SVD half-infinite environment - QQ = halfinfinite_environment( + QQ = half_infinite_environment( enlarged_envs[dir, r, c], enlarged_envs[_next(dir, 4), next_rc...] ) diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 9243067e..b6ca163b 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -139,7 +139,7 @@ end Instantiate half-infinite environment as `TensorMap` explicitly. """ function TensorKit.TensorMap(env::HalfInfiniteEnv) # Dense operator - return halfinfinite_environment( + return half_infinite_environment( env.C_1, env.C_2, env.E_1, @@ -161,7 +161,7 @@ Contract half-infinite environment with a vector `x`, such that the environment linear map or adjoint linear map on `x` if `Val(true)` or `Val(false)` is passed, respectively. """ function (env::HalfInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x - return halfinfinite_environment( + return half_infinite_environment( env.C_1, env.C_2, env.E_1, @@ -176,7 +176,7 @@ function (env::HalfInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x ) end function (env::HalfInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * x - return halfinfinite_environment( + return half_infinite_environment( x, env.C_1, env.C_2, @@ -191,8 +191,8 @@ function (env::HalfInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * ) end -# Wrapper around halfinfinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) -function halfinfinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) +# Wrapper around half_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) +function half_infinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) return HalfInfiniteEnv( ec_1.C, ec_2.C, From 582e3aeef75eeb92757056671ca3a8d90b36d820 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 3 Dec 2024 18:04:54 +0100 Subject: [PATCH 03/19] Add FullInfiniteEnviroment sparse struct --- .../contractions/ctmrg_contractions.jl | 7 +- src/algorithms/ctmrg/sparse_environments.jl | 198 ++++++++++++++++-- 2 files changed, 179 insertions(+), 26 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 1ef726a6..b7deba6b 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -375,11 +375,8 @@ Alternatively, contract the environment with a vector `x` acting on it or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. """ function full_infinite_environment( - quadrant1::AbstractTensorMap{S,3,3}, - quadrant2::AbstractTensorMap{S,3,3}, - quadrant3::AbstractTensorMap{S,3,3}, - quadrant4::AbstractTensorMap{S,3,3}, -) where {S} end + quadrant1::T, quadrant2::T, quadrant3::T, quadrant4::T +) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} end function full_infinite_environment( C_1, C_2, diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index b6ca163b..4926a7f7 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -97,9 +97,9 @@ function renormalize_southwest_corner(ec::EnlargedCorner, P_left, P_right) ) end -# ------------------ -# Sparse environment -# ------------------ +# -------------------------------- +# Sparse half-infinite environment +# -------------------------------- """ struct HalfInfiniteEnv{C,E,A,A′} @@ -114,8 +114,8 @@ struct HalfInfiniteEnv{C,E,A,A′} # TODO: subtype as AbstractTensorMap once Te E_3::E E_4::E ket_1::A - bra_1::A′ ket_2::A + bra_1::A′ bra_2::A′ end @@ -128,8 +128,10 @@ function HalfInfiniteEnv(quadrant1::EnlargedCorner, quadrant2::EnlargedCorner) quadrant1.E_2, quadrant2.E_1, quadrant2.E_2, - quadrant1.ket_bra, - quadrant2.ket_bra, + quadrant1.ket, + quadrant2.ket, + quadrant1.bra, + quadrant2.bra, ) end @@ -193,24 +195,10 @@ end # Wrapper around half_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) function half_infinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) - return HalfInfiniteEnv( - ec_1.C, - ec_2.C, - ec_1.E_1, - ec_1.E_2, - ec_2.E_1, - ec_2.E_2, - ec_1.ket, - ec_2.ket, - ec_1.bra, - ec_2.bra, - ) + return HalfInfiniteEnv(ec_1, ec_2) end -# ----------------------------------------------------- # AbstractTensorMap subtyping and IterSVD compatibility -# ----------------------------------------------------- - function TensorKit.domain(env::HalfInfiniteEnv) return domain(env.E_4) * domain(env.ket_2)[3] * domain(env.bra_2)[3]' end @@ -222,3 +210,171 @@ end function random_start_vector(env::HalfInfiniteEnv) return Tensor(randn, domain(env)) end + +# -------------------------------- +# Sparse full-infinite environment +# -------------------------------- + +""" + struct FullInfiniteEnv{C,E,A,A′} + +Full-infinite CTMRG environment tensor storage. +""" +struct FullInfiniteEnv{C,E,A,A′} # TODO: subtype as AbstractTensorMap once TensorKit is updated + C_1::C + C_2::C + C_3::C + C_4::C + E_1::E + E_2::E + E_3::E + E_4::E + E_5::E + E_6::E + E_7::E + E_8::E + ket_1::A + ket_2::A + ket_3::A + ket_4::A + bra_1::A′ + bra_2::A′ + bra_3::A′ + bra_4::A′ +end + +# Construct environment from two enlarged corners +function FullInfiniteEnv( + quadrant1::E, quadrant2::E, quadrant3::E, quadrant4::E +) where {E<:EnlargedCorner} + return FullInfiniteEnv( + quadrant1.C, + quadrant2.C, + quadrant3.C, + quadrant4.C, + quadrant1.E_1, + quadrant1.E_2, + quadrant2.E_1, + quadrant2.E_2, + quadrant3.E_1, + quadrant3.E_2, + quadrant4.E_1, + quadrant4.E_2, + quadrant1.ket, + quadrant2.ket, + quadrant3.ket, + quadrant4.ket, + quadrant1.bra, + quadrant2.bra, + quadrant3.bra, + quadrant4.bra, + ) +end + +""" + TensorKit.TensorMap(env::FullInfiniteEnv) + +Instantiate full-infinite environment as `TensorMap` explicitly. +""" +function TensorKit.TensorMap(env::FullInfiniteEnv) # Dense operator + return full_infinite_environment( + env.C_1, + env.C_2, + env.C_3, + env.C_4, + env.E_1, + env.E_2, + env.E_3, + env.E_4, + env.E_2, + env.E_3, + env.E_4, + env.E_5, + env.ket_1, + env.ket_2, + env.ket_3, + env.ket_4, + env.bra_1, + env.bra_2, + env.bra_3, + env.bra_4, + ) +end + +""" + (env::FullInfiniteEnv)(x, ::Val{false}) + (env::FullInfiniteEnv)(x, ::Val{true}) + +Contract full-infinite environment with a vector `x`, such that the environment acts as a +linear map or adjoint linear map on `x` if `Val(true)` or `Val(false)` is passed, respectively. +""" +function (env::FullInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x + return full_infinite_environment( + env.C_1, + env.C_2, + env.C_3, + env.C_4, + env.E_1, + env.E_2, + env.E_3, + env.E_4, + env.E_5, + env.E_6, + env.E_7, + env.E_8, + x, + env.ket_1, + env.ket_2, + env.ket_3, + env.ket_4, + env.bra_1, + env.bra_2, + env.bra_3, + env.bra_4, + ) +end +function (env::FullInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * x + return full_infinite_environment( + x, + env.C_1, + env.C_2, + env.C_3, + env.C_4, + env.E_1, + env.E_2, + env.E_3, + env.E_4, + env.E_5, + env.E_6, + env.E_7, + env.E_8, + env.ket_1, + env.ket_2, + env.ket_3, + env.ket_4, + env.bra_1, + env.bra_2, + env.bra_3, + env.bra_4, + ) +end + +# Wrapper around full_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) +function full_infinite_environment( + ec_1::E, ec_2::E, ec_3::E, ec_4::E +) where {E<:EnlargedCorner} + return FullInfiniteEnv(ec_1, ec_2, ec_3, ec_4) +end + +# AbstractTensorMap subtyping and IterSVD compatibility +function TensorKit.domain(env::FullInfiniteEnv) + return domain(env.E_8) * domain(env.ket_4)[3] * domain(env.bra_4)[3]' +end + +function TensorKit.codomain(env::FullInfiniteEnv) + return codomain(env.E_1)[1] * domain(env.ket_1)[3]' * domain(env.bra_1)[3] +end + +function random_start_vector(env::FullInfiniteEnv) + return Tensor(randn, domain(env)) +end From 1016c6b75a3a26052219dd5558bbdb563fcc0bf6 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 3 Dec 2024 18:15:39 +0100 Subject: [PATCH 04/19] Update full-inf docstrings --- .../contractions/ctmrg_contractions.jl | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index b7deba6b..17b7b667 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -212,11 +212,11 @@ end """ half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} half_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} Contract two quadrants (enlarged corners) to form a half-infinite environment. @@ -325,9 +325,19 @@ function half_infinite_environment( end """ + full_infinite_environment( + quadrant1::T, quadrant2::T, quadrant3::T, quadrant4::T + ) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} + full_infinite_environment(C_1, C_2, C_3, C_4, E_1, E_2, E_3, E_4, E_5, E_6, E_7, E_8, + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + full_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, + ket_1::P, ket_2::P, ket_3::P, ket_4::P, + bra_1::P=ket_1, bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4) where {P<:PEPSTensor} + full_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, + ket_1::P, ket_2::P, ket_3::P, ket_4::P, + bra_1::P=ket_1, bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4) where {P<:PEPSTensor} -Contract four quadrants (enlarged corners) to form a full-infinite environment where the -open bonds can be on any side (diagrams show open bonds in the west). +Contract four quadrants (enlarged corners) to form a full-infinite environment. ``` |~~~~~~~~~| -- |~~~~~~~~~| From 8c8eec11322f92a7b25a051e61c32511d3867b83 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Wed, 4 Dec 2024 10:50:34 +0100 Subject: [PATCH 05/19] Add actual full-infinite contractions --- .../contractions/ctmrg_contractions.jl | 86 +++++++++++++++++-- 1 file changed, 79 insertions(+), 7 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 17b7b667..61185f85 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -244,7 +244,6 @@ Alternatively, contract the environment with a vector `x` acting on it E_1 == ket_bra_1 == ket_bra_2 == E_4 | || || | [~~~~~~x~~~~~~] - || | ``` or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. @@ -310,7 +309,7 @@ function half_infinite_environment( bra_1::P=ket_1, bra_2::P=ket_2, ) where {S,P<:PEPSTensor} - return @autoopt @tensor env_x[χ_in D_inabove D_inbelow] := + return @autoopt @tensor x_env[χ_in D_inabove D_inbelow] := x[χ1 D1 D2] * conj(E_1[χ1 D3 D4; χ2]) * conj(C_1[χ2; χ3]) * @@ -374,7 +373,6 @@ Alternatively, contract the environment with a vector `x` acting on it | || || | || | || | - | || || | [~~~~x~~~~~] || | | || || | E_8 == ket_bra_4 == ket_bra_3 == E_5 @@ -386,7 +384,13 @@ or contract the adjoint environment with `x`, e.g. as needed for iterative solve """ function full_infinite_environment( quadrant1::T, quadrant2::T, quadrant3::T, quadrant4::T -) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} end +) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} + return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := + quadrant1[χ_in D_inabove D_inbelow; χ1 D1 D2] * + quadrant2[χ1 D1 D2; χ2 D3 D4] * + quadrant3[χ2 D3 D4; χ3 D5 D6] * + quadrant4[χ3 D5 D6; χ_out D_outabove D_outbelow] +end function full_infinite_environment( C_1, C_2, @@ -408,7 +412,29 @@ function full_infinite_environment( bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4, -) where {P<:PEPSTensor} end +) where {P<:PEPSTensor} + return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := + E_1[χ_in D1 D2; χ1] * + C_1[χ1; χ2] * + E_2[χ2 D3 D4; χ3] * + ket_1[d1; D3 D11 D_inabove D1] * + conj(bra_1[d1; D4 D12 D_inbelow D2]) * + ket_2[d2; D5 D7 D9 D11] * + conj(bra_2[d2; D6 D8 D10 D12]) * + E_3[χ3 D5 D6; χ4] * + C_2[χ4; χ5] * + E_4[χ5 D7 D8; χ6] * + E_5[χ6 D13 D14; χ7] * + C_3[χ7; χ8] * + E_6[χ8 D15 D16; χ9] * + ket_3[d3; D9 D13 D15 D17] * + conj(bra_3[d3; D10 D14 D16 D18]) * + ket_4[d4; D_outabove D17 D19 D21] * + conj(bra_4[d4; D_outbelow D18 D20 D22]) * + E_7[χ8 D19 D20; χ9] * + C_4[χ9; χ10] * + E_8[χ10 D21 D22; χ_out] +end function full_infinite_environment( C_1, C_2, @@ -431,7 +457,30 @@ function full_infinite_environment( bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4, -) where {S,P<:PEPSTensor} end +) where {S,P<:PEPSTensor} + return @autoopt @tensor env_x[χ_in D_inabove D_inbelow] := + E_1[χ_in D1 D2; χ1] * + C_1[χ1; χ2] * + E_2[χ2 D3 D4; χ3] * + ket_1[d1; D3 D11 D_inabove D1] * + conj(bra_1[d1; D4 D12 D_inbelow D2]) * + ket_2[d2; D5 D7 D9 D11] * + conj(bra_2[d2; D6 D8 D10 D12]) * + E_3[χ3 D5 D6; χ4] * + C_2[χ4; χ5] * + E_4[χ5 D7 D8; χ6] * + E_5[χ6 D13 D14; χ7] * + C_3[χ7; χ8] * + E_6[χ8 D15 D16; χ9] * + ket_3[d3; D9 D13 D15 D17] * + conj(bra_3[d3; D10 D14 D16 D18]) * + ket_4[d4; D_xabove D17 D19 D21] * + conj(bra_4[d4; D_xbelow D18 D20 D22]) * + E_7[χ8 D19 D20; χ9] * + C_4[χ9; χ10] * + E_8[χ10 D21 D22; χ_x] * + x[χ_x D_xabove D_xbelow] +end function full_infinite_environment( x::AbstractTensor{S,3}, C_1, @@ -454,7 +503,30 @@ function full_infinite_environment( bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4, -) where {S,P<:PEPSTensor} end +) where {S,P<:PEPSTensor} + return @autoopt @tensor x_env[χ_in D_inabove D_inbelow] := + x[χ_x D_xabove D_xbelow] * + E_1[χ_x D1 D2; χ1] * + C_1[χ1; χ2] * + E_2[χ2 D3 D4; χ3] * + ket_1[d1; D3 D11 D_xabove D1] * + conj(bra_1[d1; D4 D12 D_xbelow D2]) * + ket_2[d2; D5 D7 D9 D11] * + conj(bra_2[d2; D6 D8 D10 D12]) * + E_3[χ3 D5 D6; χ4] * + C_2[χ4; χ5] * + E_4[χ5 D7 D8; χ6] * + E_5[χ6 D13 D14; χ7] * + C_3[χ7; χ8] * + E_6[χ8 D15 D16; χ9] * + ket_3[d3; D9 D13 D15 D17] * + conj(bra_3[d3; D10 D14 D16 D18]) * + ket_4[d4; D_inabove D17 D19 D21] * + conj(bra_4[d4; D_inbelow D18 D20 D22]) * + E_7[χ8 D19 D20; χ9] * + C_4[χ9; χ10] * + E_8[χ10 D21 D22; χ_in] +end # Renormalization contractions # ---------------------------- From 4d0db6e0fcb9bb3db3141761c3c8047c7ceb1bf6 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Wed, 4 Dec 2024 18:23:20 +0100 Subject: [PATCH 06/19] Add FullInfiniteProjector projector computation and remove redundant previous code --- .../contractions/ctmrg_contractions.jl | 205 -------------- src/algorithms/ctmrg/ctmrg.jl | 258 +++++++++--------- src/algorithms/ctmrg/sparse_environments.jl | 168 ------------ 3 files changed, 133 insertions(+), 498 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 61185f85..2c95af6f 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -323,211 +323,6 @@ function half_infinite_environment( conj(E_4[χ6 D9 D10; χ_in]) end -""" - full_infinite_environment( - quadrant1::T, quadrant2::T, quadrant3::T, quadrant4::T - ) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} - full_infinite_environment(C_1, C_2, C_3, C_4, E_1, E_2, E_3, E_4, E_5, E_6, E_7, E_8, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} - full_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, - ket_1::P, ket_2::P, ket_3::P, ket_4::P, - bra_1::P=ket_1, bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4) where {P<:PEPSTensor} - full_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, - ket_1::P, ket_2::P, ket_3::P, ket_4::P, - bra_1::P=ket_1, bra_2::P=ket_2, bra_3::P=ket_3, bra_4::P=ket_4) where {P<:PEPSTensor} - -Contract four quadrants (enlarged corners) to form a full-infinite environment. - -``` - |~~~~~~~~~| -- |~~~~~~~~~| - |quadrant1| |quadrant2| - |~~~~~~~~~| == |~~~~~~~~~| - | || || | - || | - | || || | - |~~~~~~~~~| -- |~~~~~~~~~| - |quadrant4| |quadrant3| - |~~~~~~~~~| == |~~~~~~~~~| -``` - -The environment can also be contracted directly from all its constituent tensors. - -``` - C_1 -- E_2 -- E_3 -- C_2 - | || || | - E_1 == ket_bra_1 == ket_bra_2 == E_4 - | || || | - || | - | || || | - E_8 == ket_bra_4 == ket_bra_3 == E_5 - | || || | - C_4 -- E_7 -- E_6 -- C_3 -``` - -Alternatively, contract the environment with a vector `x` acting on it - -``` - C_1 -- E_2 -- E_3 -- C_2 - | || || | - E_1 == ket_bra_1 == ket_bra_2 == E_4 - | || || | - || | - || | - [~~~~x~~~~~] || | - | || || | - E_8 == ket_bra_4 == ket_bra_3 == E_5 - | || || | - C_4 -- E_7 -- E_6 -- C_3 - -``` -or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. -""" -function full_infinite_environment( - quadrant1::T, quadrant2::T, quadrant3::T, quadrant4::T -) where {T<:AbstractTensorMap{<:ElementarySpace,3,3}} - return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := - quadrant1[χ_in D_inabove D_inbelow; χ1 D1 D2] * - quadrant2[χ1 D1 D2; χ2 D3 D4] * - quadrant3[χ2 D3 D4; χ3 D5 D6] * - quadrant4[χ3 D5 D6; χ_out D_outabove D_outbelow] -end -function full_infinite_environment( - C_1, - C_2, - C_3, - C_4, - E_1, - E_2, - E_3, - E_4, - E_5, - E_6, - E_7, - E_8, - ket_1::P, - ket_2::P, - ket_3::P, - ket_4::P, - bra_1::P=ket_1, - bra_2::P=ket_2, - bra_3::P=ket_3, - bra_4::P=ket_4, -) where {P<:PEPSTensor} - return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := - E_1[χ_in D1 D2; χ1] * - C_1[χ1; χ2] * - E_2[χ2 D3 D4; χ3] * - ket_1[d1; D3 D11 D_inabove D1] * - conj(bra_1[d1; D4 D12 D_inbelow D2]) * - ket_2[d2; D5 D7 D9 D11] * - conj(bra_2[d2; D6 D8 D10 D12]) * - E_3[χ3 D5 D6; χ4] * - C_2[χ4; χ5] * - E_4[χ5 D7 D8; χ6] * - E_5[χ6 D13 D14; χ7] * - C_3[χ7; χ8] * - E_6[χ8 D15 D16; χ9] * - ket_3[d3; D9 D13 D15 D17] * - conj(bra_3[d3; D10 D14 D16 D18]) * - ket_4[d4; D_outabove D17 D19 D21] * - conj(bra_4[d4; D_outbelow D18 D20 D22]) * - E_7[χ8 D19 D20; χ9] * - C_4[χ9; χ10] * - E_8[χ10 D21 D22; χ_out] -end -function full_infinite_environment( - C_1, - C_2, - C_3, - C_4, - E_1, - E_2, - E_3, - E_4, - E_5, - E_6, - E_7, - E_8, - x::AbstractTensor{S,3}, - ket_1::P, - ket_2::P, - ket_3::P, - ket_4::P, - bra_1::P=ket_1, - bra_2::P=ket_2, - bra_3::P=ket_3, - bra_4::P=ket_4, -) where {S,P<:PEPSTensor} - return @autoopt @tensor env_x[χ_in D_inabove D_inbelow] := - E_1[χ_in D1 D2; χ1] * - C_1[χ1; χ2] * - E_2[χ2 D3 D4; χ3] * - ket_1[d1; D3 D11 D_inabove D1] * - conj(bra_1[d1; D4 D12 D_inbelow D2]) * - ket_2[d2; D5 D7 D9 D11] * - conj(bra_2[d2; D6 D8 D10 D12]) * - E_3[χ3 D5 D6; χ4] * - C_2[χ4; χ5] * - E_4[χ5 D7 D8; χ6] * - E_5[χ6 D13 D14; χ7] * - C_3[χ7; χ8] * - E_6[χ8 D15 D16; χ9] * - ket_3[d3; D9 D13 D15 D17] * - conj(bra_3[d3; D10 D14 D16 D18]) * - ket_4[d4; D_xabove D17 D19 D21] * - conj(bra_4[d4; D_xbelow D18 D20 D22]) * - E_7[χ8 D19 D20; χ9] * - C_4[χ9; χ10] * - E_8[χ10 D21 D22; χ_x] * - x[χ_x D_xabove D_xbelow] -end -function full_infinite_environment( - x::AbstractTensor{S,3}, - C_1, - C_2, - C_3, - C_4, - E_1, - E_2, - E_3, - E_4, - E_5, - E_6, - E_7, - E_8, - ket_1::P, - ket_2::P, - ket_3::P, - ket_4::P, - bra_1::P=ket_1, - bra_2::P=ket_2, - bra_3::P=ket_3, - bra_4::P=ket_4, -) where {S,P<:PEPSTensor} - return @autoopt @tensor x_env[χ_in D_inabove D_inbelow] := - x[χ_x D_xabove D_xbelow] * - E_1[χ_x D1 D2; χ1] * - C_1[χ1; χ2] * - E_2[χ2 D3 D4; χ3] * - ket_1[d1; D3 D11 D_xabove D1] * - conj(bra_1[d1; D4 D12 D_xbelow D2]) * - ket_2[d2; D5 D7 D9 D11] * - conj(bra_2[d2; D6 D8 D10 D12]) * - E_3[χ3 D5 D6; χ4] * - C_2[χ4; χ5] * - E_4[χ5 D7 D8; χ6] * - E_5[χ6 D13 D14; χ7] * - C_3[χ7; χ8] * - E_6[χ8 D15 D16; χ9] * - ket_3[d3; D9 D13 D15 D17] * - conj(bra_3[d3; D10 D14 D16 D18]) * - ket_4[d4; D_inabove D17 D19 D21] * - conj(bra_4[d4; D_inbelow D18 D20 D22]) * - E_7[χ8 D19 D20; χ9] * - C_4[χ9; χ10] * - E_8[χ10 D21 D22; χ_in] -end - # Renormalization contractions # ---------------------------- diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index ea4a0805..8d6f6642 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -44,9 +44,9 @@ function svd_algorithm(alg::ProjectorAlgs, (dir, r, c)) end end -function truncation_scheme(alg::ProjectorAlgs, edge) +function truncation_scheme(alg::ProjectorAlgs, Espace) if alg.trscheme isa FixedSpaceTruncation - return truncspace(space(edge, 1)) + return truncspace(Espace) else return alg.trscheme end @@ -213,139 +213,147 @@ end # Projector step # ======================================================================================== # -""" - ctmrg_projectors(enlarged_envs, env, alg::CTMRG{M}) - -Compute the CTMRG projectors based from enlarged environments. -In the `:simultaneous` mode, the environment SVD is run in parallel. -""" -function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:HalfInfiniteProjector} -) where {C,E} - projector_alg = alg.projector_alg - ϵ = zero(real(scalartype(envs))) - - coordinates = eachcoordinate(envs) - projectors = dtmap(coordinates) do (r, c) - # SVD half-infinite environment - r′ = _prev(r, size(envs.corners, 2)) - QQ = half_infinite_environment(enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) - - trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) - svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) - U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - ϵ = max(ϵ, ϵ_local / norm(S)) - - # Compute SVD truncation error and check for degenerate singular values - 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 - - # Compute projectors - return build_projectors(U, S, V, enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) - end - - return (map(first, projectors), map(last, projectors)), (; err=ϵ) -end -function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:FullInfiniteProjector} -) where {C,E} - projector_alg = alg.projector_alg - ϵ = zero(real(scalartype(envs))) - - coordinates = eachcoordinate(envs) - projectors = dtmap(coordinates) do (r, c) - # SVD half-infinite environment - r′ = _prev(r, size(envs.corners, 2)) - c′ = _next(c, size(envs.corners, 3)) - QQ = fullinfinite_environment( - enlarged_envs[1, r, c], - enlarged_envs[2, r′, c], - enlarged_envs[3, r′, c′], - enlarged_envs[4, r, c′], # TODO: break up the CTMRG expand/project/renormalize scheme! - ) - - trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) - svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) - U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - ϵ = max(ϵ, ϵ_local / norm(S)) - - # Compute SVD truncation error and check for degenerate singular values - 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 +# TODO: embed this into new functions after refactor +function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) + # SVD half-infinite environment + halfinf = half_infinite_environment(enlarged_corners...) + trscheme = truncation_scheme(projector_alg, space(enlarged_corners[2], 1)) + svd_alg = svd_algorithm(projector_alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=trscheme) + + # Compute SVD truncation error and check for degenerate singular values + 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 - - # Compute projectors - return build_projectors(U, S, V, enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) end - return (map(first, projectors), map(last, projectors)), (; err=ϵ) + P_left, P_right = build_projectors(U, S, V, enlarged_corners...) + return (P_left, P_right), (; err, U, S, V) end -function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG{<:HalfInfiniteProjector} -) where {C,E} - projector_alg = alg.projector_alg - # pre-allocation - U, V = Zygote.Buffer.(projector_type(envs.edges)) - # Corner type but with real numbers - S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) - - ϵ = zero(real(scalartype(envs))) - coordinates = eachcoordinate(envs, 1:4) - projectors = dtmap(coordinates) do (dir, r, c) - # Row-column index of next enlarged corner - next_rc = if dir == 1 - (r, _next(c, size(envs.corners, 3))) - elseif dir == 2 - (_next(r, size(envs.corners, 2)), c) - elseif dir == 3 - (r, _prev(c, size(envs.corners, 3))) - elseif dir == 4 - (_prev(r, size(envs.corners, 2)), c) - end - # SVD half-infinite environment - QQ = half_infinite_environment( - enlarged_envs[dir, r, c], enlarged_envs[_next(dir, 4), next_rc...] - ) - - trscheme = truncation_scheme(projector_alg, envs.edges[dir, next_rc...]) - svd_alg = svd_algorithm(projector_alg, (dir, r, c)) - U_local, S_local, V_local, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - U[dir, r, c] = U_local - S[dir, r, c] = S_local - V[dir, r, c] = V_local - ϵ = max(ϵ, ϵ_local / norm(S_local)) - - # Compute SVD truncation error and check for degenerate singular values - Zygote.isderiving() && ignore_derivatives() do - if alg.verbosity > 0 && is_degenerate_spectrum(S_local) - svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S_local)) - @warn("degenerate singular values detected: ", svals) - end +function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) + # QR top and bottom half-infinite environments + halfinf_top = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) + halfinf_bot = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) + _, R_top = leftorth!(halfinf_top) + _, R_bot = leftorth!(halfinf_bot) + R_fullinf = R_top * R_bot + + # SVD product of R's + trscheme = truncation_scheme(projector_alg, space(enlarged_corners[4], 1)) + svd_alg = svd_algorithm(projector_alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(R_fullinf, svd_alg; trunc=trscheme) + + # Compute SVD truncation error and check for degenerate singular values + 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 - - # Compute projectors - return build_projectors( - U_local, - S_local, - V_local, - enlarged_envs[dir, r, c], - enlarged_envs[_next(dir, 4), next_rc...], - ) end - P_left = map(first, projectors) - P_right = map(last, projectors) - return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) + P_left, P_right = build_projectors(U, S, V, R_top, R_bot) + return (P_left, P_right), (; err, U, S, V) end +# """ +# ctmrg_projectors(enlarged_envs, env, alg::CTMRG{M}) + +# Compute the CTMRG projectors based from enlarged environments. +# In the `:simultaneous` mode, the environment SVD is run in parallel. +# """ +# function ctmrg_projectors( +# enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:HalfInfiniteProjector} +# ) where {C,E} +# projector_alg = alg.projector_alg +# ϵ = zero(real(scalartype(envs))) + +# coordinates = eachcoordinate(envs) +# projectors = dtmap(coordinates) do (r, c) +# # SVD half-infinite environment +# r′ = _prev(r, size(envs.corners, 2)) +# QQ = half_infinite_environment(enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) + +# trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) +# svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) +# U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) +# ϵ = max(ϵ, ϵ_local / norm(S)) + +# # Compute SVD truncation error and check for degenerate singular values +# 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 + +# # Compute projectors +# return build_projectors(U, S, V, enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) +# end + +# return (map(first, projectors), map(last, projectors)), (; err=ϵ) +# end +# function ctmrg_projectors( +# enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG{<:HalfInfiniteProjector} +# ) where {C,E} +# projector_alg = alg.projector_alg +# # pre-allocation +# U, V = Zygote.Buffer.(projector_type(envs.edges)) +# # Corner type but with real numbers +# S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) + +# ϵ = zero(real(scalartype(envs))) +# coordinates = eachcoordinate(envs, 1:4) +# projectors = dtmap(coordinates) do (dir, r, c) +# # Row-column index of next enlarged corner +# next_rc = if dir == 1 +# (r, _next(c, size(envs.corners, 3))) +# elseif dir == 2 +# (_next(r, size(envs.corners, 2)), c) +# elseif dir == 3 +# (r, _prev(c, size(envs.corners, 3))) +# elseif dir == 4 +# (_prev(r, size(envs.corners, 2)), c) +# end + +# # SVD half-infinite environment +# QQ = half_infinite_environment( +# enlarged_envs[dir, r, c], enlarged_envs[_next(dir, 4), next_rc...] +# ) + +# trscheme = truncation_scheme(projector_alg, envs.edges[dir, next_rc...]) +# svd_alg = svd_algorithm(projector_alg, (dir, r, c)) +# U_local, S_local, V_local, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) +# U[dir, r, c] = U_local +# S[dir, r, c] = S_local +# V[dir, r, c] = V_local +# ϵ = max(ϵ, ϵ_local / norm(S_local)) + +# # Compute SVD truncation error and check for degenerate singular values +# Zygote.isderiving() && ignore_derivatives() do +# if alg.verbosity > 0 && is_degenerate_spectrum(S_local) +# svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S_local)) +# @warn("degenerate singular values detected: ", svals) +# end +# end + +# # Compute projectors +# return build_projectors( +# U_local, +# S_local, +# V_local, +# enlarged_envs[dir, r, c], +# enlarged_envs[_next(dir, 4), next_rc...], +# ) +# end + +# P_left = map(first, projectors) +# P_right = map(last, projectors) +# return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) +# end + """ build_projectors(U::AbstractTensorMap{E,3,1}, S::AbstractTensorMap{E,1,1}, V::AbstractTensorMap{E,1,3}, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3}) where {E<:ElementarySpace} diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 4926a7f7..4f9fbad2 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -210,171 +210,3 @@ end function random_start_vector(env::HalfInfiniteEnv) return Tensor(randn, domain(env)) end - -# -------------------------------- -# Sparse full-infinite environment -# -------------------------------- - -""" - struct FullInfiniteEnv{C,E,A,A′} - -Full-infinite CTMRG environment tensor storage. -""" -struct FullInfiniteEnv{C,E,A,A′} # TODO: subtype as AbstractTensorMap once TensorKit is updated - C_1::C - C_2::C - C_3::C - C_4::C - E_1::E - E_2::E - E_3::E - E_4::E - E_5::E - E_6::E - E_7::E - E_8::E - ket_1::A - ket_2::A - ket_3::A - ket_4::A - bra_1::A′ - bra_2::A′ - bra_3::A′ - bra_4::A′ -end - -# Construct environment from two enlarged corners -function FullInfiniteEnv( - quadrant1::E, quadrant2::E, quadrant3::E, quadrant4::E -) where {E<:EnlargedCorner} - return FullInfiniteEnv( - quadrant1.C, - quadrant2.C, - quadrant3.C, - quadrant4.C, - quadrant1.E_1, - quadrant1.E_2, - quadrant2.E_1, - quadrant2.E_2, - quadrant3.E_1, - quadrant3.E_2, - quadrant4.E_1, - quadrant4.E_2, - quadrant1.ket, - quadrant2.ket, - quadrant3.ket, - quadrant4.ket, - quadrant1.bra, - quadrant2.bra, - quadrant3.bra, - quadrant4.bra, - ) -end - -""" - TensorKit.TensorMap(env::FullInfiniteEnv) - -Instantiate full-infinite environment as `TensorMap` explicitly. -""" -function TensorKit.TensorMap(env::FullInfiniteEnv) # Dense operator - return full_infinite_environment( - env.C_1, - env.C_2, - env.C_3, - env.C_4, - env.E_1, - env.E_2, - env.E_3, - env.E_4, - env.E_2, - env.E_3, - env.E_4, - env.E_5, - env.ket_1, - env.ket_2, - env.ket_3, - env.ket_4, - env.bra_1, - env.bra_2, - env.bra_3, - env.bra_4, - ) -end - -""" - (env::FullInfiniteEnv)(x, ::Val{false}) - (env::FullInfiniteEnv)(x, ::Val{true}) - -Contract full-infinite environment with a vector `x`, such that the environment acts as a -linear map or adjoint linear map on `x` if `Val(true)` or `Val(false)` is passed, respectively. -""" -function (env::FullInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x - return full_infinite_environment( - env.C_1, - env.C_2, - env.C_3, - env.C_4, - env.E_1, - env.E_2, - env.E_3, - env.E_4, - env.E_5, - env.E_6, - env.E_7, - env.E_8, - x, - env.ket_1, - env.ket_2, - env.ket_3, - env.ket_4, - env.bra_1, - env.bra_2, - env.bra_3, - env.bra_4, - ) -end -function (env::FullInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * x - return full_infinite_environment( - x, - env.C_1, - env.C_2, - env.C_3, - env.C_4, - env.E_1, - env.E_2, - env.E_3, - env.E_4, - env.E_5, - env.E_6, - env.E_7, - env.E_8, - env.ket_1, - env.ket_2, - env.ket_3, - env.ket_4, - env.bra_1, - env.bra_2, - env.bra_3, - env.bra_4, - ) -end - -# Wrapper around full_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) -function full_infinite_environment( - ec_1::E, ec_2::E, ec_3::E, ec_4::E -) where {E<:EnlargedCorner} - return FullInfiniteEnv(ec_1, ec_2, ec_3, ec_4) -end - -# AbstractTensorMap subtyping and IterSVD compatibility -function TensorKit.domain(env::FullInfiniteEnv) - return domain(env.E_8) * domain(env.ket_4)[3] * domain(env.bra_4)[3]' -end - -function TensorKit.codomain(env::FullInfiniteEnv) - return codomain(env.E_1)[1] * domain(env.ket_1)[3]' * domain(env.bra_1)[3] -end - -function random_start_vector(env::FullInfiniteEnv) - return Tensor(randn, domain(env)) -end From 3ea633eb2d829b617beba3fddda3efbb871047f4 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 5 Dec 2024 17:40:37 +0100 Subject: [PATCH 07/19] Split up CTMRG file and refactor simultaneous and sequential modes, incorporate new projectors algs --- src/PEPSKit.jl | 6 +- .../contractions/ctmrg_contractions.jl | 25 ++ src/algorithms/ctmrg/ctmrg.jl | 335 ++---------------- src/algorithms/ctmrg/sequential.jl | 85 +++++ src/algorithms/ctmrg/simultaneous.jl | 81 +++++ src/algorithms/peps_opt.jl | 2 +- src/environments/ctmrg_environments.jl | 33 +- src/utility/util.jl | 26 +- test/ctmrg/{ctmrgschemes.jl => flavors.jl} | 12 +- test/ctmrg/gaugefix.jl | 19 +- test/ctmrg/gradients.jl | 3 +- test/ctmrg/gradparts.jl | 5 +- test/runtests.jl | 2 +- 13 files changed, 284 insertions(+), 350 deletions(-) create mode 100644 src/algorithms/ctmrg/sequential.jl create mode 100644 src/algorithms/ctmrg/simultaneous.jl rename test/ctmrg/{ctmrgschemes.jl => flavors.jl} (88%) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 388b6efe..a1875fc1 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -54,9 +54,9 @@ include("utility/symmetrization.jl") const ctmrg_maxiter = 100 const ctmrg_miniter = 4 const ctmrg_tol = 1e-8 + const ctmrg_flavor = :simultaneous const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 - const ctmrgscheme = :simultaneous const reuse_env = true const trscheme = FixedSpaceTruncation() const fwd_alg = TensorKit.SDD() @@ -77,9 +77,9 @@ Module containing default values that represent typical algorithm parameters. - `ctmrg_maxiter`: Maximal number of CTMRG iterations per run - `ctmrg_miniter`: Minimal number of CTMRG carried out - `ctmrg_tol`: Tolerance checking singular value and norm convergence +- `ctmrg_flavor`: Scheme for growing, projecting and renormalizing CTMRG environments - `fpgrad_maxiter`: Maximal number of iterations for computing the CTMRG fixed-point gradient - `fpgrad_tol`: Convergence tolerance for the fixed-point gradient iteration -- `ctmrgscheme`: Scheme for growing, projecting and renormalizing CTMRG environments - `reuse_env`: If `true`, the current optimization step is initialized on the previous environment - `trscheme`: Truncation scheme for SVDs and other decompositions - `fwd_alg`: SVD algorithm that is used in the forward pass @@ -98,9 +98,9 @@ module Defaults const ctmrg_maxiter = 100 const ctmrg_miniter = 4 const ctmrg_tol = 1e-8 + const ctmrg_flavor = :simultaneous const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 - const ctmrgscheme = :simultaneous const sparse = false const reuse_env = true const trscheme = FixedSpaceTruncation() diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index c014265d..0a84ae06 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -163,6 +163,31 @@ end # Projector contractions # ---------------------- +""" + left_and_right_projector(U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} + left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner) + +Compute left and right projectors based on a SVD and quadrant tensors, specified either as +`AbstractTensorMap`s or sparsely as `EnlargedCorner`s such that the quadrants are never +constructed explicitly. +""" +function left_and_right_projector( + U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} +) where {E<:ElementarySpace} + isqS = sdiag_inv_sqrt(S) + P_left = Q_next * V' * isqS + P_right = isqS * U' * Q + return P_left, P_right +end +function left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner) + isqS = sdiag_inv_sqrt(S) + P_left = left_projector(Q.E_1, Q.C, Q.E_2, V, isqS, Q.ket, Q.bra) + P_right = right_projector( + Q_next.E_1, Q_next.C, Q_next.E_2, U, isqS, Q_next.ket, Q_next.bra + ) + return P_left, P_right +end + """ left_projector(E_1, C, E_2, V, isqS, ket::PEPSTensor, bra::PEPSTensor=ket) diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index bc478870..445ee14e 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -54,60 +54,54 @@ end """ CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, verbosity=0, - svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation(), - ctmrgscheme=Defaults.ctmrgscheme) + miniter=Defaults.ctmrg_miniter, flavor=Defaults.ctmrg_flavor, verbosity=0, + svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) Algorithm struct that represents the CTMRG algorithm for contracting infinite PEPS. Each CTMRG run is converged up to `tol` where the singular value convergence of the corners as well as the norm is checked. The maximal and minimal number of CTMRG iterations -is set with `maxiter` and `miniter`. Different levels of output information are printed -depending on `verbosity`, where `0` suppresses all output, `1` only prints warnings, `2` -gives information at the start and end, and `3` prints information every iteration. +is set with `maxiter` and `miniter`. + +In general, two different flavors of CTMRG can be selected with `flavor` which determine how +CTMRG is implemented. It can either be `:sequential`, where the projectors are succesively +computed on the west side, and then applied and rotated. Or with `:simultaneous` all projectors +are computed and applied simultaneously on all sides, where the corners get contracted with +two projectors at the same time. + +Different levels of output information are printed depending on `verbosity`, where `0` +suppresses all output, `1` only prints warnings, `2` gives information at the start and +end, and `3` prints information every iteration. The projectors are computed from `svd_alg` SVDs where the truncation scheme is set via `trscheme`. - -In general, two different schemes can be selected with `ctmrgscheme` which determine how -CTMRG is implemented. It can either be `:sequential`, where the projectors are succesively -computed on the western side, and then applied and rotated. Or with `:simultaneous` all projectors -are computed and applied simultaneously on all sides, where in particular the corners get -contracted with two projectors at the same time. """ -struct CTMRG{S,P<:ProjectorAlgs} +struct CTMRG tol::Float64 maxiter::Int miniter::Int + flavor::Symbol verbosity::Int - projector_alg::P + projector_alg::ProjectorAlgs end function CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, miniter=Defaults.ctmrg_miniter, + flavor=Defaults.ctmrg_flavor, verbosity=2, svd_alg=Defaults.svd_alg, trscheme=Defaults.trscheme, - ctmrgscheme::Symbol=Defaults.ctmrgscheme, ) - return CTMRG{ctmrgscheme}( + return CTMRG( tol, maxiter, miniter, + flavor, verbosity, Defaults.projector_alg(; svd_alg, trscheme, verbosity), ) end -ctmrgscheme(::CTMRG{S}) where {S} = S - -# aliases for the different CTMRG schemes -const SequentialCTMRG{P} = CTMRG{:sequential,P} -const SimultaneousCTMRG{P} = CTMRG{:simultaneous,P} - -# supply correct constructor for Accessors.@set -Accessors.constructorof(::Type{CTMRG{S}}) where {S} = CTMRG{S} - """ MPSKit.leading_boundary([envinit], state, alg::CTMRG) @@ -126,10 +120,15 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) env = deepcopy(envinit) log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) + f = if alg.flavor == :sequential + sequential_ctmrg_iter + elseif alg.flavor == :simultaneous + simultaneous_ctmrg_iter + end return LoggingExtras.withlevel(; alg.verbosity) do ctmrg_loginit!(log, η, N) for iter in 1:(alg.maxiter) - env, = ctmrg_iter(state, env, alg) # Grow and renormalize in all 4 directions + env, = f(state, env, alg) # Grow and renormalize in all 4 directions η, CS, TS = calc_convergence(env, CS, TS) N = norm(state, env) @@ -147,33 +146,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) end end -""" - ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info - -Perform one iteration of CTMRG that maps the `state` and `envs` to a new environment, -and also returns the `info` `NamedTuple`. -""" -function ctmrg_iter(state, envs::CTMRGEnv, alg::SequentialCTMRG{}) - ϵ = zero(real(scalartype(state))) - for _ in 1:4 # rotate - for col in 1:size(state, 2) # left move column-wise - projectors, info = ctmrg_projectors(col, state, envs, alg) - envs = ctmrg_renormalize(col, projectors, state, envs, alg) - ϵ = max(ϵ, info.err) - end - state = rotate_north(state, EAST) - envs = rotate_north(envs, EAST) - end - - return envs, (; err=ϵ) -end -function ctmrg_iter(state, envs::CTMRGEnv, alg::SimultaneousCTMRG) - enlarged_envs = ctmrg_expand(eachcoordinate(state, 1:4), state, envs) - projectors, info = ctmrg_projectors(enlarged_envs, envs, alg) - envs′ = ctmrg_renormalize(enlarged_envs, projectors, state, envs, alg) - return envs′, info -end - +# custom CTMRG logging ctmrg_loginit!(log, η, N) = @infov 2 loginit!(log, η, N) ctmrg_logiter!(log, iter, η, N) = @infov 3 logiter!(log, iter, η, N) ctmrg_logfinish!(log, iter, η, N) = @infov 2 logfinish!(log, iter, η, N) @@ -184,24 +157,12 @@ ctmrg_logcancel!(log, iter, η, N) = @warnv 1 logcancel!(log, iter, η, N) @non_differentiable ctmrg_logfinish!(args...) @non_differentiable ctmrg_logcancel!(args...) -# ======================================================================================== # -# Expansion step -# ======================================================================================== # - """ - ctmrg_expand(coordinates, state, envs) + compute_projector(enlarged_corners, coordinate, alg::ProjectorAlgs) -Expand the environment by absorbing a new PEPS tensor on the given coordinates. +Determine left and right projectors at the bond given determined by the enlarged corners +and the given coordinate using the specified `alg`. """ -function ctmrg_expand(coordinates, state, envs::CTMRGEnv) - return dtmap(idx -> TensorMap(EnlargedCorner(state, envs, idx), idx[1]), coordinates) -end - -# ======================================================================================== # -# Projector step -# ======================================================================================== # - -# TODO: embed this into new functions after refactor function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) # SVD half-infinite environment halfinf = half_infinite_environment(enlarged_corners...) @@ -217,50 +178,9 @@ function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjec end end - P_left, P_right = build_projectors(U, S, V, enlarged_corners...) + P_left, P_right = left_and_right_projector(U, S, V, enlarged_corners...) return (P_left, P_right), (; err, U, S, V) end - -""" - ctmrg_projectors(col::Int, enlarged_envs, env, alg::CTMRG{:sequential}) - ctmrg_projectors(enlarged_envs, env, alg::CTMRG{:simultaneous}) - -Compute the CTMRG projectors based on enlarged environments. -In the `:sequential` mode the projectors are computed for the column `col`, whereas -in the `:simultaneous` mode, all projectors (and corresponding SVDs) are computed in parallel. -""" -function ctmrg_projectors( - col::Int, state::InfinitePEPS, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG -) where {C,E} - projector_alg = alg.projector_alg - ϵ = zero(real(scalartype(envs))) - - # SVD half-infinite environment - coordinates = eachcoordinate(envs)[:, col] - projectors = dtmap(coordinates) do (r, c) - r′ = _prev(r, size(envs.corners, 2)) - Q1 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r, c)), SOUTHWEST) - Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST) - QQ = halfinfinite_environment(Q1, Q2) - trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) - svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) - U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - ϵ = max(ϵ, ϵ_local / norm(S)) - - # Compute SVD truncation error and check for degenerate singular values - 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 - - # Compute projectors - return build_projectors(U, S, V, Q1, Q2) - end - return (map(first, projectors), map(last, projectors)), (; err=ϵ) -end - function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) # QR top and bottom half-infinite environments halfinf_top = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) @@ -282,199 +202,6 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec end end - P_left, P_right = build_projectors(U, S, V, R_top, R_bot) + P_left, P_right = left_and_right_projector(U, S, V, R_top, R_bot) return (P_left, P_right), (; err, U, S, V) end - -# """ -# ctmrg_projectors(enlarged_envs, env, alg::CTMRG{M}) - -# Compute the CTMRG projectors based from enlarged environments. -# In the `:simultaneous` mode, the environment SVD is run in parallel. -# """ -# function ctmrg_projectors( -# enlarged_envs, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG{<:HalfInfiniteProjector} -# ) where {C,E} -# projector_alg = alg.projector_alg -# ϵ = zero(real(scalartype(envs))) - -# coordinates = eachcoordinate(envs) -# projectors = dtmap(coordinates) do (r, c) -# # SVD half-infinite environment -# r′ = _prev(r, size(envs.corners, 2)) -# QQ = half_infinite_environment(enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) - -# trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) -# svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) -# U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) -# ϵ = max(ϵ, ϵ_local / norm(S)) - -# # Compute SVD truncation error and check for degenerate singular values -# 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 - -# # Compute projectors -# return build_projectors(U, S, V, enlarged_envs[1, r, c], enlarged_envs[2, r′, c]) -# end - -# return (map(first, projectors), map(last, projectors)), (; err=ϵ) -# end -# function ctmrg_projectors( -# enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG{<:HalfInfiniteProjector} -# ) where {C,E} -# projector_alg = alg.projector_alg -# # pre-allocation -# U, V = Zygote.Buffer.(projector_type(envs.edges)) -# # Corner type but with real numbers -# S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) - -# ϵ = zero(real(scalartype(envs))) -# coordinates = eachcoordinate(envs, 1:4) -# projectors = dtmap(coordinates) do (dir, r, c) -# # Row-column index of next enlarged corner -# next_rc = if dir == 1 -# (r, _next(c, size(envs.corners, 3))) -# elseif dir == 2 -# (_next(r, size(envs.corners, 2)), c) -# elseif dir == 3 -# (r, _prev(c, size(envs.corners, 3))) -# elseif dir == 4 -# (_prev(r, size(envs.corners, 2)), c) -# end - -# # SVD half-infinite environment -# QQ = half_infinite_environment( -# enlarged_envs[dir, r, c], enlarged_envs[_next(dir, 4), next_rc...] -# ) - -# trscheme = truncation_scheme(projector_alg, envs.edges[dir, next_rc...]) -# svd_alg = svd_algorithm(projector_alg, (dir, r, c)) -# U_local, S_local, V_local, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) -# U[dir, r, c] = U_local -# S[dir, r, c] = S_local -# V[dir, r, c] = V_local -# ϵ = max(ϵ, ϵ_local / norm(S_local)) - -# # Compute SVD truncation error and check for degenerate singular values -# Zygote.isderiving() && ignore_derivatives() do -# if alg.verbosity > 0 && is_degenerate_spectrum(S_local) -# svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S_local)) -# @warn("degenerate singular values detected: ", svals) -# end -# end - -# # Compute projectors -# return build_projectors( -# U_local, -# S_local, -# V_local, -# enlarged_envs[dir, r, c], -# enlarged_envs[_next(dir, 4), next_rc...], -# ) -# end - -# P_left = map(first, projectors) -# P_right = map(last, projectors) -# return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) -# end - -""" - build_projectors(U::AbstractTensorMap{E,3,1}, S::AbstractTensorMap{E,1,1}, V::AbstractTensorMap{E,1,3}, - Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3}) where {E<:ElementarySpace} - build_projectors(U::AbstractTensorMap{E,3,1}, S::AbstractTensorMap{E,1,1}, V::AbstractTensorMap{E,1,3}, - Q::EnlargedCorner, Q_next::EnlargedCorner) where {E<:ElementarySpace} - -Construct left and right projectors where the higher-dimensional is facing left and right, respectively. -""" -function build_projectors( - U::AbstractTensorMap{E,3,1}, - S::AbstractTensorMap{E,1,1}, - V::AbstractTensorMap{E,1,3}, - Q::AbstractTensorMap{E,3,3}, - Q_next::AbstractTensorMap{E,3,3}, -) where {E<:ElementarySpace} - isqS = sdiag_inv_sqrt(S) - P_left = Q_next * V' * isqS - P_right = isqS * U' * Q - return P_left, P_right -end -function build_projectors( - U::AbstractTensorMap{E,3,1}, - S::AbstractTensorMap{E,1,1}, - V::AbstractTensorMap{E,1,3}, - Q::EnlargedCorner, - Q_next::EnlargedCorner, -) where {E<:ElementarySpace} - isqS = sdiag_inv_sqrt(S) - P_left = left_projector(Q.E_1, Q.C, Q.E_2, V, isqS, Q.ket, Q.bra) - P_right = right_projector( - Q_next.E_1, Q_next.C, Q_next.E_2, U, isqS, Q_next.ket, Q_next.bra - ) - return P_left, P_right -end - -# ======================================================================================== # -# Renormalization step -# ======================================================================================== # - -""" - ctmrg_renormalize(col::Int, projectors, state, envs, ::CTMRG{:sequential}) - ctmrg_renormalize(enlarged_envs, projectors, state, envs, ::CTMRG{:simultaneous}) - -Apply projectors to renormalize corners and edges. -The `:sequential` mode renormalizes the environment on the column `col`, where as the -`:simultaneous` mode renormalizes all environment tensors simultaneously. -""" -function ctmrg_renormalize(col::Int, projectors, state, envs, ::SequentialCTMRG) - corners = Zygote.Buffer(envs.corners) - edges = Zygote.Buffer(envs.edges) - - for (dir, r, c) in eachcoordinate(state, 1:4) - (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue - corners[dir, r, c] = envs.corners[dir, r, c] - end - for (dir, r, c) in eachcoordinate(state, 1:4) - (c == col && dir == WEST) && continue - edges[dir, r, c] = envs.edges[dir, r, c] - end - - # Apply projectors to renormalize corners and edge - for row in axes(envs.corners, 2) - C_southwest = renormalize_bottom_corner((row, col), envs, projectors) - corners[SOUTHWEST, row, col] = C_southwest / norm(C_southwest) - - C_northwest = renormalize_top_corner((row, col), envs, projectors) - corners[NORTHWEST, row, col] = C_northwest / norm(C_northwest) - - E_west = renormalize_west_edge((row, col), envs, projectors, state) - edges[WEST, row, col] = E_west / norm(E_west) - end - - return CTMRGEnv(copy(corners), copy(edges)) -end -function ctmrg_renormalize(enlarged_envs, projectors, state, envs, ::SimultaneousCTMRG) - P_left, P_right = projectors - coordinates = eachcoordinate(envs, 1:4) - corners_edges = dtmap(coordinates) do (dir, r, c) - if dir == NORTH - corner = renormalize_northwest_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_north_edge((r, c), envs, P_left, P_right, state) - elseif dir == EAST - corner = renormalize_northeast_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_east_edge((r, c), envs, P_left, P_right, state) - elseif dir == SOUTH - corner = renormalize_southeast_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_south_edge((r, c), envs, P_left, P_right, state) - elseif dir == WEST - corner = renormalize_southwest_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_west_edge((r, c), envs, P_left, P_right, state) - end - return corner / norm(corner), edge / norm(edge) - end - - return CTMRGEnv(map(first, corners_edges), map(last, corners_edges)) -end diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl new file mode 100644 index 00000000..7631fefe --- /dev/null +++ b/src/algorithms/ctmrg/sequential.jl @@ -0,0 +1,85 @@ +""" + sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info + +Perform one sequential iteration of CTMRG, where one iteration consists of four expansion, +renormalization and rotation steps that are performed sequentially. +""" +function sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) + ϵ = zero(real(scalartype(state))) + for _ in 1:4 # rotate + for col in 1:size(state, 2) # left move column-wise + projectors, info = sequential_projectors(col, state, envs, alg.projector_alg) + envs = renormalize_sequentially(col, projectors, state, envs) + ϵ = max(ϵ, info.err) + end + state = rotate_north(state, EAST) + envs = rotate_north(envs, EAST) + end + + return envs, (; err=ϵ) +end + +function sequential_projectors( + col::Int, state::InfinitePEPS, envs::CTMRGEnv{C,E}, alg::ProjectorAlgs +) where {C,E} + ϵ = zero(real(scalartype(envs))) + + # SVD half-infinite environment column-wise + coordinates = eachcoordinate(envs)[:, col] + projectors = dtmap(coordinates) do (r, c) + proj, info = simultaneous_projector(state, envs, (WEST, r, c), alg) + ϵ = max(ϵ, info.err / norm(info.S)) + return proj + end + + return (map(first, projectors), map(last, projectors)), (; err=ϵ) +end +function sequential_projector( + state::InfinitePEPS, envs::CTMRGEnv, coordinate, alg::HalfInfiniteProjector +) + _, r, c = coordinate + r′ = _prev(r, size(envs, 2)) + Q1 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r, c)), SOUTHWEST) + Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST) + return compute_projector((Q1, Q2), coordinate, alg) +end +function sequential_projector( + state::InfinitePEPS, envs::CTMRGEnv, coordinate, alg::FullInfiniteProjector +) + _, r, c = coordinate + r′ = _next(r, size(envs, 2)) + c′ = _next(c, size(envs, 3)) + Q1 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r, c)), NORTHWEST) + Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHEAST, r, c′)), NORTHEAST) + Q3 = TensorMap(EnlargedCorner(state, envs, (SOUTHEAST, r′, c′)), SOUTHEAST) + Q4 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r′, c)), SOUTHWEST) + return compute_projector((Q1, Q2, Q3, Q4), coordinate, alg) +end + +function renormalize_sequentially(col::Int, projectors, state, envs) + corners = Zygote.Buffer(envs.corners) + edges = Zygote.Buffer(envs.edges) + + for (dir, r, c) in eachcoordinate(state, 1:4) + (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue + corners[dir, r, c] = envs.corners[dir, r, c] + end + for (dir, r, c) in eachcoordinate(state, 1:4) + (c == col && dir == WEST) && continue + edges[dir, r, c] = envs.edges[dir, r, c] + end + + # Apply projectors to renormalize corners and edge + for row in axes(envs.corners, 2) + C_southwest = renormalize_bottom_corner((row, col), envs, projectors) + corners[SOUTHWEST, row, col] = C_southwest / norm(C_southwest) + + C_northwest = renormalize_top_corner((row, col), envs, projectors) + corners[NORTHWEST, row, col] = C_northwest / norm(C_northwest) + + E_west = renormalize_west_edge((row, col), envs, projectors, state) + edges[WEST, row, col] = E_west / norm(E_west) + end + + return CTMRGEnv(copy(corners), copy(edges)) +end diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl new file mode 100644 index 00000000..43f34fc9 --- /dev/null +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -0,0 +1,81 @@ +""" + simultaneous_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info + +Perform one simultaneous iteration of CTMRG, in which the environment is expanded and +renormalized in all directions at the same time. +""" +function simultaneous_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) + enlarged_corners = dtmap(eachcoordinate(state, 1:4)) do idx + return TensorMap(EnlargedCorner(state, envs, idx), idx[1]) + end # expand environment + projectors, info = simultaneous_projectors(enlarged_corners, envs, alg.projector_alg) # compute projectors on all coordinates + envs′ = renormalize_simultaneously(enlarged_corners, projectors, state, envs) # renormalize enlarged corners + return envs′, info +end + +function simultaneous_projectors( + enlarged_corners, envs::CTMRGEnv{C,E}, alg::ProjectorAlgs +) where {C,E} + # pre-allocation + U, V = Zygote.Buffer.(projector_type(envs.edges)) + S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) # Corner type but with real numbers + ϵ = zero(real(scalartype(envs))) + + projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate + proj, info = compute_projector(enlarged_corners, coordinate, alg) + U[dir, r, c] = info.U + S[dir, r, c] = info.S + V[dir, r, c] = info.V + ϵ = max(ϵ, info.err / norm(info.S)) + return proj + end + + P_left = map(first, projectors) + P_right = map(last, projectors) + return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) +end +function simultaneous_projector( + enlarged_corners::Array{E,3}, coordinate, alg::HalfInfiniteProjector +) where {E} + coordinate′ = _next_projector_coordinate(coordinate, size(enlarged_corners)[2:3]...) + ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) + return compute_projector(ec, coordinate, alg) +end +function simultaneous_projector( + enlarged_corners::Array{E,3}, coordinate, alg::FullInfiniteProjector +) where {E} + rowsize, colsize = size(enlarged_corners)[2:3] + coordinate2 = _next_projector_coordinate(coordinate, rowsize, colsize) + coordinate3 = _next_projector_coordinate(coordinate2, rowsize, colsize) + coordinate4 = _next_projector_coordinate(coordinate3, rowsize, colsize) + ec = ( + enlarged_corners[coordinate...], + enlarged_corners[coordinate2...], + enlarged_corners[coordinate3...], + enlarged_corners[coordinate4...], + ) + return compute_projector(ec, coordinate, alg) +end + +function renormalize_simultaneously(enlarged_corners, projectors, state, envs) + P_left, P_right = projectors + coordinates = eachcoordinate(envs, 1:4) + corners_edges = dtmap(coordinates) do (dir, r, c) + if dir == NORTH + corner = renormalize_northwest_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_north_edge((r, c), envs, P_left, P_right, state) + elseif dir == EAST + corner = renormalize_northeast_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_east_edge((r, c), envs, P_left, P_right, state) + elseif dir == SOUTH + corner = renormalize_southeast_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_south_edge((r, c), envs, P_left, P_right, state) + elseif dir == WEST + corner = renormalize_southwest_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_west_edge((r, c), envs, P_left, P_right, state) + end + return corner / norm(corner), edge / norm(edge) + end + + return CTMRGEnv(map(first, corners_edges), map(last, corners_edges)) +end diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 39e138a1..00eac45f 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -265,7 +265,7 @@ function _rrule( fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) alg_fixed = CTMRG(; - svd_alg=svd_alg_fixed, trscheme=notrunc(), ctmrgscheme=:simultaneous + svd_alg=svd_alg_fixed, trscheme=notrunc(), flavor=:simultaneous ) function leading_boundary_fixed_pullback(Δenvs′) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index c026762b..29a823cc 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -354,23 +354,7 @@ function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnv, name::Symbol) throw(ArgumentError("No rrule for getproperty of $name")) end end - -# Rotate corners & edges counter-clockwise -function Base.rotl90(env::CTMRGEnv{C,T}) where {C,T} - # Initialize rotated corners & edges with rotated sizes - corners′ = Zygote.Buffer( - Array{C,3}(undef, 4, size(env.corners, 3), size(env.corners, 2)) - ) - edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env.edges, 3), size(env.edges, 2))) - - for dir in 1:4 - corners′[_prev(dir, 4), :, :] = rotl90(env.corners[dir, :, :]) - edges′[_prev(dir, 4), :, :] = rotl90(env.edges[dir, :, :]) - end - - return CTMRGEnv(copy(corners′), copy(edges′)) -end - +Base.size(env::CTMRG, args...) = size(env.corners, args...) Base.eltype(env::CTMRGEnv) = eltype(env.corners[1]) Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) function eachcoordinate(x::CTMRGEnv) @@ -379,7 +363,6 @@ end function eachcoordinate(x::CTMRGEnv, dirs) return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) end - Base.real(env::CTMRGEnv) = CTMRGEnv(real.(env.corners), real.(env.edges)) Base.complex(env::CTMRGEnv) = CTMRGEnv(complex.(env.corners), complex.(env.edges)) @@ -390,6 +373,20 @@ function update!(env::CTMRGEnv{C,T}, env´::CTMRGEnv{C,T}) where {C,T} return env end +# Rotate corners & edges counter-clockwise +function Base.rotl90(env::CTMRGEnv{C,T}) where {C,T} + # Initialize rotated corners & edges with rotated sizes + corners′ = Zygote.Buffer(Array{C,3}(undef, 4, size(env, 3), size(env, 2))) + edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env, 3), size(env, 2))) + + for dir in 1:4 + corners′[_prev(dir, 4), :, :] = rotl90(env.corners[dir, :, :]) + edges′[_prev(dir, 4), :, :] = rotl90(env.edges[dir, :, :]) + end + + return CTMRGEnv(copy(corners′), copy(edges′)) +end + # Functions used for FP differentiation and by KrylovKit.linsolve function Base.:+(e₁::CTMRGEnv, e₂::CTMRGEnv) return CTMRGEnv(e₁.corners + e₂.corners, e₁.edges + e₂.edges) diff --git a/src/utility/util.jl b/src/utility/util.jl index 0b9fd4c4..93cb9577 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -1,7 +1,31 @@ -# Get next and previous directional CTM enviroment index, respecting periodicity +# Get next and previous directional CTMRG environment index, respecting periodicity _next(i, total) = mod1(i + 1, total) _prev(i, total) = mod1(i - 1, total) +# Get next and previous coordinate (direction, row, column), given a direction and going around the environment clockwise +function _next_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_next(dir, 4), row, _next(col, colsize)) + elseif dir == 2 + return (_next(dir, 4), _next(row, rowsize), col) + elseif dir == 3 + return (_next(dir, 4), row, _prev(col, colsize)) + elseif dir == 4 + return (_next(dir, 4), _prev(row, rowsize), col) + end +end +function _prev_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_prev(dir, 4), _next(row, rowsize), col) + elseif dir == 2 + return (_prev(dir, 4), row, _prev(col, colsize)) + elseif dir == 3 + return (_prev(dir, 4), _prev(row, rowsize), col) + elseif dir == 4 + return (_prev(dir, 4), row, _next(col, colsize)) + end +end + # iterator over each coordinates """ eachcoordinate(x, dirs=1:4) diff --git a/test/ctmrg/ctmrgschemes.jl b/test/ctmrg/flavors.jl similarity index 88% rename from test/ctmrg/ctmrgschemes.jl rename to test/ctmrg/flavors.jl index 406d4c1f..c1c848a9 100644 --- a/test/ctmrg/ctmrgschemes.jl +++ b/test/ctmrg/flavors.jl @@ -7,8 +7,8 @@ using PEPSKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg_sequential = CTMRG(; ctmrgscheme=:sequential) -ctm_alg_simultaneous = CTMRG(; ctmrgscheme=:simultaneous) +ctm_alg_sequential = CTMRG(; flavor=:sequential) +ctm_alg_simultaneous = CTMRG(; flavor=:simultaneous) unitcells = [(1, 1), (3, 4)] @testset "$(unitcell) unit cell" for unitcell in unitcells @@ -52,13 +52,9 @@ unitcells = [(1, 1), (3, 4)] end # test fixedspace actually fixes space -@testset "Fixedspace truncation ($scheme)" for scheme in [:sequential, :simultaneous] +@testset "Fixedspace truncation ($flavor)" for flavor in [:sequential, :simultaneous] ctm_alg = CTMRG(; - tol=1e-6, - maxiter=1, - verbosity=0, - ctmrgscheme=scheme, - trscheme=FixedSpaceTruncation(), + tol=1e-6, maxiter=1, verbosity=0, flavor, trscheme=FixedSpaceTruncation() ) Ds = fill(2, 3, 3) χs = [16 17 18; 15 20 21; 14 19 22] diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index ef36943c..748c5b05 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -8,7 +8,7 @@ using PEPSKit: ctmrg_iter, gauge_fix, calc_elementwise_convergence scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] maxiter = 200 -schemes = [:simultaneous, :sequential] +ctmrg_flavors = [:simultaneous, :sequential] χ = 6 atol = 1e-4 @@ -34,10 +34,9 @@ function _semi_random_peps!(psi::InfinitePEPS) return InfinitePEPS(A′) end -@testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrgscheme)" for ( - T, unitcell, ctmrgscheme -) in Iterators.product( - scalartypes, unitcells, schemes +@testset "Trivial symmetry ($T) - ($unitcell) - ($flavor)" for (T, unitcell, flavor) in + Iterators.product( + scalartypes, unitcells, ctmrg_flavors ) physical_space = ComplexSpace(2) peps_space = ComplexSpace(2) @@ -50,7 +49,7 @@ end Random.seed!(987654321) # Seed RNG to make random environment consistent ctm = CTMRGEnv(psi, ctm_space) - alg = CTMRG(; maxiter, ctmrgscheme) + alg = CTMRG(; maxiter, flavor) ctm = leading_boundary(ctm, psi, alg) ctm2, = ctmrg_iter(psi, ctm, alg) @@ -58,9 +57,9 @@ end @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol end -@testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrgscheme)" for (T, unitcell, ctmrgscheme) in - Iterators.product( - scalartypes, unitcells, schemes +@testset "Z2 symmetry ($T) - ($unitcell) - ($flavor)" for (T, unitcell, flavor) in + Iterators.product( + scalartypes, unitcells, ctrmg_flavors ) physical_space = Z2Space(0 => 1, 1 => 1) peps_space = Z2Space(0 => 1, 1 => 1) @@ -74,7 +73,7 @@ end psi = InfinitePEPS(physical_space, peps_space; unitcell) ctm = CTMRGEnv(psi, ctm_space) - alg = CTMRG(; maxiter, ctmrgscheme) + alg = CTMRG(; maxiter, flavor) ctm = leading_boundary(ctm, psi, alg) ctm2, = ctmrg_iter(psi, ctm, alg) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 2c18ef90..46677810 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -18,8 +18,7 @@ names = ["Heisenberg", "p-wave superconductor"] gradtol = 1e-4 boundary_algs = [ - CTMRG(; verbosity=0, ctmrgscheme=:simultaneous), - CTMRG(; verbosity=0, ctmrgscheme=:sequential), + CTMRG(; verbosity=0, flavor=:simultaneous), CTMRG(; verbosity=0, flavor=:sequential) ] gradmodes = [ [ diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl index 846ea90a..4af4627a 100644 --- a/test/ctmrg/gradparts.jl +++ b/test/ctmrg/gradparts.jl @@ -30,9 +30,10 @@ Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbon Espaces = [ComplexSpace(χenv), Vect[FermionParity](0 => χenv / 2, 1 => χenv / 2)] tol = 1e-10 atol = 1e-6 +verbosity = 0 boundary_algs = [ - CTMRG(; tol, verbosity=0, ctmrgscheme=:simultaneous), - CTMRG(; tol, verbosity=0, ctmrgscheme=:sequential), + CTMRG(; tol, verbosity, flavor=:simultaneous), + CTMRG(; tol, verbosity, flavor=:sequential), ] ## Gauge invariant function of the environment diff --git a/test/runtests.jl b/test/runtests.jl index 07a4cd2f..49d2fc0c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,7 +27,7 @@ end include("ctmrg/unitcell.jl") end @time @safetestset "CTMRG schemes" begin - include("ctmrg/ctmrgschemes.jl") + include("ctmrg/flavors.jl") end end if GROUP == "ALL" || GROUP == "BOUNDARYMPS" From d9d0b05a606dfebe05f01a147a3df2f45a6d9756 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Fri, 6 Dec 2024 12:11:57 +0100 Subject: [PATCH 08/19] Fix typos, make runnable, replace pullback with rrule_via_ad --- src/PEPSKit.jl | 9 ++- .../contractions/ctmrg_contractions.jl | 56 ++++++++++--------- src/algorithms/ctmrg/ctmrg.jl | 19 ++++--- src/algorithms/ctmrg/sequential.jl | 2 +- src/algorithms/ctmrg/simultaneous.jl | 16 +++--- src/algorithms/ctmrg/sparse_environments.jl | 20 +++++-- src/algorithms/peps_opt.jl | 48 +++++++--------- src/environments/ctmrg_environments.jl | 2 +- test/ctmrg/fixed_iterscheme.jl | 14 ++--- test/ctmrg/gaugefix.jl | 6 +- test/ctmrg/gradparts.jl | 10 ++-- test/ctmrg/unitcell.jl | 4 +- 12 files changed, 112 insertions(+), 94 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index a1875fc1..0b2d11d4 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -41,6 +41,8 @@ include("algorithms/contractions/ctmrg_contractions.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") +include("algorithms/ctmrg/simultaneous.jl") +include("algorithms/ctmrg/sequential.jl") include("algorithms/ctmrg/gaugefix.jl") include("algorithms/toolbox.jl") @@ -94,7 +96,7 @@ Module containing default values that represent typical algorithm parameters. """ module Defaults using TensorKit, KrylovKit, OptimKit, OhMyThreads - using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint + using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint, HalfInfiniteProjector const ctmrg_maxiter = 100 const ctmrg_miniter = 4 const ctmrg_tol = 1e-8 @@ -107,6 +109,7 @@ module Defaults const fwd_alg = TensorKit.SDD() const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1) const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg) + const projector_alg = HalfInfiniteProjector const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -162,9 +165,9 @@ end using .Defaults: set_scheduler! export set_scheduler! export SVDAdjoint, IterSVD, NonTruncSVDAdjoint -export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length +export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector, CTMRG, CTMRGEnv export LocalOperator -export expectation_value, costfun, product_peps +export expectation_value, costfun, product_peps, correlation_length export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver export fixedpoint diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 0a84ae06..69921868 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -163,31 +163,6 @@ end # Projector contractions # ---------------------- -""" - left_and_right_projector(U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} - left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner) - -Compute left and right projectors based on a SVD and quadrant tensors, specified either as -`AbstractTensorMap`s or sparsely as `EnlargedCorner`s such that the quadrants are never -constructed explicitly. -""" -function left_and_right_projector( - U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} -) where {E<:ElementarySpace} - isqS = sdiag_inv_sqrt(S) - P_left = Q_next * V' * isqS - P_right = isqS * U' * Q - return P_left, P_right -end -function left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner) - isqS = sdiag_inv_sqrt(S) - P_left = left_projector(Q.E_1, Q.C, Q.E_2, V, isqS, Q.ket, Q.bra) - P_right = right_projector( - Q_next.E_1, Q_next.C, Q_next.E_2, U, isqS, Q_next.ket, Q_next.bra - ) - return P_left, P_right -end - """ left_projector(E_1, C, E_2, V, isqS, ket::PEPSTensor, bra::PEPSTensor=ket) @@ -234,6 +209,37 @@ function right_projector(E_1, C, E_2, U, isqS, ket::PEPSTensor, bra::PEPSTensor= E_1[χ4 D5 D6; χ_out] end +""" + left_and_right_projector(U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} + +Compute left and right projectors based on a SVD and quadrant tensors, where the inverse +square root `isqS` of the singular values is computed. + +Left projector: +``` + -- |~~~~~~~~| -- |~~| + | Q_next | |V'| -- isqS -- + == |~~~~~~~~| == |~~| +``` + +Right projector: +``` + |~~| -- |~~~| -- + -- isqS -- |U'| | Q | + |~~| == |~~~| == +``` +""" +function left_and_right_projector( + U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} +) where {E<:ElementarySpace} + isqS = sdiag_inv_sqrt(S) + @autoopt @tensor P_left[χ_in D_inabove D_inbelow; χ_out] := + Q_next[χ_in D_inabove D_inbelow; χ1 D1 D2] * conj(V[χ2; χ1 D1 D2]) * isqS[χ2; χ_out] + @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := + isqS[χ_in; χ1] * conj(U[χ2 D1 D2; χ1]) * Q[χ2 D1 D2; χ_out D_outabove D_outbelow] + return P_left, P_right +end + """ half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 445ee14e..c9aecc1d 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -102,6 +102,14 @@ function CTMRG(; ) end +function ctmrg_iteration(state, env, alg::CTMRG) + if alg.flavor == :simultaneous + return simultaneous_ctmrg_iter(state, env, alg) + elseif alg.flavor == :sequential + return sequential_ctmrg_iter(state, env, alg) + end +end + """ MPSKit.leading_boundary([envinit], state, alg::CTMRG) @@ -120,15 +128,10 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) env = deepcopy(envinit) log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) - f = if alg.flavor == :sequential - sequential_ctmrg_iter - elseif alg.flavor == :simultaneous - simultaneous_ctmrg_iter - end return LoggingExtras.withlevel(; alg.verbosity) do ctmrg_loginit!(log, η, N) for iter in 1:(alg.maxiter) - env, = f(state, env, alg) # Grow and renormalize in all 4 directions + env, = ctmrg_iteration(state, env, alg) # Grow and renormalize in all 4 directions η, CS, TS = calc_convergence(env, CS, TS) N = norm(state, env) @@ -166,8 +169,8 @@ and the given coordinate using the specified `alg`. function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) # SVD half-infinite environment halfinf = half_infinite_environment(enlarged_corners...) - trscheme = truncation_scheme(projector_alg, space(enlarged_corners[2], 1)) - svd_alg = svd_algorithm(projector_alg, coordinate) + trscheme = truncation_scheme(alg, space(enlarged_corners[2], 1)) + svd_alg = svd_algorithm(alg, coordinate) U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=trscheme) # Compute SVD truncation error and check for degenerate singular values diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 7631fefe..efebc594 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -27,7 +27,7 @@ function sequential_projectors( # SVD half-infinite environment column-wise coordinates = eachcoordinate(envs)[:, col] projectors = dtmap(coordinates) do (r, c) - proj, info = simultaneous_projector(state, envs, (WEST, r, c), alg) + proj, info = sequential_projector(state, envs, (WEST, r, c), alg) ϵ = max(ϵ, info.err / norm(info.S)) return proj end diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index 43f34fc9..f94bef38 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -22,10 +22,10 @@ function simultaneous_projectors( ϵ = zero(real(scalartype(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate - proj, info = compute_projector(enlarged_corners, coordinate, alg) - U[dir, r, c] = info.U - S[dir, r, c] = info.S - V[dir, r, c] = info.V + proj, info = simultaneous_projector(enlarged_corners, coordinate, alg) + U[coordinate...] = info.U + S[coordinate...] = info.S + V[coordinate...] = info.V ϵ = max(ϵ, info.err / norm(info.S)) return proj end @@ -37,7 +37,7 @@ end function simultaneous_projector( enlarged_corners::Array{E,3}, coordinate, alg::HalfInfiniteProjector ) where {E} - coordinate′ = _next_projector_coordinate(coordinate, size(enlarged_corners)[2:3]...) + coordinate′ = _next_coordinate(coordinate, size(enlarged_corners)[2:3]...) ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) return compute_projector(ec, coordinate, alg) end @@ -45,9 +45,9 @@ function simultaneous_projector( enlarged_corners::Array{E,3}, coordinate, alg::FullInfiniteProjector ) where {E} rowsize, colsize = size(enlarged_corners)[2:3] - coordinate2 = _next_projector_coordinate(coordinate, rowsize, colsize) - coordinate3 = _next_projector_coordinate(coordinate2, rowsize, colsize) - coordinate4 = _next_projector_coordinate(coordinate3, rowsize, colsize) + coordinate2 = _next_coordinate(coordinate, rowsize, colsize) + coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) + coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) ec = ( enlarged_corners[coordinate...], enlarged_corners[coordinate2...], diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 4f9fbad2..7b9edab9 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -97,6 +97,21 @@ function renormalize_southwest_corner(ec::EnlargedCorner, P_left, P_right) ) end +# Wrapper around half_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) +function half_infinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) + return HalfInfiniteEnv(ec_1, ec_2) +end + +# Compute left and right projectors sparsely without constructing enlarged corners explicitly +function left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner) + isqS = sdiag_inv_sqrt(S) + P_left = left_projector(Q.E_1, Q.C, Q.E_2, V, isqS, Q.ket, Q.bra) + P_right = right_projector( + Q_next.E_1, Q_next.C, Q_next.E_2, U, isqS, Q_next.ket, Q_next.bra + ) + return P_left, P_right +end + # -------------------------------- # Sparse half-infinite environment # -------------------------------- @@ -193,11 +208,6 @@ function (env::HalfInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * ) end -# Wrapper around half_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) -function half_infinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) - return HalfInfiniteEnv(ec_1, ec_2) -end - # AbstractTensorMap subtyping and IterSVD compatibility function TensorKit.domain(env::HalfInfiniteEnv) return domain(env.E_4) * domain(env.ket_2)[3] * domain(env.bra_2)[3]' diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 00eac45f..91e5adc8 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -95,13 +95,13 @@ struct PEPSOptimize{G} gradient_alg::G function PEPSOptimize( # Inner constructor to prohibit illegal setting combinations - boundary_alg::CTMRG{S}, + boundary_alg::CTMRG, optimizer, reuse_env, gradient_alg::G, - ) where {S,G} + ) where {G} if gradient_alg isa GradMode - if S === :sequential && iterscheme(gradient_alg) === :fixed + if boundary_alg.flavor === :sequential && iterscheme(gradient_alg) === :fixed throw(ArgumentError(":sequential and :fixed are not compatible")) end end @@ -216,26 +216,24 @@ Evaluating the gradient of the cost function for CTMRG: function _rrule( gradmode::GradMode{:diffgauge}, - ::RuleConfig, + config::RuleConfig, ::typeof(MPSKit.leading_boundary), envinit, state, alg::CTMRG, ) - envs = leading_boundary(envinit, state, alg) #TODO: fixed space for unit cells + envs = leading_boundary(envinit, state, alg) function leading_boundary_diffgauge_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) # find partial gradients of gauge_fixed single CTMRG iteration - # TODO: make this rrule_via_ad so it's zygote-agnostic - _, env_vjp = pullback(state, envs) do A, x - return gauge_fix(x, ctmrg_iter(A, x, alg)[1])[1] - end + f(A, x) = gauge_fix(x, ctmrg_iteration(A, x, alg)[1])[1] + _, env_vjp = rrule_via_ad(config, f, state, envs) # evaluate the geometric sum - ∂f∂A(x)::typeof(state) = env_vjp(x)[1] - ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂f∂A(x)::typeof(state) = env_vjp(x)[2] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[3] ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() @@ -247,44 +245,40 @@ end # Here f is differentiated from an pre-computed SVD with fixed U, S and V function _rrule( gradmode::GradMode{:fixed}, - ::RuleConfig, + config::RuleConfig, ::typeof(MPSKit.leading_boundary), envinit, state, - alg::CTMRG{C}, -) where {C} - @assert C === :simultaneous + alg::CTMRG, +) + @assert alg.flavor === :simultaneous @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg) envs = leading_boundary(envinit, state, alg) - envsconv, info = ctmrg_iter(state, envs, alg) - envsfix, signs = gauge_fix(envs, envsconv) + envsconv, info = ctmrg_iteration(state, envs, alg) + envs_fixed, signs = gauge_fix(envs, envsconv) # Fix SVD Ufix, Vfix = fix_relative_phases(info.U, info.V, signs) svd_alg_fixed = SVDAdjoint(; fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) - alg_fixed = CTMRG(; - svd_alg=svd_alg_fixed, trscheme=notrunc(), flavor=:simultaneous - ) + alg_fixed = CTMRG(; svd_alg=svd_alg_fixed, trscheme=notrunc(), flavor=:simultaneous) function leading_boundary_fixed_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) - _, env_vjp = pullback(state, envsfix) do A, x - e, = ctmrg_iter(A, x, alg_fixed) - return fix_global_phases(x, e) - end + f(A, x) = fix_global_phases(x, ctmrg_iteration(A, x, alg_fixed)[1]) + _, env_vjp = rrule_via_ad(config, f, state, envs_fixed) # evaluate the geometric sum - ∂f∂A(x)::typeof(state) = env_vjp(x)[1] - ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂f∂A(x)::typeof(state) = env_vjp(x)[2] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[3] ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() end - return envsfix, leading_boundary_fixed_pullback + return envs_fixed, leading_boundary_fixed_pullback end @doc """ diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 29a823cc..6819159e 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -354,7 +354,7 @@ function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnv, name::Symbol) throw(ArgumentError("No rrule for getproperty of $name")) end end -Base.size(env::CTMRG, args...) = size(env.corners, args...) +Base.size(env::CTMRGEnv, args...) = size(env.corners, args...) Base.eltype(env::CTMRGEnv) = eltype(env.corners[1]) Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) function eachcoordinate(x::CTMRGEnv) diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl index 3f005986..b5f9a291 100644 --- a/test/ctmrg/fixed_iterscheme.jl +++ b/test/ctmrg/fixed_iterscheme.jl @@ -5,7 +5,7 @@ using TensorKit, KrylovKit using PEPSKit using PEPSKit: FixedSVD, - ctmrg_iter, + ctmrg_iteration, gauge_fix, fix_relative_phases, fix_global_phases, @@ -30,7 +30,7 @@ unitcells = [(1, 1), (3, 4)] env_conv1 = leading_boundary(CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg) # do extra iteration to get SVD - env_conv2, info = ctmrg_iter(psi, env_conv1, ctm_alg) + env_conv2, info = ctmrg_iteration(psi, env_conv1, ctm_alg) env_fix, signs = gauge_fix(env_conv1, env_conv2) @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = 1e-6 @@ -40,7 +40,7 @@ unitcells = [(1, 1), (3, 4)] ctm_alg_fix = CTMRG(; svd_alg=svd_alg_fix, trscheme=notrunc()) # do iteration with FixedSVD - env_fixedsvd, = ctmrg_iter(psi, env_conv1, ctm_alg_fix) + env_fixedsvd, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix) env_fixedsvd = fix_global_phases(env_conv1, env_fixedsvd) @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = 1e-6 end @@ -59,11 +59,11 @@ end env_conv1 = leading_boundary(env_init, psi, ctm_alg_iter) # do extra iteration to get SVD - env_conv2_iter, info_iter = ctmrg_iter(psi, env_conv1, ctm_alg_iter) + env_conv2_iter, info_iter = ctmrg_iteration(psi, env_conv1, ctm_alg_iter) env_fix_iter, signs_iter = gauge_fix(env_conv1, env_conv2_iter) @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = 1e-6 - env_conv2_full, info_full = ctmrg_iter(psi, env_conv1, ctm_alg_full) + env_conv2_full, info_full = ctmrg_iteration(psi, env_conv1, ctm_alg_full) env_fix_full, signs_full = gauge_fix(env_conv1, env_conv2_full) @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = 1e-6 @@ -77,11 +77,11 @@ end ctm_alg_fix_full = CTMRG(; svd_alg=svd_alg_fix_full, trscheme=notrunc()) # do iteration with FixedSVD - env_fixedsvd_iter, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_iter) + env_fixedsvd_iter, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_iter) env_fixedsvd_iter = fix_global_phases(env_conv1, env_fixedsvd_iter) @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = 1e-6 # This doesn't work for x₀ = rand(size(b, 1))? - env_fixedsvd_full, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_full) + env_fixedsvd_full, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_full) env_fixedsvd_full = fix_global_phases(env_conv1, env_fixedsvd_full) @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = 1e-6 diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 748c5b05..8757b231 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -3,7 +3,7 @@ using Random using PEPSKit using TensorKit -using PEPSKit: ctmrg_iter, gauge_fix, calc_elementwise_convergence +using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] @@ -52,7 +52,7 @@ end alg = CTMRG(; maxiter, flavor) ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iter(psi, ctm, alg) + ctm2, = ctmrg_iteration(psi, ctm, alg) ctm_fixed, = gauge_fix(ctm, ctm2) @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol end @@ -76,7 +76,7 @@ end alg = CTMRG(; maxiter, flavor) ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iter(psi, ctm, alg) + ctm2, = ctmrg_iteration(psi, ctm, alg) ctm_fixed, = gauge_fix(ctm, ctm2) @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol end diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl index 4af4627a..795d52e0 100644 --- a/test/ctmrg/gradparts.jl +++ b/test/ctmrg/gradparts.jl @@ -13,7 +13,7 @@ using PEPSKit: SOUTHWEST, rotate_north, left_move, - ctmrg_iter + ctmrg_iteration using Zygote using ChainRulesCore @@ -53,15 +53,17 @@ end ## Tests # ------ -@testset "Reverse rules for ctmrg_iter with spacetype $(Vspaces[i])" for i in - eachindex(Pspaces) +@testset "Reverse rules for ctmrg_iteration with spacetype $(Vspaces[i])" for i in + eachindex( + Pspaces +) Random.seed!(42039482030) psi = InfinitePEPS(Pspaces[i], Vspaces[i], Vspaces[i]) env = CTMRGEnv(psi, Espaces[i]) @testset "$alg" for alg in boundary_algs @info "$(typeof(alg)) on $(Vspaces[i])" - f(state, env) = rho(ctmrg_iter(state, env, alg)[1]) + f(state, env) = rho(ctmrg_iteration(state, env, alg)[1]) # use rrule testing functionality but compute rrule via Zygote test_rrule( diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 0ffbdc09..1ba2ea40 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -1,7 +1,7 @@ using Test using Random using PEPSKit -using PEPSKit: _prev, _next, ctmrg_iter +using PEPSKit: _prev, _next, ctmrg_iteration using TensorKit # settings @@ -16,7 +16,7 @@ function test_unitcell( env = CTMRGEnv(randn, stype, peps, chis_north, chis_east, chis_south, chis_west) # apply one CTMRG iteration with fixeds - env′, = ctmrg_iter(peps, env, ctm_alg) + env′, = ctmrg_iteration(peps, env, ctm_alg) # compute random expecation value to test matching bonds random_op = LocalOperator( From 85c2cf88859ae2ee84dbef2d923d0fbe706befda Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Fri, 6 Dec 2024 18:17:33 +0100 Subject: [PATCH 09/19] Make FullInfiniteProjector runnable --- .../contractions/ctmrg_contractions.jl | 18 ++++++++++--- src/algorithms/ctmrg/ctmrg.jl | 27 ++++++++++--------- src/algorithms/ctmrg/simultaneous.jl | 25 ++++++++++++----- src/utility/util.jl | 22 --------------- 4 files changed, 47 insertions(+), 45 deletions(-) diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 69921868..8e784dab 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -217,9 +217,9 @@ square root `isqS` of the singular values is computed. Left projector: ``` - -- |~~~~~~~~| -- |~~| - | Q_next | |V'| -- isqS -- - == |~~~~~~~~| == |~~| + -- |~~~~~~| -- |~~| + |Q_next| |V'| -- isqS -- + == |~~~~~~| == |~~| ``` Right projector: @@ -235,10 +235,20 @@ function left_and_right_projector( isqS = sdiag_inv_sqrt(S) @autoopt @tensor P_left[χ_in D_inabove D_inbelow; χ_out] := Q_next[χ_in D_inabove D_inbelow; χ1 D1 D2] * conj(V[χ2; χ1 D1 D2]) * isqS[χ2; χ_out] - @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := + @autoopt @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := isqS[χ_in; χ1] * conj(U[χ2 D1 D2; χ1]) * Q[χ2 D1 D2; χ_out D_outabove D_outbelow] return P_left, P_right end +function left_and_right_projector( + U, S, V, R_left::AbstractTensorMap{E,1,3}, L_right::AbstractTensorMap{E,3,1} +) where {E<:ElementarySpace} + isqS = sdiag_inv_sqrt(S) + @tensor P_left[χ_in D_inabove D_inbelow; χ_out] := + L_right[χ_in D_inabove D_inbelow; full] * conj(V[χ2; full]) * isqS[χ2; χ_out] + @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := + isqS[χ_in; χ1] * conj(U[full; χ1]) * R_left[full; χ_out D_outabove D_outbelow] + return P_left, P_right +end """ half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index c9aecc1d..ad4648ec 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -89,6 +89,7 @@ function CTMRG(; miniter=Defaults.ctmrg_miniter, flavor=Defaults.ctmrg_flavor, verbosity=2, + projector_alg=Defaults.projector_alg, svd_alg=Defaults.svd_alg, trscheme=Defaults.trscheme, ) @@ -98,7 +99,7 @@ function CTMRG(; miniter, flavor, verbosity, - Defaults.projector_alg(; svd_alg, trscheme, verbosity), + projector_alg(; svd_alg, trscheme, verbosity), ) end @@ -185,17 +186,17 @@ function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjec return (P_left, P_right), (; err, U, S, V) end function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) - # QR top and bottom half-infinite environments - halfinf_top = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) - halfinf_bot = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) - _, R_top = leftorth!(halfinf_top) - _, R_bot = leftorth!(halfinf_bot) - R_fullinf = R_top * R_bot - - # SVD product of R's - trscheme = truncation_scheme(projector_alg, space(enlarged_corners[4], 1)) - svd_alg = svd_algorithm(projector_alg, coordinate) - U, S, V, err = PEPSKit.tsvd!(R_fullinf, svd_alg; trunc=trscheme) + # QR left and right half-infinite environments (cut placement is consistent with half-infinite proj!) + halfinf_left = half_infinite_environment(enlarged_corners[4], enlarged_corners[1]) + halfinf_right = half_infinite_environment(enlarged_corners[2], enlarged_corners[3]) + _, R_left = leftorth!(halfinf_left) + L_right, _ = rightorth!(halfinf_right) + + # SVD product of QRs + fullinf = R_left * L_right + trscheme = truncation_scheme(alg, space(enlarged_corners[4], 1)) + svd_alg = svd_algorithm(alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=trscheme) # Compute SVD truncation error and check for degenerate singular values Zygote.isderiving() && ignore_derivatives() do @@ -205,6 +206,6 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec end end - P_left, P_right = left_and_right_projector(U, S, V, R_top, R_bot) + P_left, P_right = left_and_right_projector(U, S, V, R_left, L_right) return (P_left, P_right), (; err, U, S, V) end diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index f94bef38..e69bc51c 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -13,12 +13,25 @@ function simultaneous_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) return envs′, info end -function simultaneous_projectors( - enlarged_corners, envs::CTMRGEnv{C,E}, alg::ProjectorAlgs -) where {C,E} - # pre-allocation - U, V = Zygote.Buffer.(projector_type(envs.edges)) - S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) # Corner type but with real numbers +# Pre-allocate U, S, and V tensor as Zygote buffers to make it differentiable +function _prealloc_svd(edges::Array{E,N}, ::HalfInfiniteProjector) where {E,N} + Sc = scalartype(E) + U = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, space(e)), edges)) + V = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, domain(e), codomain(e)), edges)) + S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc))) # Corner type but with real numbers + return U, S, V +end +function _prealloc_svd(edges::Array{E,N}, ::FullInfiniteProjector) where {E,N} + Sc = scalartype(E) + Rspace(x) = spacetype(E)(dim(codomain(x))) + U = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, Rspace(e), domain(e)), edges)) + V = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, domain(e), Rspace(e)), edges)) + S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc))) # Corner type but with real numbers + return U, S, V +end + +function simultaneous_projectors(enlarged_corners, envs::CTMRGEnv, alg::ProjectorAlgs) + U, S, V = _prealloc_svd(envs.edges, alg) ϵ = zero(real(scalartype(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate diff --git a/src/utility/util.jl b/src/utility/util.jl index 93cb9577..858f8e86 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -92,28 +92,6 @@ function is_degenerate_spectrum( return false end -""" - projector_type(T::DataType, size) - projector_type(edges::Array{<:AbstractTensorMap}) - -Create two arrays of specified `size` that contain undefined tensors representing -left and right acting projectors, respectively. The projector types are inferred -from the TensorMap type `T` which avoids having to recompute transpose tensors. -Alternatively, supply an array of edge tensors from which left and right projectors -are intialized explicitly with zeros. -""" -function projector_type(T::DataType, size) - P_left = Array{T,length(size)}(undef, size) - Prtype = tensormaptype(spacetype(T), numin(T), numout(T), storagetype(T)) - P_right = Array{Prtype,length(size)}(undef, size) - return P_left, P_right -end -function projector_type(edges::Array{<:AbstractTensorMap}) - P_left = map(e -> TensorMap(zeros, scalartype(e), space(e)), edges) - P_right = map(e -> TensorMap(zeros, scalartype(e), domain(e), codomain(e)), edges) - return P_left, P_right -end - # There are no rrules for rotl90 and rotr90 in ChainRules.jl function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix) function rotl90_pullback(x) From 1d75b48adb889accb9531ff7de00415e0c0303d2 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Fri, 6 Dec 2024 18:35:54 +0100 Subject: [PATCH 10/19] Add docstrings --- src/algorithms/ctmrg/sequential.jl | 26 +++++++++++++++++++------- src/algorithms/ctmrg/simultaneous.jl | 22 ++++++++++++++++++---- 2 files changed, 37 insertions(+), 11 deletions(-) diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index efebc594..0632b7ea 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -19,23 +19,30 @@ function sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) return envs, (; err=ϵ) end +""" + sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) + sequential_projectors(coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) + +Compute CTMRG projectors in the `:sequential` scheme either for an entire column `col` or +for a specific `coordinate` (where `dir=WEST` is already implied in the `:sequential` scheme). +""" function sequential_projectors( - col::Int, state::InfinitePEPS, envs::CTMRGEnv{C,E}, alg::ProjectorAlgs -) where {C,E} + col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs +) ϵ = zero(real(scalartype(envs))) # SVD half-infinite environment column-wise coordinates = eachcoordinate(envs)[:, col] projectors = dtmap(coordinates) do (r, c) - proj, info = sequential_projector(state, envs, (WEST, r, c), alg) + proj, info = sequential_projectors((WEST, r, c), state, envs, alg) ϵ = max(ϵ, info.err / norm(info.S)) return proj end return (map(first, projectors), map(last, projectors)), (; err=ϵ) end -function sequential_projector( - state::InfinitePEPS, envs::CTMRGEnv, coordinate, alg::HalfInfiniteProjector +function sequential_projectors( + coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::HalfInfiniteProjector ) _, r, c = coordinate r′ = _prev(r, size(envs, 2)) @@ -43,8 +50,8 @@ function sequential_projector( Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST) return compute_projector((Q1, Q2), coordinate, alg) end -function sequential_projector( - state::InfinitePEPS, envs::CTMRGEnv, coordinate, alg::FullInfiniteProjector +function sequential_projectors( + coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::FullInfiniteProjector ) _, r, c = coordinate r′ = _next(r, size(envs, 2)) @@ -56,6 +63,11 @@ function sequential_projector( return compute_projector((Q1, Q2, Q3, Q4), coordinate, alg) end +""" + renormalize_sequentially(col::Int, projectors, state, envs) + +Renormalize one column of the CTMRG environment. +""" function renormalize_sequentially(col::Int, projectors, state, envs) corners = Zygote.Buffer(envs.corners) edges = Zygote.Buffer(envs.edges) diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index e69bc51c..19828e4f 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -30,12 +30,21 @@ function _prealloc_svd(edges::Array{E,N}, ::FullInfiniteProjector) where {E,N} return U, S, V end -function simultaneous_projectors(enlarged_corners, envs::CTMRGEnv, alg::ProjectorAlgs) +""" + simultaneous_projectors(enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgs) + simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgs) + +Compute CTMRG projectors in the `:simultaneous` scheme either for all provided +enlarged corners or on a specific `coordinate`. +""" +function simultaneous_projectors( + enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgs +) where {E} U, S, V = _prealloc_svd(envs.edges, alg) ϵ = zero(real(scalartype(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate - proj, info = simultaneous_projector(enlarged_corners, coordinate, alg) + proj, info = simultaneous_projectors(enlarged_corners, coordinate, alg) U[coordinate...] = info.U S[coordinate...] = info.S V[coordinate...] = info.V @@ -47,14 +56,14 @@ function simultaneous_projectors(enlarged_corners, envs::CTMRGEnv, alg::Projecto P_right = map(last, projectors) return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) end -function simultaneous_projector( +function simultaneous_projectors( enlarged_corners::Array{E,3}, coordinate, alg::HalfInfiniteProjector ) where {E} coordinate′ = _next_coordinate(coordinate, size(enlarged_corners)[2:3]...) ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) return compute_projector(ec, coordinate, alg) end -function simultaneous_projector( +function simultaneous_projectors( enlarged_corners::Array{E,3}, coordinate, alg::FullInfiniteProjector ) where {E} rowsize, colsize = size(enlarged_corners)[2:3] @@ -70,6 +79,11 @@ function simultaneous_projector( return compute_projector(ec, coordinate, alg) end +""" + renormalize_simultaneously(enlarged_corners, projectors, state, envs) + +Renormalize all enlarged corners and edges simultaneously. +""" function renormalize_simultaneously(enlarged_corners, projectors, state, envs) P_left, P_right = projectors coordinates = eachcoordinate(envs, 1:4) From 16a9ee1a829b0d98ae8bc9e09e4f9e2e57b9d28b Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Fri, 6 Dec 2024 20:36:30 +0100 Subject: [PATCH 11/19] Fix FullInfiniteProjector for :sequential, add tests, remove old gradparts test --- src/algorithms/ctmrg/ctmrg.jl | 11 ++-- src/algorithms/ctmrg/sequential.jl | 34 ++++++++---- src/algorithms/ctmrg/simultaneous.jl | 10 ++-- test/ctmrg/fixed_iterscheme.jl | 26 ++++----- test/ctmrg/flavors.jl | 6 ++- test/ctmrg/gaugefix.jl | 43 +++++++-------- test/ctmrg/gradients.jl | 5 +- test/ctmrg/gradparts.jl | 79 ---------------------------- 8 files changed, 81 insertions(+), 133 deletions(-) delete mode 100644 test/ctmrg/gradparts.jl diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index ad4648ec..2485be29 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -103,6 +103,11 @@ function CTMRG(; ) end +""" + ctmrg_iteration(state, env, alg::CTMRG) + +Perform a single CTMRG iteration in which all directions are being grown and renormalized. +""" function ctmrg_iteration(state, env, alg::CTMRG) if alg.flavor == :simultaneous return simultaneous_ctmrg_iter(state, env, alg) @@ -186,9 +191,9 @@ function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjec return (P_left, P_right), (; err, U, S, V) end function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) - # QR left and right half-infinite environments (cut placement is consistent with half-infinite proj!) - halfinf_left = half_infinite_environment(enlarged_corners[4], enlarged_corners[1]) - halfinf_right = half_infinite_environment(enlarged_corners[2], enlarged_corners[3]) + # QR left and right half-infinite environments + halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) + halfinf_right = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) _, R_left = leftorth!(halfinf_left) L_right, _ = rightorth!(halfinf_right) diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 0632b7ea..60c27a4d 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -21,7 +21,7 @@ end """ sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) - sequential_projectors(coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) + sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) Compute CTMRG projectors in the `:sequential` scheme either for an entire column `col` or for a specific `coordinate` (where `dir=WEST` is already implied in the `:sequential` scheme). @@ -42,7 +42,10 @@ function sequential_projectors( return (map(first, projectors), map(last, projectors)), (; err=ϵ) end function sequential_projectors( - coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::HalfInfiniteProjector + coordinate::NTuple{3,Int}, + state::InfinitePEPS, + envs::CTMRGEnv, + alg::HalfInfiniteProjector, ) _, r, c = coordinate r′ = _prev(r, size(envs, 2)) @@ -51,16 +54,25 @@ function sequential_projectors( return compute_projector((Q1, Q2), coordinate, alg) end function sequential_projectors( - coordinate, state::InfinitePEPS, envs::CTMRGEnv, alg::FullInfiniteProjector + coordinate::NTuple{3,Int}, + state::InfinitePEPS, + envs::CTMRGEnv, + alg::FullInfiniteProjector, ) - _, r, c = coordinate - r′ = _next(r, size(envs, 2)) - c′ = _next(c, size(envs, 3)) - Q1 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r, c)), NORTHWEST) - Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHEAST, r, c′)), NORTHEAST) - Q3 = TensorMap(EnlargedCorner(state, envs, (SOUTHEAST, r′, c′)), SOUTHEAST) - Q4 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r′, c)), SOUTHWEST) - return compute_projector((Q1, Q2, Q3, Q4), coordinate, alg) + # _, r, c = coordinate + # r′ = _next(r, size(envs, 2)) + # c′ = _next(c, size(envs, 3)) + rowsize, colsize = size(envs)[2:3] + coordinate_nw = _next_coordinate(coordinate, rowsize, colsize) + coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize) + coordinate_se = _next_coordinate(coordinate_ne, rowsize, colsize) + ec = ( + TensorMap(EnlargedCorner(state, envs, coordinate_se), SOUTHEAST), + TensorMap(EnlargedCorner(state, envs, coordinate), SOUTHWEST), + TensorMap(EnlargedCorner(state, envs, coordinate_nw), NORTHWEST), + TensorMap(EnlargedCorner(state, envs, coordinate_ne), NORTHEAST), + ) + return compute_projector(ec, coordinate, alg) end """ diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index 19828e4f..b6992fbb 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -23,7 +23,7 @@ function _prealloc_svd(edges::Array{E,N}, ::HalfInfiniteProjector) where {E,N} end function _prealloc_svd(edges::Array{E,N}, ::FullInfiniteProjector) where {E,N} Sc = scalartype(E) - Rspace(x) = spacetype(E)(dim(codomain(x))) + Rspace(x) = fuse(codomain(x)) U = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, Rspace(e), domain(e)), edges)) V = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, domain(e), Rspace(e)), edges)) S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc))) # Corner type but with real numbers @@ -44,7 +44,7 @@ function simultaneous_projectors( ϵ = zero(real(scalartype(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate - proj, info = simultaneous_projectors(enlarged_corners, coordinate, alg) + proj, info = simultaneous_projectors(coordinate, enlarged_corners, alg) U[coordinate...] = info.U S[coordinate...] = info.S V[coordinate...] = info.V @@ -57,24 +57,24 @@ function simultaneous_projectors( return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) end function simultaneous_projectors( - enlarged_corners::Array{E,3}, coordinate, alg::HalfInfiniteProjector + coordinate, enlarged_corners::Array{E,3}, alg::HalfInfiniteProjector ) where {E} coordinate′ = _next_coordinate(coordinate, size(enlarged_corners)[2:3]...) ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) return compute_projector(ec, coordinate, alg) end function simultaneous_projectors( - enlarged_corners::Array{E,3}, coordinate, alg::FullInfiniteProjector + coordinate, enlarged_corners::Array{E,3}, alg::FullInfiniteProjector ) where {E} rowsize, colsize = size(enlarged_corners)[2:3] coordinate2 = _next_coordinate(coordinate, rowsize, colsize) coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) ec = ( + enlarged_corners[coordinate4...], enlarged_corners[coordinate...], enlarged_corners[coordinate2...], enlarged_corners[coordinate3...], - enlarged_corners[coordinate4...], ) return compute_projector(ec, coordinate, alg) end diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl index b5f9a291..f13b04a6 100644 --- a/test/ctmrg/fixed_iterscheme.jl +++ b/test/ctmrg/fixed_iterscheme.jl @@ -15,14 +15,17 @@ using PEPSKit: χbond = 2 χenv = 16 svd_algs = [SVDAdjoint(; fwd_alg=TensorKit.SDD()), SVDAdjoint(; fwd_alg=IterSVD())] +projector_algs = [HalfInfiniteProjector] #, FullInfiniteProjector] unitcells = [(1, 1), (3, 4)] +atol = 1e-5 # test for element-wise convergence after application of fixed step -@testset "$unitcell unit cell with $(typeof(svd_alg.fwd_alg))" for (unitcell, svd_alg) in - Iterators.product( - unitcells, svd_algs +@testset "$unitcell unit cell with $(typeof(svd_alg.fwd_alg)) and $projector_alg" for ( + unitcell, svd_alg, projector_alg +) in Iterators.product( + unitcells, svd_algs, projector_algs ) - ctm_alg = CTMRG(; svd_alg) + ctm_alg = CTMRG(; svd_alg, projector_alg) # initialize states Random.seed!(2394823842) @@ -32,17 +35,17 @@ unitcells = [(1, 1), (3, 4)] # do extra iteration to get SVD env_conv2, info = ctmrg_iteration(psi, env_conv1, ctm_alg) env_fix, signs = gauge_fix(env_conv1, env_conv2) - @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = atol # fix gauge of SVD U_fix, V_fix = fix_relative_phases(info.U, info.V, signs) svd_alg_fix = SVDAdjoint(; fwd_alg=FixedSVD(U_fix, info.S, V_fix)) - ctm_alg_fix = CTMRG(; svd_alg=svd_alg_fix, trscheme=notrunc()) + ctm_alg_fix = CTMRG(; projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc()) # do iteration with FixedSVD env_fixedsvd, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix) env_fixedsvd = fix_global_phases(env_conv1, env_fixedsvd) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = atol end @testset "Element-wise consistency of TensorKit.SDD and IterSVD" begin @@ -61,11 +64,11 @@ end # do extra iteration to get SVD env_conv2_iter, info_iter = ctmrg_iteration(psi, env_conv1, ctm_alg_iter) env_fix_iter, signs_iter = gauge_fix(env_conv1, env_conv2_iter) - @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = atol env_conv2_full, info_full = ctmrg_iteration(psi, env_conv1, ctm_alg_full) env_fix_full, signs_full = gauge_fix(env_conv1, env_conv2_full) - @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = atol # fix gauge of SVD U_fix_iter, V_fix_iter = fix_relative_phases(info_iter.U, info_iter.V, signs_iter) @@ -79,14 +82,13 @@ end # do iteration with FixedSVD env_fixedsvd_iter, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_iter) env_fixedsvd_iter = fix_global_phases(env_conv1, env_fixedsvd_iter) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = 1e-6 # This doesn't work for x₀ = rand(size(b, 1))? + @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = atol # This doesn't work for x₀ = rand(size(b, 1))? env_fixedsvd_full, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_full) env_fixedsvd_full = fix_global_phases(env_conv1, env_fixedsvd_full) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = atol # check matching decompositions - atol = 1e-12 decomposition_check = all( zip(info_iter.U, info_iter.S, info_iter.V, info_full.U, info_full.S, info_full.V), ) do (U_iter, S_iter, V_iter, U_full, S_full, V_full) diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl index c1c848a9..9337bc08 100644 --- a/test/ctmrg/flavors.jl +++ b/test/ctmrg/flavors.jl @@ -10,8 +10,12 @@ using PEPSKit ctm_alg_sequential = CTMRG(; flavor=:sequential) ctm_alg_simultaneous = CTMRG(; flavor=:simultaneous) unitcells = [(1, 1), (3, 4)] +projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] -@testset "$(unitcell) unit cell" for unitcell in unitcells +@testset "$(unitcell) unit cell with $projector_alg" for (unitcell, projector_alg) in + Iterators.product( + unitcells, projector_algs +) # compute environments Random.seed!(32350283290358) psi = InfinitePEPS(2, χbond; unitcell) diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 8757b231..040d72ed 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -7,8 +7,9 @@ using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] -maxiter = 200 +maxiter = 400 ctmrg_flavors = [:simultaneous, :sequential] +projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] χ = 6 atol = 1e-4 @@ -34,9 +35,10 @@ function _semi_random_peps!(psi::InfinitePEPS) return InfinitePEPS(A′) end -@testset "Trivial symmetry ($T) - ($unitcell) - ($flavor)" for (T, unitcell, flavor) in - Iterators.product( - scalartypes, unitcells, ctmrg_flavors +@testset "Trivial symmetry ($T) - ($unitcell) - ($flavor) - ($projector_alg)" for ( + T, unitcell, flavor, projector_alg +) in Iterators.product( + scalartypes, unitcells, ctmrg_flavors, projector_algs ) physical_space = ComplexSpace(2) peps_space = ComplexSpace(2) @@ -47,19 +49,19 @@ end _make_symmetric!(psi) Random.seed!(987654321) # Seed RNG to make random environment consistent - ctm = CTMRGEnv(psi, ctm_space) + env = CTMRGEnv(psi, ctm_space) + alg = CTMRG(; maxiter, flavor, projector_alg) - alg = CTMRG(; maxiter, flavor) - - ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iteration(psi, ctm, alg) - ctm_fixed, = gauge_fix(ctm, ctm2) - @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol + env = leading_boundary(env, psi, alg) + env′, = ctmrg_iteration(psi, env, alg) + env_fixed, = gauge_fix(env, env′) + @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol end -@testset "Z2 symmetry ($T) - ($unitcell) - ($flavor)" for (T, unitcell, flavor) in - Iterators.product( - scalartypes, unitcells, ctrmg_flavors +@testset "Z2 symmetry ($T) - ($unitcell) - ($flavor) - ($projector_alg)" for ( + T, unitcell, flavor, projector_alg +) in Iterators.product( + scalartypes, unitcells, ctmrg_flavors, projector_algs ) physical_space = Z2Space(0 => 1, 1 => 1) peps_space = Z2Space(0 => 1, 1 => 1) @@ -71,12 +73,11 @@ end Random.seed!(987654321) # Seed RNG to make random environment consistent psi = InfinitePEPS(physical_space, peps_space; unitcell) - ctm = CTMRGEnv(psi, ctm_space) - - alg = CTMRG(; maxiter, flavor) + env = CTMRGEnv(psi, ctm_space) + alg = CTMRG(; maxiter, flavor, projector_alg) - ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iteration(psi, ctm, alg) - ctm_fixed, = gauge_fix(ctm, ctm2) - @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol + env = leading_boundary(env, psi, alg) + env′, = ctmrg_iteration(psi, env, alg) + env_fixed, = gauge_fix(env, env′) + @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol end diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 46677810..a4ceafea 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -18,7 +18,10 @@ names = ["Heisenberg", "p-wave superconductor"] gradtol = 1e-4 boundary_algs = [ - CTMRG(; verbosity=0, flavor=:simultaneous), CTMRG(; verbosity=0, flavor=:sequential) + # CTMRG(; verbosity=0, flavor=:simultaneous, projector_alg=HalfInfiniteProjector), + CTMRG(; verbosity=0, flavor=:simultaneous, projector_alg=FullInfiniteProjector), + # CTMRG(; verbosity=0, flavor=:sequential, projector_alg=HalfInfiniteProjector), + CTMRG(; verbosity=0, flavor=:sequential, projector_alg=FullInfiniteProjector), ] gradmodes = [ [ diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl deleted file mode 100644 index 795d52e0..00000000 --- a/test/ctmrg/gradparts.jl +++ /dev/null @@ -1,79 +0,0 @@ -using Test -using Random -using PEPSKit -using TensorKit -using PEPSKit: - NORTH, - SOUTH, - WEST, - EAST, - NORTHWEST, - NORTHEAST, - SOUTHEAST, - SOUTHWEST, - rotate_north, - left_move, - ctmrg_iteration - -using Zygote -using ChainRulesCore -using ChainRulesTestUtils - -include(joinpath("..", "utility.jl")) - -## Test spaces, tested functions and CTMRG algorithm -# -------------------------------------------------- -χbond = 2 -χenv = 4 -Pspaces = [ComplexSpace(2), Vect[FermionParity](0 => 1, 1 => 1)] -Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbond / 2)] -Espaces = [ComplexSpace(χenv), Vect[FermionParity](0 => χenv / 2, 1 => χenv / 2)] -tol = 1e-10 -atol = 1e-6 -verbosity = 0 -boundary_algs = [ - CTMRG(; tol, verbosity, flavor=:simultaneous), - CTMRG(; tol, verbosity, flavor=:sequential), -] - -## Gauge invariant function of the environment -# -------------------------------------------- -function rho(env) - @tensor ρ[-1 -2 -3 -4 -5 -6 -7 -8] := - env.edges[WEST, 1, 1][1 -1 -2; 4] * - env.corners[NORTHWEST, 1, 1][4; 5] * - env.edges[NORTH, 1, 1][5 -3 -4; 8] * - env.corners[NORTHEAST, 1, 1][8; 9] * - env.edges[EAST, 1, 1][9 -5 -6; 12] * - env.corners[SOUTHEAST, 1, 1][12; 13] * - env.edges[SOUTH, 1, 1][13 -7 -8; 16] * - env.corners[SOUTHWEST, 1, 1][16; 1] - return ρ -end - -## Tests -# ------ -@testset "Reverse rules for ctmrg_iteration with spacetype $(Vspaces[i])" for i in - eachindex( - Pspaces -) - Random.seed!(42039482030) - psi = InfinitePEPS(Pspaces[i], Vspaces[i], Vspaces[i]) - env = CTMRGEnv(psi, Espaces[i]) - - @testset "$alg" for alg in boundary_algs - @info "$(typeof(alg)) on $(Vspaces[i])" - f(state, env) = rho(ctmrg_iteration(state, env, alg)[1]) - - # use rrule testing functionality but compute rrule via Zygote - test_rrule( - Zygote.ZygoteRuleConfig(), - f, - psi, - env; - check_inferred=false, - atol, - rrule_f=rrule_via_ad, - ) - end -end From 231cf927c2c49466e15992a0a3a0476d9aa91410 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Mon, 9 Dec 2024 17:36:30 +0100 Subject: [PATCH 12/19] Apply suggestions, update gaugefix test --- README.md | 2 +- docs/src/index.md | 2 +- examples/boundary_mps.jl | 4 +- examples/heisenberg.jl | 2 +- src/PEPSKit.jl | 30 ++- .../contractions/ctmrg_contractions.jl | 26 +-- src/algorithms/ctmrg/ctmrg.jl | 201 +++++------------- src/algorithms/ctmrg/gaugefix.jl | 34 --- src/algorithms/ctmrg/projectors.jl | 106 +++++++++ src/algorithms/ctmrg/sequential.jl | 47 +++- src/algorithms/ctmrg/simultaneous.jl | 46 +++- src/algorithms/peps_opt.jl | 19 +- test/boundarymps/vumps.jl | 8 +- test/ctmrg/fixed_iterscheme.jl | 12 +- test/ctmrg/flavors.jl | 11 +- test/ctmrg/gaugefix.jl | 46 +--- test/ctmrg/gradients.jl | 8 +- test/ctmrg/unitcell.jl | 2 +- test/heisenberg.jl | 2 +- test/j1j2_model.jl | 2 +- test/pwave.jl | 2 +- test/tf_ising.jl | 2 +- 22 files changed, 305 insertions(+), 309 deletions(-) create mode 100644 src/algorithms/ctmrg/projectors.jl diff --git a/README.md b/README.md index ebb55861..25f9a924 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/docs/src/index.md b/docs/src/index.md index 465b150f..709f4222 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -26,7 +26,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl index eb699b21..4f00915e 100644 --- a/examples/boundary_mps.jl +++ b/examples/boundary_mps.jl @@ -32,7 +32,7 @@ mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) N = abs(prod(expectation_value(mps, T))) # This can be compared to the result obtained using the CTMRG algorithm -ctm = leading_boundary(peps, CTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))) +ctm = leading_boundary(peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))) N´ = abs(norm(peps, ctm)) @show abs(N - N´) / N @@ -53,7 +53,7 @@ mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2)) mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS()) N2 = abs(prod(expectation_value(mps2, T2))) -ctm2 = leading_boundary(peps2, CTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))) +ctm2 = leading_boundary(peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))) N2´ = abs(norm(peps2, ctm2)) @show abs(N2 - N2´) / N2 diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index 25b30c02..34baed18 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -11,7 +11,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # Parameters χbond = 2 χenv = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=2) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, verbosity=2) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 0b2d11d4..334846d4 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -41,6 +41,7 @@ include("algorithms/contractions/ctmrg_contractions.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") +include("algorithms/ctmrg/projectors.jl") include("algorithms/ctmrg/simultaneous.jl") include("algorithms/ctmrg/sequential.jl") include("algorithms/ctmrg/gaugefix.jl") @@ -56,7 +57,6 @@ include("utility/symmetrization.jl") const ctmrg_maxiter = 100 const ctmrg_miniter = 4 const ctmrg_tol = 1e-8 - const ctmrg_flavor = :simultaneous const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 const reuse_env = true @@ -64,7 +64,11 @@ include("utility/symmetrization.jl") const fwd_alg = TensorKit.SDD() const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1) const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg) - const projector_alg = HalfInfiniteProjector + const projector_alg_type = HalfInfiniteProjector + const projector_alg = projector_alg_type(svd_alg, trscheme, 2) + const ctmrg_alg = SimultaneousCTMRG( + ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg + ) const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -79,7 +83,6 @@ Module containing default values that represent typical algorithm parameters. - `ctmrg_maxiter`: Maximal number of CTMRG iterations per run - `ctmrg_miniter`: Minimal number of CTMRG carried out - `ctmrg_tol`: Tolerance checking singular value and norm convergence -- `ctmrg_flavor`: Scheme for growing, projecting and renormalizing CTMRG environments - `fpgrad_maxiter`: Maximal number of iterations for computing the CTMRG fixed-point gradient - `fpgrad_tol`: Convergence tolerance for the fixed-point gradient iteration - `reuse_env`: If `true`, the current optimization step is initialized on the previous environment @@ -87,7 +90,9 @@ Module containing default values that represent typical algorithm parameters. - `fwd_alg`: SVD algorithm that is used in the forward pass - `rrule_alg`: Reverse-rule for differentiating that SVD - `svd_alg`: Combination of `fwd_alg` and `rrule_alg` +- `projector_alg_type`: Default type of projector algorithm - `projector_alg`: Algorithm to compute CTMRG projectors +- `ctmrg_alg`: Algorithm for performing CTMRG runs - `optimizer`: Optimization algorithm for PEPS ground-state optimization - `gradient_linsolver`: Default linear solver for the `LinSolver` gradient algorithm - `iterscheme`: Scheme for differentiating one CTMRG iteration @@ -96,11 +101,15 @@ Module containing default values that represent typical algorithm parameters. """ module Defaults using TensorKit, KrylovKit, OptimKit, OhMyThreads - using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint, HalfInfiniteProjector + using PEPSKit: + LinSolver, + FixedSpaceTruncation, + SVDAdjoint, + HalfInfiniteProjector, + SimultaneousCTMRG + const ctmrg_tol = 1e-8 const ctmrg_maxiter = 100 const ctmrg_miniter = 4 - const ctmrg_tol = 1e-8 - const ctmrg_flavor = :simultaneous const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 const sparse = false @@ -109,7 +118,11 @@ module Defaults const fwd_alg = TensorKit.SDD() const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1) const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg) - const projector_alg = HalfInfiniteProjector + const projector_alg_type = HalfInfiniteProjector + const projector_alg = projector_alg_type(svd_alg, trscheme, 2) + const ctmrg_alg = SimultaneousCTMRG( + ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg + ) const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -165,7 +178,8 @@ end using .Defaults: set_scheduler! export set_scheduler! export SVDAdjoint, IterSVD, NonTruncSVDAdjoint -export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector, CTMRG, CTMRGEnv +export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG +export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector export LocalOperator export expectation_value, costfun, product_peps, correlation_length export leading_boundary diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index 8e784dab..227e5a7b 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -210,10 +210,10 @@ function right_projector(E_1, C, E_2, U, isqS, ket::PEPSTensor, bra::PEPSTensor= end """ - left_and_right_projector(U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} + contract_projectors(U, S, V, Q, Q_next) -Compute left and right projectors based on a SVD and quadrant tensors, where the inverse -square root `isqS` of the singular values is computed. +Compute projectors based on a SVD of `Q * Q_next`, where the inverse square root +`isqS` of the singular values is computed. Left projector: ``` @@ -229,24 +229,10 @@ Right projector: |~~| == |~~~| == ``` """ -function left_and_right_projector( - U, S, V, Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3} -) where {E<:ElementarySpace} +function contract_projectors(U, S, V, Q, Q_next) isqS = sdiag_inv_sqrt(S) - @autoopt @tensor P_left[χ_in D_inabove D_inbelow; χ_out] := - Q_next[χ_in D_inabove D_inbelow; χ1 D1 D2] * conj(V[χ2; χ1 D1 D2]) * isqS[χ2; χ_out] - @autoopt @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := - isqS[χ_in; χ1] * conj(U[χ2 D1 D2; χ1]) * Q[χ2 D1 D2; χ_out D_outabove D_outbelow] - return P_left, P_right -end -function left_and_right_projector( - U, S, V, R_left::AbstractTensorMap{E,1,3}, L_right::AbstractTensorMap{E,3,1} -) where {E<:ElementarySpace} - isqS = sdiag_inv_sqrt(S) - @tensor P_left[χ_in D_inabove D_inbelow; χ_out] := - L_right[χ_in D_inabove D_inbelow; full] * conj(V[χ2; full]) * isqS[χ2; χ_out] - @tensor P_right[χ_in; χ_out D_outabove D_outbelow] := - isqS[χ_in; χ1] * conj(U[full; χ1]) * R_left[full; χ_out D_outabove D_outbelow] + P_left = Q_next * V' * isqS # use * to respect fermionic case + P_right = isqS * U' * Q return P_left, P_right end diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 2485be29..bfad86c0 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -1,131 +1,37 @@ -""" - FixedSpaceTruncation <: TensorKit.TruncationScheme -CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD -is performed fixed. Since different environment directions and unit cell entries might -have different spaces, this truncation style is different from `TruncationSpace`. """ -struct FixedSpaceTruncation <: TensorKit.TruncationScheme end + CTMRGAlgorithm +Abstract super type for the corner transfer matrix renormalization group (CTMRG) algorithm +for contracting infinite PEPS. """ - struct HalfInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, - trscheme=Defaults.trscheme, verbosity=0) +abstract type CTMRGAlgorithm end -Projector algorithm implementing projectors from SVDing the half-infinite CTMRG environment. """ -@kwdef struct HalfInfiniteProjector{S<:SVDAdjoint,T} - svd_alg::S = Defaults.svd_alg - trscheme::T = Defaults.trscheme - verbosity::Int = 0 -end + ctmrg_iteration(state, env, alg::CTMRG) -> env′, info +Perform a single CTMRG iteration in which all directions are being grown and renormalized. """ - struct FullInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, - trscheme=Defaults.trscheme, verbosity=0) +function ctmrg_iteration(state, env, alg::CTMRGAlgorithm) end -Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG environment. """ -@kwdef struct FullInfiniteProjector{S<:SVDAdjoint,T} - svd_alg::S = Defaults.svd_alg - trscheme::T = Defaults.trscheme - verbosity::Int = 0 -end + MPSKit.leading_boundary([envinit], state, alg::CTMRGAlgorithm) -# TODO: do AbstractProjectorAlg type instead? -> would make it easier for users to implement custom projector alg -const ProjectorAlgs = Union{HalfInfiniteProjector,FullInfiniteProjector} +Contract `state` using CTMRG and return the CTM environment. Per default, a random +initial environment is used. -function svd_algorithm(alg::ProjectorAlgs, (dir, r, c)) - if alg.svd_alg isa SVDAdjoint{<:FixedSVD} - fwd_alg = alg.svd_alg.fwd_alg - fix_svd = FixedSVD(fwd_alg.U[dir, r, c], fwd_alg.S[dir, r, c], fwd_alg.V[dir, r, c]) - return SVDAdjoint(; fwd_alg=fix_svd, rrule_alg=alg.svd_alg.rrule_alg) - else - return alg.svd_alg - end -end +Each CTMRG run is converged up to `alg.tol` where the singular value convergence +of the corners and edges is checked. The maximal and minimal number of CTMRG +iterations is set with `alg.maxiter` and `alg.miniter`. -function truncation_scheme(alg::ProjectorAlgs, Espace) - if alg.trscheme isa FixedSpaceTruncation - return truncspace(Espace) - else - return alg.trscheme - end -end - -""" - CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, flavor=Defaults.ctmrg_flavor, verbosity=0, - svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) - -Algorithm struct that represents the CTMRG algorithm for contracting infinite PEPS. -Each CTMRG run is converged up to `tol` where the singular value convergence of the -corners as well as the norm is checked. The maximal and minimal number of CTMRG iterations -is set with `maxiter` and `miniter`. - -In general, two different flavors of CTMRG can be selected with `flavor` which determine how -CTMRG is implemented. It can either be `:sequential`, where the projectors are succesively -computed on the west side, and then applied and rotated. Or with `:simultaneous` all projectors -are computed and applied simultaneously on all sides, where the corners get contracted with -two projectors at the same time. - -Different levels of output information are printed depending on `verbosity`, where `0` +Different levels of output information are printed depending on `alg.verbosity`, where `0` suppresses all output, `1` only prints warnings, `2` gives information at the start and end, and `3` prints information every iteration. - -The projectors are computed from `svd_alg` SVDs where the truncation scheme is set via -`trscheme`. """ -struct CTMRG - tol::Float64 - maxiter::Int - miniter::Int - flavor::Symbol - verbosity::Int - projector_alg::ProjectorAlgs -end -function CTMRG(; - tol=Defaults.ctmrg_tol, - maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, - flavor=Defaults.ctmrg_flavor, - verbosity=2, - projector_alg=Defaults.projector_alg, - svd_alg=Defaults.svd_alg, - trscheme=Defaults.trscheme, -) - return CTMRG( - tol, - maxiter, - miniter, - flavor, - verbosity, - projector_alg(; svd_alg, trscheme, verbosity), - ) -end - -""" - ctmrg_iteration(state, env, alg::CTMRG) - -Perform a single CTMRG iteration in which all directions are being grown and renormalized. -""" -function ctmrg_iteration(state, env, alg::CTMRG) - if alg.flavor == :simultaneous - return simultaneous_ctmrg_iter(state, env, alg) - elseif alg.flavor == :sequential - return sequential_ctmrg_iter(state, env, alg) - end -end - -""" - MPSKit.leading_boundary([envinit], state, alg::CTMRG) - -Contract `state` using CTMRG and return the CTM environment. -Per default, a random initial environment is used. -""" -function MPSKit.leading_boundary(state, alg::CTMRG) +function MPSKit.leading_boundary(state, alg::CTMRGAlgorithm) return MPSKit.leading_boundary(CTMRGEnv(state, oneunit(spacetype(state))), state, alg) end -function MPSKit.leading_boundary(envinit, state, alg::CTMRG) +function MPSKit.leading_boundary(envinit, state, alg::CTMRGAlgorithm) CS = map(x -> tsvd(x)[2], envinit.corners) TS = map(x -> tsvd(x)[2], envinit.edges) @@ -166,51 +72,42 @@ ctmrg_logcancel!(log, iter, η, N) = @warnv 1 logcancel!(log, iter, η, N) @non_differentiable ctmrg_logfinish!(args...) @non_differentiable ctmrg_logcancel!(args...) +#= +In order to compute an error measure, we compare the singular values of the current iteration with the previous one. +However, when the virtual spaces change, this comparison is not directly possible. +Instead, we project both tensors into the smaller space and then compare the difference. + +TODO: we might want to consider embedding the smaller tensor into the larger space and then compute the difference +=# +function _singular_value_distance((S₁, S₂)) + V₁ = space(S₁, 1) + V₂ = space(S₂, 1) + if V₁ == V₂ + return norm(S₁ - S₂) + else + V = infimum(V₁, V₂) + e1 = isometry(V₁, V) + e2 = isometry(V₂, V) + return norm(e1' * S₁ * e1 - e2' * S₂ * e2) + end +end + """ - compute_projector(enlarged_corners, coordinate, alg::ProjectorAlgs) + calc_convergence(envs, CSold, TSold) -Determine left and right projectors at the bond given determined by the enlarged corners -and the given coordinate using the specified `alg`. +Given a new environment `envs` and the singular values of previous corners and edges +`CSold` and `TSold`, compute the maximal singular value distance. """ -function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) - # SVD half-infinite environment - halfinf = half_infinite_environment(enlarged_corners...) - trscheme = truncation_scheme(alg, space(enlarged_corners[2], 1)) - svd_alg = svd_algorithm(alg, coordinate) - U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=trscheme) - - # Compute SVD truncation error and check for degenerate singular values - 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 +function calc_convergence(envs, CSold, TSold) + CSnew = map(x -> tsvd(x)[2], envs.corners) + ΔCS = maximum(_singular_value_distance, zip(CSold, CSnew)) - P_left, P_right = left_and_right_projector(U, S, V, enlarged_corners...) - return (P_left, P_right), (; err, U, S, V) -end -function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) - # QR left and right half-infinite environments - halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) - halfinf_right = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) - _, R_left = leftorth!(halfinf_left) - L_right, _ = rightorth!(halfinf_right) - - # SVD product of QRs - fullinf = R_left * L_right - trscheme = truncation_scheme(alg, space(enlarged_corners[4], 1)) - svd_alg = svd_algorithm(alg, coordinate) - U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=trscheme) - - # Compute SVD truncation error and check for degenerate singular values - 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 + TSnew = map(x -> tsvd(x)[2], envs.edges) + ΔTS = maximum(_singular_value_distance, zip(TSold, TSnew)) - P_left, P_right = left_and_right_projector(U, S, V, R_left, L_right) - return (P_left, P_right), (; err, U, S, V) + @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS" + + return max(ΔCS, ΔTS), CSnew, TSnew end + +@non_differentiable calc_convergence(args...) diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index 5e63d5ef..bedc784d 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -141,40 +141,6 @@ function fix_global_phases(envprev::CTMRGEnv, envfix::CTMRGEnv) return CTMRGEnv(cornersgfix, edgesgfix) end -#= -In order to compute an error measure, we compare the singular values of the current iteration with the previous one. -However, when the virtual spaces change, this comparison is not directly possible. -Instead, we project both tensors into the smaller space and then compare the difference. - -TODO: we might want to consider embedding the smaller tensor into the larger space and then compute the difference -=# -function _singular_value_distance((S₁, S₂)) - V₁ = space(S₁, 1) - V₂ = space(S₂, 1) - if V₁ == V₂ - return norm(S₁ - S₂) - else - V = infimum(V₁, V₂) - e1 = isometry(V₁, V) - e2 = isometry(V₂, V) - return norm(e1' * S₁ * e1 - e2' * S₂ * e2) - end -end - -function calc_convergence(envs, CSold, TSold) - CSnew = map(x -> tsvd(x)[2], envs.corners) - ΔCS = maximum(_singular_value_distance, zip(CSold, CSnew)) - - TSnew = map(x -> tsvd(x)[2], envs.edges) - ΔTS = maximum(_singular_value_distance, zip(TSold, TSnew)) - - @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS" - - return max(ΔCS, ΔTS), CSnew, TSnew -end - -@non_differentiable calc_convergence(args...) - """ calc_elementwise_convergence(envfinal, envfix; atol=1e-6) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl new file mode 100644 index 00000000..cf098609 --- /dev/null +++ b/src/algorithms/ctmrg/projectors.jl @@ -0,0 +1,106 @@ +""" + FixedSpaceTruncation <: TensorKit.TruncationScheme + +CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD +is performed fixed. Since different environment directions and unit cell entries might +have different spaces, this truncation style is different from `TruncationSpace`. +""" +struct FixedSpaceTruncation <: TensorKit.TruncationScheme end + +""" + ProjectorAlgorithm + +Abstract super type for all CTMRG projector algorithms. +""" +abstract type ProjectorAlgorithm end + +function svd_algorithm(alg::ProjectorAlgorithm, (dir, r, c)) + if alg.svd_alg isa SVDAdjoint{<:FixedSVD} + fwd_alg = alg.svd_alg.fwd_alg + fix_svd = FixedSVD(fwd_alg.U[dir, r, c], fwd_alg.S[dir, r, c], fwd_alg.V[dir, r, c]) + return SVDAdjoint(; fwd_alg=fix_svd, rrule_alg=alg.svd_alg.rrule_alg) + else + return alg.svd_alg + end +end + +function truncation_scheme(alg::ProjectorAlgorithm, Espace) + if alg.trscheme isa FixedSpaceTruncation + return truncspace(Espace) + else + return alg.trscheme + end +end + +""" + struct HalfInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the half-infinite CTMRG environment. +""" +@kwdef struct HalfInfiniteProjector{S<:SVDAdjoint,T} <: ProjectorAlgorithm + svd_alg::S = Defaults.svd_alg + trscheme::T = Defaults.trscheme + verbosity::Int = 0 +end + +""" + struct FullInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG environment. +""" +@kwdef struct FullInfiniteProjector{S<:SVDAdjoint,T} <: ProjectorAlgorithm + svd_alg::S = Defaults.svd_alg + trscheme::T = Defaults.trscheme + verbosity::Int = 0 +end + +""" + compute_projector(enlarged_corners, coordinate, alg::ProjectorAlgorithm) + +Determine left and right projectors at the bond given determined by the enlarged corners +and the given coordinate using the specified `alg`. +""" +function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) + # SVD half-infinite environment + halfinf = half_infinite_environment(enlarged_corners...) + trscheme = truncation_scheme(alg, space(enlarged_corners[2], 1)) + svd_alg = svd_algorithm(alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=trscheme) + + # Compute SVD truncation error and check for degenerate singular values + 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, enlarged_corners...) + return (P_left, P_right), (; err, U, S, V) +end +function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) + # QR left and right half-infinite environments + halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) + halfinf_right = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) + _, R_left = leftorth!(halfinf_left) + L_right, _ = rightorth!(halfinf_right) + + # SVD product of QRs + fullinf = R_left * L_right + trscheme = truncation_scheme(alg, space(enlarged_corners[4], 1)) + svd_alg = svd_algorithm(alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=trscheme) + + # Compute SVD truncation error and check for degenerate singular values + 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, R_left, L_right) + return (P_left, P_right), (; err, U, S, V) +end \ No newline at end of file diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 60c27a4d..bcb147e6 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -1,10 +1,36 @@ """ - sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info + SequentialCTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, verbosity=0, + projector_alg=typeof(Defaults.projector_alg), + svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) -Perform one sequential iteration of CTMRG, where one iteration consists of four expansion, -renormalization and rotation steps that are performed sequentially. +CTMRG algorithm where the expansions and renormalization is performed sequentially +column-wise. This is implemented as a growing and projecting step to the left, followed by +a clockwise rotation (performed four times). The projectors are computed using +`projector_alg` from `svd_alg` SVDs where the truncation scheme is set via `trscheme`. """ -function sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) +struct SequentialCTMRG <: CTMRGAlgorithm + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::ProjectorAlgorithm +end +function SequentialCTMRG(; + tol=Defaults.ctmrg_tol, + maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, + verbosity=2, + projector_alg=Defaults.projector_alg_type, + svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, +) + return SequentialCTMRG( + tol, maxiter, miniter, verbosity, projector_alg(; svd_alg, trscheme, verbosity) + ) +end + +function ctmrg_iteration(state, envs::CTMRGEnv, alg::SequentialCTMRG) ϵ = zero(real(scalartype(state))) for _ in 1:4 # rotate for col in 1:size(state, 2) # left move column-wise @@ -20,26 +46,25 @@ function sequential_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) end """ - sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) - sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs) + sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm) + sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm) Compute CTMRG projectors in the `:sequential` scheme either for an entire column `col` or for a specific `coordinate` (where `dir=WEST` is already implied in the `:sequential` scheme). """ function sequential_projectors( - col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgs + col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm ) - ϵ = zero(real(scalartype(envs))) - # SVD half-infinite environment column-wise + ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs, 2))) coordinates = eachcoordinate(envs)[:, col] projectors = dtmap(coordinates) do (r, c) proj, info = sequential_projectors((WEST, r, c), state, envs, alg) - ϵ = max(ϵ, info.err / norm(info.S)) + ϵ[c] = max(ϵ, info.err / norm(info.S)) return proj end - return (map(first, projectors), map(last, projectors)), (; err=ϵ) + return (map(first, projectors), map(last, projectors)), (; err=maximum(copy(ϵ))) end function sequential_projectors( coordinate::NTuple{3,Int}, diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index b6992fbb..75d1e72f 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -1,10 +1,36 @@ """ - simultaneous_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info + SimultaneousCTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, verbosity=0, + projector_alg=Defaults.projector_alg, + svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) -Perform one simultaneous iteration of CTMRG, in which the environment is expanded and -renormalized in all directions at the same time. +CTMRG algorithm where all sides are grown and renormalized at the same time. In particular, +the projectors are applied to the corners from two sides simultaneously. The projectors are +computed using `projector_alg` from `svd_alg` SVDs where the truncation scheme is set via +`trscheme`. """ -function simultaneous_ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) +struct SimultaneousCTMRG <: CTMRGAlgorithm + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::ProjectorAlgorithm +end +function SimultaneousCTMRG(; + tol=Defaults.ctmrg_tol, + maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, + verbosity=2, + projector_alg=Defaults.projector_alg_type, + svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, +) + return SimultaneousCTMRG( + tol, maxiter, miniter, verbosity, projector_alg(; svd_alg, trscheme, verbosity) + ) +end + +function ctmrg_iteration(state, envs::CTMRGEnv, alg::SimultaneousCTMRG) enlarged_corners = dtmap(eachcoordinate(state, 1:4)) do idx return TensorMap(EnlargedCorner(state, envs, idx), idx[1]) end # expand environment @@ -31,30 +57,30 @@ function _prealloc_svd(edges::Array{E,N}, ::FullInfiniteProjector) where {E,N} end """ - simultaneous_projectors(enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgs) - simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgs) + simultaneous_projectors(enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm) + simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgorithm) Compute CTMRG projectors in the `:simultaneous` scheme either for all provided enlarged corners or on a specific `coordinate`. """ function simultaneous_projectors( - enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgs + enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm ) where {E} U, S, V = _prealloc_svd(envs.edges, alg) - ϵ = zero(real(scalartype(envs))) + ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate proj, info = simultaneous_projectors(coordinate, enlarged_corners, alg) U[coordinate...] = info.U S[coordinate...] = info.S V[coordinate...] = info.V - ϵ = max(ϵ, info.err / norm(info.S)) + ϵ[coordinate...] = info.err / norm(info.S) return proj end P_left = map(first, projectors) P_right = map(last, projectors) - return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) + return (P_left, P_right), (; err=maximum(copy(ϵ)), U=copy(U), S=copy(S), V=copy(V)) end function simultaneous_projectors( coordinate, enlarged_corners::Array{E,3}, alg::HalfInfiniteProjector diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 91e5adc8..83ffb7ac 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -77,8 +77,8 @@ function LinSolver(; end """ - PEPSOptimize{G}(; boundary_alg=CTMRG(), optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer - reuse_env::Bool=true, gradient_alg::G=LinSolver()) + PEPSOptimize{G}(; boundary_alg=Defaults.ctmrg_alg, optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer + reuse_env::Bool=true, gradient_alg::G=Defaults.gradient_alg) Algorithm struct that represent PEPS ground-state optimization using AD. Set the algorithm to contract the infinite PEPS in `boundary_alg`; @@ -89,19 +89,19 @@ step by setting `reuse_env` to true. Otherwise a random environment is used at e step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm. """ struct PEPSOptimize{G} - boundary_alg::CTMRG + boundary_alg::CTMRGAlgorithm optimizer::OptimKit.OptimizationAlgorithm reuse_env::Bool gradient_alg::G function PEPSOptimize( # Inner constructor to prohibit illegal setting combinations - boundary_alg::CTMRG, + boundary_alg::CTMRGAlgorithm, optimizer, reuse_env, gradient_alg::G, ) where {G} if gradient_alg isa GradMode - if boundary_alg.flavor === :sequential && iterscheme(gradient_alg) === :fixed + if boundary_alg isa SequentialCTMRG && iterscheme(gradient_alg) === :fixed throw(ArgumentError(":sequential and :fixed are not compatible")) end end @@ -109,7 +109,7 @@ struct PEPSOptimize{G} end end function PEPSOptimize(; - boundary_alg=CTMRG(), + boundary_alg=Defaults.ctmrg_alg, optimizer=Defaults.optimizer, reuse_env=Defaults.reuse_env, gradient_alg=Defaults.gradient_alg, @@ -220,7 +220,7 @@ function _rrule( ::typeof(MPSKit.leading_boundary), envinit, state, - alg::CTMRG, + alg::CTMRGAlgorithm, ) envs = leading_boundary(envinit, state, alg) @@ -249,9 +249,8 @@ function _rrule( ::typeof(MPSKit.leading_boundary), envinit, state, - alg::CTMRG, + alg::SimultaneousCTMRG, ) - @assert alg.flavor === :simultaneous @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg) envs = leading_boundary(envinit, state, alg) envsconv, info = ctmrg_iteration(state, envs, alg) @@ -262,7 +261,7 @@ function _rrule( svd_alg_fixed = SVDAdjoint(; fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) - alg_fixed = CTMRG(; svd_alg=svd_alg_fixed, trscheme=notrunc(), flavor=:simultaneous) + alg_fixed = SimultaneousCTMRG(; svd_alg=svd_alg_fixed, trscheme=notrunc()) function leading_boundary_fixed_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl index a161c073..2ec6265c 100644 --- a/test/boundarymps/vumps.jl +++ b/test/boundarymps/vumps.jl @@ -17,7 +17,9 @@ const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitia mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(sum(expectation_value(mps, T))) - ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) + ctm = leading_boundary( + CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1) + ) N´ = abs(norm(psi, ctm)) @test N ≈ N´ atol = 1e-3 @@ -31,7 +33,9 @@ end mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(prod(expectation_value(mps, T))) - ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) + ctm = leading_boundary( + CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1) + ) N´ = abs(norm(psi, ctm)) @test N ≈ N´ rtol = 1e-2 diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl index f13b04a6..c0062abf 100644 --- a/test/ctmrg/fixed_iterscheme.jl +++ b/test/ctmrg/fixed_iterscheme.jl @@ -25,7 +25,7 @@ atol = 1e-5 ) in Iterators.product( unitcells, svd_algs, projector_algs ) - ctm_alg = CTMRG(; svd_alg, projector_alg) + ctm_alg = SimultaneousCTMRG(; svd_alg, projector_alg) # initialize states Random.seed!(2394823842) @@ -40,7 +40,7 @@ atol = 1e-5 # fix gauge of SVD U_fix, V_fix = fix_relative_phases(info.U, info.V, signs) svd_alg_fix = SVDAdjoint(; fwd_alg=FixedSVD(U_fix, info.S, V_fix)) - ctm_alg_fix = CTMRG(; projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc()) + ctm_alg_fix = SimultaneousCTMRG(; projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc()) # do iteration with FixedSVD env_fixedsvd, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix) @@ -49,11 +49,11 @@ atol = 1e-5 end @testset "Element-wise consistency of TensorKit.SDD and IterSVD" begin - ctm_alg_iter = CTMRG(; + ctm_alg_iter = SimultaneousCTMRG(; maxiter=200, svd_alg=SVDAdjoint(; fwd_alg=IterSVD(; alg=GKL(; tol=1e-14, krylovdim=χenv + 10))), ) - ctm_alg_full = CTMRG(; svd_alg=SVDAdjoint(; fwd_alg=TensorKit.SDD())) + ctm_alg_full = SimultaneousCTMRG(; svd_alg=SVDAdjoint(; fwd_alg=TensorKit.SDD())) # initialize states Random.seed!(91283219347) @@ -73,11 +73,11 @@ end # fix gauge of SVD U_fix_iter, V_fix_iter = fix_relative_phases(info_iter.U, info_iter.V, signs_iter) svd_alg_fix_iter = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_iter, info_iter.S, V_fix_iter)) - ctm_alg_fix_iter = CTMRG(; svd_alg=svd_alg_fix_iter, trscheme=notrunc()) + ctm_alg_fix_iter = SimultaneousCTMRG(; svd_alg=svd_alg_fix_iter, trscheme=notrunc()) U_fix_full, V_fix_full = fix_relative_phases(info_full.U, info_full.V, signs_full) svd_alg_fix_full = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_full, info_full.S, V_fix_full)) - ctm_alg_fix_full = CTMRG(; svd_alg=svd_alg_fix_full, trscheme=notrunc()) + ctm_alg_fix_full = SimultaneousCTMRG(; svd_alg=svd_alg_fix_full, trscheme=notrunc()) # do iteration with FixedSVD env_fixedsvd_iter, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_iter) diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl index 9337bc08..3dd26156 100644 --- a/test/ctmrg/flavors.jl +++ b/test/ctmrg/flavors.jl @@ -7,8 +7,8 @@ using PEPSKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg_sequential = CTMRG(; flavor=:sequential) -ctm_alg_simultaneous = CTMRG(; flavor=:simultaneous) +ctm_alg_sequential = SequentialCTMRG() +ctm_alg_simultaneous = SimultaneousCTMRG() unitcells = [(1, 1), (3, 4)] projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] @@ -56,10 +56,9 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] end # test fixedspace actually fixes space -@testset "Fixedspace truncation ($flavor)" for flavor in [:sequential, :simultaneous] - ctm_alg = CTMRG(; - tol=1e-6, maxiter=1, verbosity=0, flavor, trscheme=FixedSpaceTruncation() - ) +@testset "Fixedspace truncation ($ctmrg_alg)" for ctmrg_alg in + [SequentialCTMRG, SimultanenousCTMRG] + ctm_alg = ctmrg_alg(; tol=1e-6, maxiter=1, verbosity=0, trscheme=FixedSpaceTruncation()) Ds = fill(2, 3, 3) χs = [16 17 18; 15 20 21; 14 19 22] psi = InfinitePEPS(Ds, Ds, Ds) diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 040d72ed..bfb39fef 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -8,35 +8,13 @@ using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] maxiter = 400 -ctmrg_flavors = [:simultaneous, :sequential] +ctmrg_flavors = [SequentialCTMRG, SimultaneousCTMRG] projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] χ = 6 atol = 1e-4 -function _make_symmetric!(psi) - if ==(size(psi)...) - return symmetrize!(psi, RotateReflect()) - else - return symmetrize!(symmetrize!(psi, ReflectDepth()), ReflectWidth()) - end -end - -# If I can't make the rng seed behave, I'll just randomly define a peps somehow -function _semi_random_peps!(psi::InfinitePEPS) - i = 0 - A′ = map(psi.A) do a - for (_, b) in blocks(a) - l = length(b) - b .= reshape(collect((1:l) .+ i), size(b)) - i += l - end - return a - end - return InfinitePEPS(A′) -end - -@testset "Trivial symmetry ($T) - ($unitcell) - ($flavor) - ($projector_alg)" for ( - T, unitcell, flavor, projector_alg +@testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( + T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( scalartypes, unitcells, ctmrg_flavors, projector_algs ) @@ -44,13 +22,11 @@ end peps_space = ComplexSpace(2) ctm_space = ComplexSpace(χ) - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) + psi = InfinitePEPS(randn, T, physical_space, peps_space; unitcell) Random.seed!(987654321) # Seed RNG to make random environment consistent env = CTMRGEnv(psi, ctm_space) - alg = CTMRG(; maxiter, flavor, projector_alg) + alg = ctmrg_alg(; maxiter, projector_alg) env = leading_boundary(env, psi, alg) env′, = ctmrg_iteration(psi, env, alg) @@ -58,23 +34,21 @@ end @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol end -@testset "Z2 symmetry ($T) - ($unitcell) - ($flavor) - ($projector_alg)" for ( - T, unitcell, flavor, projector_alg +@testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( + T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( - scalartypes, unitcells, ctmrg_flavors, projector_algs + scalartypes, unitcells, ctmrg_s, projector_algs ) physical_space = Z2Space(0 => 1, 1 => 1) peps_space = Z2Space(0 => 1, 1 => 1) ctm_space = Z2Space(0 => χ ÷ 2, 1 => χ ÷ 2) - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) + psi = InfinitePEPS(randn, T, physical_space, peps_space; unitcell) Random.seed!(987654321) # Seed RNG to make random environment consistent psi = InfinitePEPS(physical_space, peps_space; unitcell) env = CTMRGEnv(psi, ctm_space) - alg = CTMRG(; maxiter, flavor, projector_alg) + alg = ctmrg_alg(; maxiter, projector_alg) env = leading_boundary(env, psi, alg) env′, = ctmrg_iteration(psi, env, alg) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index a4ceafea..93c278ce 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -18,10 +18,10 @@ names = ["Heisenberg", "p-wave superconductor"] gradtol = 1e-4 boundary_algs = [ - # CTMRG(; verbosity=0, flavor=:simultaneous, projector_alg=HalfInfiniteProjector), - CTMRG(; verbosity=0, flavor=:simultaneous, projector_alg=FullInfiniteProjector), - # CTMRG(; verbosity=0, flavor=:sequential, projector_alg=HalfInfiniteProjector), - CTMRG(; verbosity=0, flavor=:sequential, projector_alg=FullInfiniteProjector), + SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SequentialCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), + SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), ] gradmodes = [ [ diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 1ba2ea40..9601539b 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -7,7 +7,7 @@ using TensorKit # settings Random.seed!(91283219347) stype = ComplexF64 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() function test_unitcell( unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 1c05b084..d985503a 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) diff --git a/test/j1j2_model.jl b/test/j1j2_model.jl index 45acd5b1..8dd88397 100644 --- a/test/j1j2_model.jl +++ b/test/j1j2_model.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 12 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2), diff --git a/test/pwave.jl b/test/pwave.jl index d663ec0d..dcdecf7e 100644 --- a/test/pwave.jl +++ b/test/pwave.jl @@ -10,7 +10,7 @@ unitcell = (2, 2) H = pwave_superconductor(InfiniteSquare(unitcell...)) χbond = 2 χenv = 16 -ctm_alg = CTMRG(; maxiter=150) +ctm_alg = SimultaneousCTMRG(; maxiter=150) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=2) ) diff --git a/test/tf_ising.jl b/test/tf_ising.jl index 0b4d72f5..8b84696f 100644 --- a/test/tf_ising.jl +++ b/test/tf_ising.jl @@ -19,7 +19,7 @@ mˣ = 0.91 # initialize parameters χbond = 2 χenv = 16 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) From 48e41c3a03031efb1fa4b746d66db4d1489af125 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Mon, 9 Dec 2024 19:27:57 +0100 Subject: [PATCH 13/19] Fix FixedSpaceTruncation --- src/PEPSKit.jl | 2 +- src/algorithms/ctmrg/projectors.jl | 10 +++---- src/algorithms/ctmrg/sequential.jl | 10 +++---- src/algorithms/ctmrg/simultaneous.jl | 6 +++- test/ctmrg/flavors.jl | 15 ++++++++-- test/ctmrg/gaugefix.jl | 6 ++-- test/ctmrg/unitcell.jl | 45 +++++++++++++++++++++++----- test/runtests.jl | 2 +- 8 files changed, 68 insertions(+), 28 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 334846d4..4fbf6760 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -3,7 +3,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Base: @kwdef using Compat -using Accessors +using Accessors: @set using VectorInterface using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations using ChainRulesCore, Zygote diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index cf098609..a02a896c 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -24,9 +24,9 @@ function svd_algorithm(alg::ProjectorAlgorithm, (dir, r, c)) end end -function truncation_scheme(alg::ProjectorAlgorithm, Espace) +function truncation_scheme(alg::ProjectorAlgorithm, edge) if alg.trscheme isa FixedSpaceTruncation - return truncspace(Espace) + return truncspace(space(edge, 1)) else return alg.trscheme end @@ -65,9 +65,8 @@ and the given coordinate using the specified `alg`. function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) # SVD half-infinite environment halfinf = half_infinite_environment(enlarged_corners...) - trscheme = truncation_scheme(alg, space(enlarged_corners[2], 1)) svd_alg = svd_algorithm(alg, coordinate) - U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=trscheme) + U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=alg.trscheme) # Compute SVD truncation error and check for degenerate singular values Zygote.isderiving() && ignore_derivatives() do @@ -89,9 +88,8 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec # SVD product of QRs fullinf = R_left * L_right - trscheme = truncation_scheme(alg, space(enlarged_corners[4], 1)) svd_alg = svd_algorithm(alg, coordinate) - U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=trscheme) + U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=alg.trscheme) # Compute SVD truncation error and check for degenerate singular values Zygote.isderiving() && ignore_derivatives() do diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index bcb147e6..94f4846e 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -59,8 +59,11 @@ function sequential_projectors( ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs, 2))) coordinates = eachcoordinate(envs)[:, col] projectors = dtmap(coordinates) do (r, c) - proj, info = sequential_projectors((WEST, r, c), state, envs, alg) - ϵ[c] = max(ϵ, info.err / norm(info.S)) + trscheme = truncation_scheme(alg, envs.edges[WEST, _prev(r, size(envs, 2)), c]) + proj, info = sequential_projectors( + (WEST, r, c), state, envs, @set(alg.trscheme = trscheme) + ) + ϵ[r] = info.err / norm(info.S) return proj end @@ -84,9 +87,6 @@ function sequential_projectors( envs::CTMRGEnv, alg::FullInfiniteProjector, ) - # _, r, c = coordinate - # r′ = _next(r, size(envs, 2)) - # c′ = _next(c, size(envs, 3)) rowsize, colsize = size(envs)[2:3] coordinate_nw = _next_coordinate(coordinate, rowsize, colsize) coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize) diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl index 75d1e72f..fcae797b 100644 --- a/src/algorithms/ctmrg/simultaneous.jl +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -70,7 +70,11 @@ function simultaneous_projectors( ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs))) projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate - proj, info = simultaneous_projectors(coordinate, enlarged_corners, alg) + coordinate′ = _next_coordinate(coordinate, size(envs)[2:3]...) + trscheme = truncation_scheme(alg, envs.edges[coordinate[1], coordinate′[2:3]...]) + proj, info = simultaneous_projectors( + coordinate, enlarged_corners, @set(alg.trscheme = trscheme) + ) U[coordinate...] = info.U S[coordinate...] = info.S V[coordinate...] = info.V diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl index 3dd26156..dc489e9d 100644 --- a/test/ctmrg/flavors.jl +++ b/test/ctmrg/flavors.jl @@ -56,9 +56,18 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] end # test fixedspace actually fixes space -@testset "Fixedspace truncation ($ctmrg_alg)" for ctmrg_alg in - [SequentialCTMRG, SimultanenousCTMRG] - ctm_alg = ctmrg_alg(; tol=1e-6, maxiter=1, verbosity=0, trscheme=FixedSpaceTruncation()) +@testset "Fixedspace truncation using $ctmrg_alg and $projector_alg" for ( + ctmrg_alg, projector_alg +) in Iterators.product( + [SequentialCTMRG, SimultaneousCTMRG], projector_algs +) + ctm_alg = ctmrg_alg(; + tol=1e-6, + maxiter=1, + verbosity=0, + trscheme=FixedSpaceTruncation(), + projector_alg, + ) Ds = fill(2, 3, 3) χs = [16 17 18; 15 20 21; 14 19 22] psi = InfinitePEPS(Ds, Ds, Ds) diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index bfb39fef..8ffed116 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -8,7 +8,7 @@ using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] maxiter = 400 -ctmrg_flavors = [SequentialCTMRG, SimultaneousCTMRG] +ctmrg_algs = [SequentialCTMRG, SimultaneousCTMRG] projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] χ = 6 atol = 1e-4 @@ -16,7 +16,7 @@ atol = 1e-4 @testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( - scalartypes, unitcells, ctmrg_flavors, projector_algs + scalartypes, unitcells, ctmrg_algs, projector_algs ) physical_space = ComplexSpace(2) peps_space = ComplexSpace(2) @@ -37,7 +37,7 @@ end @testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( - scalartypes, unitcells, ctmrg_s, projector_algs + scalartypes, unitcells, ctmrg_algs, projector_algs ) physical_space = Z2Space(0 => 1, 1 => 1) peps_space = Z2Space(0 => 1, 1 => 1) diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 9601539b..e2f010fa 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -7,10 +7,23 @@ using TensorKit # settings Random.seed!(91283219347) stype = ComplexF64 -ctm_alg = SimultaneousCTMRG() +ctm_algs = [ + SequentialCTMRG(; projector_alg=HalfInfiniteProjector), + SequentialCTMRG(; projector_alg=FullInfiniteProjector), + SimultaneousCTMRG(; projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; projector_alg=FullInfiniteProjector), +] function test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) peps = InfinitePEPS(randn, stype, Pspaces, Nspaces, Espaces) env = CTMRGEnv(randn, stype, peps, chis_north, chis_east, chis_south, chis_west) @@ -40,7 +53,7 @@ function random_dualize!(M::AbstractMatrix{<:ElementarySpace}) return M end -@testset "Integer space specifiers" begin +@testset "Integer space specifiers with $ctm_alg" begin unitcell = (3, 3) Pspaces = rand(2:3, unitcell...) @@ -52,11 +65,19 @@ end chis_west = rand(5:10, unitcell...) test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) end -@testset "Random Cartesian spaces" begin +@testset "Random Cartesian spaces with $ctm_alg" for ctm_alg in ctm_algs unitcell = (3, 3) Pspaces = random_dualize!(ComplexSpace.(rand(2:3, unitcell...))) @@ -68,11 +89,19 @@ end chis_west = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) end -@testset "Specific U1 spaces" begin +@testset "Specific U1 spaces with $ctm_alg" for ctm_alg in ctm_algs unitcell = (2, 2) PA = U1Space(-1 => 1, 0 => 1) @@ -84,5 +113,5 @@ end Nspaces = [Vpeps Vpeps'; Vpeps' Vpeps] chis = [Venv Venv; Venv Venv] - test_unitcell(unitcell, Pspaces, Nspaces, Nspaces, chis, chis, chis, chis) + test_unitcell(ctm_alg, unitcell, Pspaces, Nspaces, Nspaces, chis, chis, chis, chis) end diff --git a/test/runtests.jl b/test/runtests.jl index 49d2fc0c..318a3f41 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,7 +26,7 @@ end @time @safetestset "Unit cells" begin include("ctmrg/unitcell.jl") end - @time @safetestset "CTMRG schemes" begin + @time @safetestset "Flavors" begin include("ctmrg/flavors.jl") end end From ea998b1f9b69b1b53cd13738e3e16653b7a9944a Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Mon, 9 Dec 2024 20:11:03 +0100 Subject: [PATCH 14/19] Format, fix gauge test, further steps towards FullInf differentiability --- examples/boundary_mps.jl | 8 ++++++-- src/algorithms/ctmrg/projectors.jl | 2 +- src/algorithms/peps_opt.jl | 3 ++- test/ctmrg/fixed_iterscheme.jl | 4 +++- test/ctmrg/flavors.jl | 6 +----- test/ctmrg/gaugefix.jl | 4 ++-- test/ctmrg/gradients.jl | 29 +++++++++++++++++------------ test/heisenberg.jl | 2 +- 8 files changed, 33 insertions(+), 25 deletions(-) diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl index 4f00915e..51af6dea 100644 --- a/examples/boundary_mps.jl +++ b/examples/boundary_mps.jl @@ -32,7 +32,9 @@ mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) N = abs(prod(expectation_value(mps, T))) # This can be compared to the result obtained using the CTMRG algorithm -ctm = leading_boundary(peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))) +ctm = leading_boundary( + peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20)) +) N´ = abs(norm(peps, ctm)) @show abs(N - N´) / N @@ -53,7 +55,9 @@ mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2)) mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS()) N2 = abs(prod(expectation_value(mps2, T2))) -ctm2 = leading_boundary(peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))) +ctm2 = leading_boundary( + peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20)) +) N2´ = abs(norm(peps2, ctm2)) @show abs(N2 - N2´) / N2 diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl index a02a896c..7d7ca4aa 100644 --- a/src/algorithms/ctmrg/projectors.jl +++ b/src/algorithms/ctmrg/projectors.jl @@ -101,4 +101,4 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec P_left, P_right = contract_projectors(U, S, V, R_left, L_right) return (P_left, P_right), (; err, U, S, V) -end \ No newline at end of file +end diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 83ffb7ac..69f17338 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -261,7 +261,8 @@ function _rrule( svd_alg_fixed = SVDAdjoint(; fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) - alg_fixed = SimultaneousCTMRG(; svd_alg=svd_alg_fixed, trscheme=notrunc()) + alg_fixed = @set alg.projector_alg.svd_alg = svd_alg_fixed + alg_fixed = @set alg.projector_alg.trscheme = notrunc() function leading_boundary_fixed_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl index c0062abf..208fdece 100644 --- a/test/ctmrg/fixed_iterscheme.jl +++ b/test/ctmrg/fixed_iterscheme.jl @@ -40,7 +40,9 @@ atol = 1e-5 # fix gauge of SVD U_fix, V_fix = fix_relative_phases(info.U, info.V, signs) svd_alg_fix = SVDAdjoint(; fwd_alg=FixedSVD(U_fix, info.S, V_fix)) - ctm_alg_fix = SimultaneousCTMRG(; projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc()) + ctm_alg_fix = SimultaneousCTMRG(; + projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc() + ) # do iteration with FixedSVD env_fixedsvd, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix) diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl index dc489e9d..5d8b440c 100644 --- a/test/ctmrg/flavors.jl +++ b/test/ctmrg/flavors.jl @@ -62,11 +62,7 @@ end [SequentialCTMRG, SimultaneousCTMRG], projector_algs ) ctm_alg = ctmrg_alg(; - tol=1e-6, - maxiter=1, - verbosity=0, - trscheme=FixedSpaceTruncation(), - projector_alg, + tol=1e-6, maxiter=1, verbosity=0, trscheme=FixedSpaceTruncation(), projector_alg ) Ds = fill(2, 3, 3) χs = [16 17 18; 15 20 21; 14 19 22] diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 8ffed116..de798ad4 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -10,7 +10,7 @@ unitcells = [(1, 1), (2, 2), (3, 2)] maxiter = 400 ctmrg_algs = [SequentialCTMRG, SimultaneousCTMRG] projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] -χ = 6 +χ = 8 atol = 1e-4 @testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( @@ -45,7 +45,7 @@ end psi = InfinitePEPS(randn, T, physical_space, peps_space; unitcell) - Random.seed!(987654321) # Seed RNG to make random environment consistent + Random.seed!(29385293852) # Seed RNG to make random environment consistent psi = InfinitePEPS(physical_space, peps_space; unitcell) env = CTMRGEnv(psi, ctm_space) alg = ctmrg_alg(; maxiter, projector_alg) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 93c278ce..1e805f9b 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -17,11 +17,17 @@ models = [heisenberg_XYZ(InfiniteSquare()), pwave_superconductor(InfiniteSquare( names = ["Heisenberg", "p-wave superconductor"] gradtol = 1e-4 -boundary_algs = [ - SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), - SequentialCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), - SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), - SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), +ctmrg_algs = [ + [ + SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), + ], + [ + SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), + SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SequentialCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), + ], ] gradmodes = [ [ @@ -52,14 +58,15 @@ steps = -0.01:0.005:0.01 Vspace = Pspaces[i] Espace = Espaces[i] gms = gradmodes[i] - boundary_alg = boundary_algs[i] + calgs = ctmrg_algs[i] psi_init = InfinitePEPS(Pspace, Vspace, Vspace) - @testset "$alg_rrule" for alg_rrule in gms - @info "optimtest of $alg_rrule on $(names[i])" + @testset "$ctmrg_alg and $alg_rrule" for (ctmrg_alg, alg_rrule) in + Iterators.product(calgs, gms) + @info "optimtest of $ctmrg_alg and $alg_rrule on $(names[i])" Random.seed!(42039482030) dir = InfinitePEPS(Pspace, Vspace, Vspace) psi = InfinitePEPS(Pspace, Vspace, Vspace) - env = leading_boundary(CTMRGEnv(psi, Espace), psi, boundary_alg) + env = leading_boundary(CTMRGEnv(psi, Espace), psi, ctmrg_alg) alphas, fs, dfs1, dfs2 = OptimKit.optimtest( (psi, env), dir; @@ -68,9 +75,7 @@ steps = -0.01:0.005:0.01 inner=PEPSKit.real_inner, ) do (peps, envs) E, g = Zygote.withgradient(peps) do psi - envs2 = PEPSKit.hook_pullback( - leading_boundary, envs, psi, boundary_alg; alg_rrule - ) + envs2 = PEPSKit.hook_pullback(leading_boundary, envs, psi, ctmrg_alg; alg_rrule) return costfun(psi, envs2, models[i]) end diff --git a/test/heisenberg.jl b/test/heisenberg.jl index d985503a..f72a31a9 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg = SimultaneousCTMRG() +ctm_alg = SimultaneousCTMRG(; projector_alg=FullInfiniteProjector) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) From 2082302af840cd37d2c277ed0788f1658370afc1 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 10 Dec 2024 13:59:48 +0100 Subject: [PATCH 15/19] Fix differentiability and gradients test --- src/algorithms/peps_opt.jl | 2 +- test/ctmrg/gaugefix.jl | 30 ++++++++++++++++++++++++++++-- test/ctmrg/gradients.jl | 11 ++++------- test/heisenberg.jl | 2 +- 4 files changed, 34 insertions(+), 11 deletions(-) diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 69f17338..a4a6cc96 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -262,7 +262,7 @@ function _rrule( fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) alg_fixed = @set alg.projector_alg.svd_alg = svd_alg_fixed - alg_fixed = @set alg.projector_alg.trscheme = notrunc() + alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc() function leading_boundary_fixed_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index de798ad4..7c2aae1e 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -13,6 +13,28 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] χ = 8 atol = 1e-4 +function _make_symmetric!(psi) + if ==(size(psi)...) + return symmetrize!(psi, RotateReflect()) + else + return symmetrize!(symmetrize!(psi, ReflectDepth()), ReflectWidth()) + end +end + +# If I can't make the rng seed behave, I'll just randomly define a peps somehow +function _semi_random_peps!(psi::InfinitePEPS) + i = 0 + A′ = map(psi.A) do a + for (_, b) in blocks(a) + l = length(b) + b .= reshape(collect((1:l) .+ i), size(b)) + i += l + end + return a + end + return InfinitePEPS(A′) +end + @testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( @@ -22,7 +44,9 @@ atol = 1e-4 peps_space = ComplexSpace(2) ctm_space = ComplexSpace(χ) - psi = InfinitePEPS(randn, T, physical_space, peps_space; unitcell) + psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) + _semi_random_peps!(psi) + _make_symmetric!(psi) Random.seed!(987654321) # Seed RNG to make random environment consistent env = CTMRGEnv(psi, ctm_space) @@ -43,7 +67,9 @@ end peps_space = Z2Space(0 => 1, 1 => 1) ctm_space = Z2Space(0 => χ ÷ 2, 1 => χ ÷ 2) - psi = InfinitePEPS(randn, T, physical_space, peps_space; unitcell) + psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) + _semi_random_peps!(psi) + _make_symmetric!(psi) Random.seed!(29385293852) # Seed RNG to make random environment consistent psi = InfinitePEPS(physical_space, peps_space; unitcell) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 1e805f9b..f39ea6cc 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -23,10 +23,7 @@ ctmrg_algs = [ SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), ], [ - SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), - SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), - SequentialCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), ], ] gradmodes = [ @@ -50,10 +47,10 @@ steps = -0.01:0.005:0.01 ## Tests # ------ -@testset "AD CTMRG energy gradients for $(names[i]) model" verbose = true for i in - eachindex( - models -) +@testset "AD CTMRG energy gradients for $(names[i]) model" verbose = true for i in [2] +# eachindex( +# models +# ) Pspace = Pspaces[i] Vspace = Pspaces[i] Espace = Espaces[i] diff --git a/test/heisenberg.jl b/test/heisenberg.jl index f72a31a9..d985503a 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg = SimultaneousCTMRG(; projector_alg=FullInfiniteProjector) +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) From 663bf84dad5e94ce926973ab2a81bf4903552e0f Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 10 Dec 2024 14:00:28 +0100 Subject: [PATCH 16/19] Restore gradients test --- test/ctmrg/gradients.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index f39ea6cc..648d48bd 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -22,9 +22,7 @@ ctmrg_algs = [ SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), ], - [ - SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), - ], + [SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector)], ] gradmodes = [ [ @@ -47,10 +45,10 @@ steps = -0.01:0.005:0.01 ## Tests # ------ -@testset "AD CTMRG energy gradients for $(names[i]) model" verbose = true for i in [2] -# eachindex( -# models -# ) +@testset "AD CTMRG energy gradients for $(names[i]) model" verbose = true for i in + eachindex( + models +) Pspace = Pspaces[i] Vspace = Pspaces[i] Espace = Espaces[i] From 45456f008b9cd4348ea97ecc1c6a742814e31ced Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 10 Dec 2024 15:13:12 +0100 Subject: [PATCH 17/19] Update gaugefix test again --- test/ctmrg/gaugefix.jl | 86 +++++++++++++----------------------------- 1 file changed, 27 insertions(+), 59 deletions(-) diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index 7c2aae1e..048d875d 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -5,78 +5,46 @@ using TensorKit using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence +spacetypes = [ComplexSpace, Z2Space] scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] -maxiter = 400 ctmrg_algs = [SequentialCTMRG, SimultaneousCTMRG] projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] -χ = 8 +χ = 6 atol = 1e-4 -function _make_symmetric!(psi) - if ==(size(psi)...) - return symmetrize!(psi, RotateReflect()) - else - return symmetrize!(symmetrize!(psi, ReflectDepth()), ReflectWidth()) - end +function _pre_converge_env( + ::Type{T}, physical_space, peps_space, ctm_space, unitcell; seed=12345 +) where {T} + Random.seed!(seed) # Seed RNG to make random environment consistent + psi = InfinitePEPS(rand, T, physical_space, peps_space; unitcell) + env₀ = CTMRGEnv(psi, ctm_space) + env_conv = leading_boundary(env₀, psi, SequentialCTMRG()) + return env_conv, psi end -# If I can't make the rng seed behave, I'll just randomly define a peps somehow -function _semi_random_peps!(psi::InfinitePEPS) - i = 0 - A′ = map(psi.A) do a - for (_, b) in blocks(a) - l = length(b) - b .= reshape(collect((1:l) .+ i), size(b)) - i += l - end - return a +# pre-converge CTMRG environments with given spacetype, scalartype and unit cell +preconv = Dict() +for (S, T, unitcell) in Iterators.product(spacetypes, scalartypes, unitcells) + if S == ComplexSpace + result = _pre_converge_env(T, S(2), S(2), S(χ), unitcell) + elseif S == Z2Space + result = _pre_converge_env( + T, S(0 => 1, 1 => 1), S(0 => 1, 1 => 1), S(0 => χ ÷ 2, 1 => χ ÷ 2), unitcell + ) end - return InfinitePEPS(A′) -end - -@testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( - T, unitcell, ctmrg_alg, projector_alg -) in Iterators.product( - scalartypes, unitcells, ctmrg_algs, projector_algs -) - physical_space = ComplexSpace(2) - peps_space = ComplexSpace(2) - ctm_space = ComplexSpace(χ) - - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) - - Random.seed!(987654321) # Seed RNG to make random environment consistent - env = CTMRGEnv(psi, ctm_space) - alg = ctmrg_alg(; maxiter, projector_alg) - - env = leading_boundary(env, psi, alg) - env′, = ctmrg_iteration(psi, env, alg) - env_fixed, = gauge_fix(env, env′) - @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol + push!(preconv, (S, T, unitcell) => result) end -@testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( - T, unitcell, ctmrg_alg, projector_alg +@testset "($S) - ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( + S, T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( - scalartypes, unitcells, ctmrg_algs, projector_algs + spacetypes, scalartypes, unitcells, ctmrg_algs, projector_algs ) - physical_space = Z2Space(0 => 1, 1 => 1) - peps_space = Z2Space(0 => 1, 1 => 1) - ctm_space = Z2Space(0 => χ ÷ 2, 1 => χ ÷ 2) - - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) - - Random.seed!(29385293852) # Seed RNG to make random environment consistent - psi = InfinitePEPS(physical_space, peps_space; unitcell) - env = CTMRGEnv(psi, ctm_space) - alg = ctmrg_alg(; maxiter, projector_alg) - - env = leading_boundary(env, psi, alg) + alg = ctmrg_alg(; projector_alg) + env_pre, psi = preconv[(S, T, unitcell)] + env_pre + env = leading_boundary(env_pre, psi, alg) env′, = ctmrg_iteration(psi, env, alg) env_fixed, = gauge_fix(env, env′) @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol From 528da01d7954715401e3d51d5805ba5877085870 Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Tue, 10 Dec 2024 15:51:11 +0100 Subject: [PATCH 18/19] Fix unitcell test --- test/ctmrg/unitcell.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index e2f010fa..126126f2 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -53,7 +53,7 @@ function random_dualize!(M::AbstractMatrix{<:ElementarySpace}) return M end -@testset "Integer space specifiers with $ctm_alg" begin +@testset "Integer space specifiers with $ctm_alg" for ctm_alg in ctm_algs unitcell = (3, 3) Pspaces = rand(2:3, unitcell...) From 1c4800c7035ee04590e6013012f7043663b66d1e Mon Sep 17 00:00:00 2001 From: Paul Brehmer Date: Thu, 12 Dec 2024 15:32:33 +0100 Subject: [PATCH 19/19] Fix formatting --- src/environments/ctmrg_environments.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index e40f2c0e..97d890de 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -418,7 +418,6 @@ function Base.rot180(env::CTMRGEnv{C,T}) where {C,T} return CTMRGEnv(copy(corners′), copy(edges′)) end - # Functions used for FP differentiation and by KrylovKit.linsolve function Base.:+(e₁::CTMRGEnv, e₂::CTMRGEnv) return CTMRGEnv(e₁.corners + e₂.corners, e₁.edges + e₂.edges)