From 5204349ed999483403452971c73373a859966ffa Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 17:27:00 +0800 Subject: [PATCH 01/32] Add local approximation of the product of two iPEPO --- src/PEPSKit.jl | 3 + src/algorithms/approximate/approx_tools.jl | 60 ++++++++++ src/algorithms/approximate/local_approx.jl | 125 +++++++++++++++++++++ 3 files changed, 188 insertions(+) create mode 100644 src/algorithms/approximate/approx_tools.jl create mode 100644 src/algorithms/approximate/local_approx.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4a39d6bd3..05c2aec5c 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -73,6 +73,9 @@ include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_truncation.jl") +include("algorithms/approximate/approx_tools.jl") +include("algorithms/approximate/local_approx.jl") + include("algorithms/time_evolution/evoltools.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl new file mode 100644 index 000000000..06c61a1c9 --- /dev/null +++ b/src/algorithms/approximate/approx_tools.jl @@ -0,0 +1,60 @@ +""" +$(TYPEDEF) + +Abstract super type for the truncation algorithms of +virtual spaces in InfinitePEPO or InfinitePEPS. +""" +abstract type ApproxAlgorithm end + +function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) + Nr, Nc = size(state) + if isa(state, InfinitePEPO) + @assert size(state, 3) == 1 + end + flip_xs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + return isdual(virtualspace(state[r, c], EAST)) + end + flip_ys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + return isdual(virtualspace(state[r, c], NORTH)) + end + return flip_xs, flip_ys +end + +function _standardize_virtual_spaces!( + state::Union{InfinitePEPS, InfinitePEPO}, + flip_xs::AbstractMatrix{Bool}, flip_ys::AbstractMatrix{Bool}; + inv::Bool = false + ) + Nr, Nc = size(flip_xs) + for r in 1:Nr, c in 1:Nc + inds = [ + flip_ys[r, c], flip_xs[r, c], + flip_ys[_next(r, Nr), c], flip_xs[r, _prev(c, Nc)] + ] + state.A[r, c] = flip(state.A[r, c], inds; inv) + end + return state +end + +""" +Flip all north and east virtual spaces of `state` to non-dual spaces. +A new state is constructed. +""" +function standardize_virtual_spaces( + state::Union{InfinitePEPS, InfinitePEPO}; inv::Bool = false + ) + flip_xs, flip_ys = _check_virtual_dualness(state) + !all(flip_xs) && !all(flip_ys) && (return state) + return _standardize_virtual_spaces!(deepcopy(state), flip_xs, flip_ys; inv) +end +""" +Flip all north and east virtual spaces of `state` to non-dual spaces. +Changes are in place. +""" +function standardize_virtual_spaces!( + state::Union{InfinitePEPS, InfinitePEPO}; inv::Bool = false + ) + flip_xs, flip_ys = _check_virtual_dualness(state) + !all(flip_xs) && !all(flip_ys) && (return state) + return _standardize_virtual_spaces!(state, flip_xs, flip_ys; inv) +end diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl new file mode 100644 index 000000000..f99625e6e --- /dev/null +++ b/src/algorithms/approximate/local_approx.jl @@ -0,0 +1,125 @@ +struct LocalApprox <: ApproxAlgorithm + trunc::TruncationStrategy +end + +""" +Calculate the QR decomposition of 2-layer PEPO tensor +``` + ↓ ╱ + ----A2-←- ┌-←- + ╱ | | + ↓ = (Q)-←--R + | ╱ | + ----A1-←- └-←- + ╱ ↓ +``` +Only `R` is calculated and returned. +""" +function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) + @assert isdual(space(A1, 4)) && isdual(space(A2, 4)) + A2′ = twistdual(A2, [2, 3, 5, 6]) + A1′ = twistdual(A1, [1, 3, 5, 6]) + @tensoropt MdagM[x2 z z′ x2′] := + conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] + @tensoropt MdagM[x1 x2; x1′ x2′] := MdagM[x2 z z′ x2′] * + conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] + # TODO: switch to eigh + _, s, R = svd_compact!(MdagM) + R = sdiag_pow(s, 0.5) * R + return R +end + +""" +Calculate the LQ decomposition of 2-layer PEPO tensor +``` + ↓ ╱ + --←-A2--- -←-┐ + ╱ | | + ↓ = L--←-(Q) + | ╱ | + --←-A1--- -←-┘ + ╱ ↓ +``` +Only `L` is calculated and returned. +""" +function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) + @assert !isdual(space(A1, 6)) && !isdual(space(A2, 6)) + A2′ = twistdual(A2, [2, 3, 4, 5]) + A1′ = twistdual(A1, [1, 3, 4, 5]) + @tensoropt MMdag[x2 z z′ x2′] := + A2[z z2; Y2 X2 y2 x2] * conj(A2′[z′ z2; Y2 X2 y2 x2′]) + @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * + A1[z1 z; Y1 X1 y1 x1] * conj(A1′[z1 z′; Y1 X1 y1 x1′]) + # TODO: switch to eigh + L, s, _ = svd_compact!(MMdag) + L = L * sdiag_pow(s, 0.5) + return L +end + +""" +Find the local projector `P1`, `P2` for the +following truncation of two layers of InfinitePEPO +``` + ↓ ╱ ↓ ╱ + ----A2-←-|╲ ╱|--←-B2--- + ╱ | | ╲ ╱ | ╱ | + ↓ |P1├-←-┤P2| ↓ + | ╱ | ╱ ╲ | | ╱ + ----A1-←-|╱ ╲|--←-B1--- + ╱ ↓ ╱ ↓ +``` +Reference: Physical Review B 100, 035449 (2019) +""" +function localapprox_projector( + A1::PEPOTensor, A2::PEPOTensor, B1::PEPOTensor, B2::PEPOTensor; + trunc::TruncationStrategy + ) + R1 = qr_twolayer(A1, A2) + R2 = lq_twolayer(B1, B2) + u, s, vh = svd_compact!(R1 * R2) + u, s, vh, ϵ = _truncate_compact((u, s, vh), trunc) + sinv_sqrt = sdiag_pow(s, -0.5) + P1 = R2 * vh' * sinv_sqrt + P2 = sinv_sqrt * u' * R1 + return P1, s, P2, ϵ +end + +""" +Compute an approximation to the product of two 1-layer InfinitePEPOs `ρ1`, `ρ2` +with virtual bond truncated with `LocalApprox`. +""" +function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalApprox) + @assert size(ρ1) == size(ρ2) + @assert size(ρ1, 3) == size(ρ2, 3) == 1 + Nr, Nc, = size(ρ1) + ρ1 = standardize_virtual_spaces(ρ1) + ρ2 = standardize_virtual_spaces(ρ2) + # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] + Pxs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + P1, s, P2, ϵ = localapprox_projector( + ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; + trunc = alg.trunc + ) + return (P1, P2) + end + # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] + Pys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + # TODO: reduce repeated rotations + P1, s, P2, ϵ = localapprox_projector( + rotr90(ρ1[r, c]), rotr90(ρ2[r, c]), + rotr90(ρ1[_prev(r, Nr), c]), rotr90(ρ2[_prev(r, Nr), c]); + trunc = alg.trunc + ) + return (P1, P2) + end + # apply projectors + As = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + Pw, Pe = Pxs[r, _prev(c, Nc)][2], Pxs[r, c][1] + Pn, Ps = Pys[r, c][1], Pys[_next(r, Nr), c][2] + @tensoropt A[p1 p2; n e s w] := + (ρ1[r, c])[p1 p; n1 e1 s1 w1] * (ρ2[r, c])[p p2; n2 e2 s2 w2] * + Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] + return A + end + return InfinitePEPO(As) +end From dc81186b63c84cf066b51352b4a8cdc5a4996298 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 17:27:09 +0800 Subject: [PATCH 02/32] Test LocalApprox projectors --- test/approximate/local_approx.jl | 47 ++++++++++++++++++++++++++++++++ test/runtests.jl | 5 ++++ 2 files changed, 52 insertions(+) create mode 100644 test/approximate/local_approx.jl diff --git a/test/approximate/local_approx.jl b/test/approximate/local_approx.jl new file mode 100644 index 000000000..4052871ed --- /dev/null +++ b/test/approximate/local_approx.jl @@ -0,0 +1,47 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using PEPSKit +using PEPSKit: localapprox_projector + +""" +Cost function of LocalApprox +``` + ↓ ╱ ↓ ╱ ↓ ╱ ↓ ╱ + ----A2---←---B2--- ----A2-←-|╲ ╱|--←-B2--- + ╱ | ╱ | ╱ | | ╲ ╱ | ╱ | + ↓ ↓ - ↓ |P1├-←-┤P2| ↓ + | ╱ | ╱ | ╱ | ╱ ╲ | | ╱ + ----A1---←---B1--- ----A1-←-|╱ ╲|--←-B1--- + ╱ ↓ ╱ ↓ ╱ ↓ ╱ ↓ +``` +For test convenience, open virtual indices are made trivial and removed. +""" +function localapprox_cost(A1, A2, B1, B2, P1, P2) + @tensor net1[pa1 pb1; pa2′ pb2′] := + A1[pa1 pa; D1] * A2[pa pa2′; D2] * B1[pb1 pb; D1] * B2[pb pb2′; D2] + @tensor net2[pa1 pb1; pa2′ pb2′] := P1[Da1 Da2; D] * P2[D; Db1 Db2] * + A1[pa1 pa; Da1] * A2[pa pa2′; Da2] * B1[pb1 pb; Db1] * B2[pb pb2′; Db2] + return norm(net1 - net2) +end + +@testset "Cost function of LocalApprox" begin + Random.seed!(0) + Vaux, Vphy, V = ℂ^1, ℂ^10, ℂ^4 + A1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ V ⊗ Vaux' ⊗ Vaux'), Inf) + A2 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ V ⊗ Vaux' ⊗ Vaux'), Inf) + B1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) + B2 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) + + P1, s, P2, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = notrunc()) + @test P1 * P2 ≈ TensorKit.id(domain(P2)) + + P1, sc, P2, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = truncrank(8)) + A1 = removeunit(removeunit(removeunit(A1, 6), 5), 3) + A2 = removeunit(removeunit(removeunit(A2, 6), 5), 3) + B1 = removeunit(removeunit(removeunit(B1, 5), 4), 3) + B2 = removeunit(removeunit(removeunit(B2, 5), 4), 3) + @info "Truncation error = $(ϵ)." + @test ϵ ≈ localapprox_cost(A1, A2, B1, B2, P1, P2) +end diff --git a/test/runtests.jl b/test/runtests.jl index 0f7a748c7..a58a97a15 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -85,6 +85,11 @@ end include("timeevol/sitedep_truncation.jl") end end + if GROUP == "ALL" || GROUP == "APPROX" + @time @safetestset "Local approximation" begin + include("approximate/local_approx.jl") + end + end if GROUP == "ALL" || GROUP == "TOOLBOX" @time @safetestset "Density matrix from double-layer PEPO" begin include("toolbox/densitymatrices.jl") From 77e4e50e99ed374918da5eb55c2f1aa911cdac3b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 17:42:56 +0800 Subject: [PATCH 03/32] Test `approximate` with time evolution --- src/PEPSKit.jl | 3 +++ src/algorithms/approximate/local_approx.jl | 2 +- test/examples/tf_ising_finiteT.jl | 8 ++++++++ 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 05c2aec5c..5ef6bbf32 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -106,6 +106,9 @@ export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint +export LocalApprox +export approximate + export absorb_weight export ALSTruncation, FullEnvTruncation export SimpleUpdate diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index f99625e6e..e60e126de 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -121,5 +121,5 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end - return InfinitePEPO(As) + return InfinitePEPO(cat(As; dims = 3)) end diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl index 7cefb841d..3686f1cfb 100644 --- a/test/examples/tf_ising_finiteT.jl +++ b/test/examples/tf_ising_finiteT.jl @@ -65,6 +65,14 @@ dt, nstep = 1.0e-3, 400 @test β ≈ info.t @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) + # use `approximate` to reach 2β + pepo2 = approximate(pepo, pepo, LocalApprox(trunc_pepo)) + normalize!.(pepo2.A) + env2 = converge_env(InfinitePartitionFunction(pepo2), 16) + result_2β = measure_mag(pepo2, env2) + @info "tr(σ(x,z)ρ) at T = $(1 / (2β)): $(result_2β)." + @test isapprox(abs.(result_2β), bm_2β, rtol = 5.0e-3) + # continue to get results at 2β, or T = 1.25 pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0 = β) env = converge_env(InfinitePartitionFunction(pepo), 16) From a740dbf0cb3623aca6ec6be14e532ef9f43079e0 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 19:26:08 +0800 Subject: [PATCH 04/32] Implement `standardize_virtual_spaces` --- src/algorithms/approximate/approx_tools.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl index 06c61a1c9..099717681 100644 --- a/src/algorithms/approximate/approx_tools.jl +++ b/src/algorithms/approximate/approx_tools.jl @@ -12,10 +12,10 @@ function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) @assert size(state, 3) == 1 end flip_xs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - return isdual(virtualspace(state[r, c], EAST)) + return !isdual(virtualspace(state[r, c], EAST)) end flip_ys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - return isdual(virtualspace(state[r, c], NORTH)) + return !isdual(virtualspace(state[r, c], NORTH)) end return flip_xs, flip_ys end @@ -27,10 +27,12 @@ function _standardize_virtual_spaces!( ) Nr, Nc = size(flip_xs) for r in 1:Nr, c in 1:Nc - inds = [ - flip_ys[r, c], flip_xs[r, c], - flip_ys[_next(r, Nr), c], flip_xs[r, _prev(c, Nc)] - ] + inds = findall( + [ + flip_ys[r, c], flip_xs[r, c], + flip_ys[_next(r, Nr), c], flip_xs[r, _prev(c, Nc)], + ] + ) .+ (isa(state, InfinitePEPS) ? 1 : 2) state.A[r, c] = flip(state.A[r, c], inds; inv) end return state From 04a032028dd434e4339aa1ddb7c8c93ee63ca22a Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 19:31:11 +0800 Subject: [PATCH 05/32] Some renaming --- src/algorithms/approximate/approx_tools.jl | 8 ++++---- src/algorithms/approximate/local_approx.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl index 099717681..f34ebc199 100644 --- a/src/algorithms/approximate/approx_tools.jl +++ b/src/algorithms/approximate/approx_tools.jl @@ -4,7 +4,7 @@ $(TYPEDEF) Abstract super type for the truncation algorithms of virtual spaces in InfinitePEPO or InfinitePEPS. """ -abstract type ApproxAlgorithm end +abstract type ApproximateAlgorithm end function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) Nr, Nc = size(state) @@ -20,7 +20,7 @@ function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) return flip_xs, flip_ys end -function _standardize_virtual_spaces!( +function _flip_virtual_spaces!( state::Union{InfinitePEPS, InfinitePEPO}, flip_xs::AbstractMatrix{Bool}, flip_ys::AbstractMatrix{Bool}; inv::Bool = false @@ -47,7 +47,7 @@ function standardize_virtual_spaces( ) flip_xs, flip_ys = _check_virtual_dualness(state) !all(flip_xs) && !all(flip_ys) && (return state) - return _standardize_virtual_spaces!(deepcopy(state), flip_xs, flip_ys; inv) + return _flip_virtual_spaces!(deepcopy(state), flip_xs, flip_ys; inv) end """ Flip all north and east virtual spaces of `state` to non-dual spaces. @@ -58,5 +58,5 @@ function standardize_virtual_spaces!( ) flip_xs, flip_ys = _check_virtual_dualness(state) !all(flip_xs) && !all(flip_ys) && (return state) - return _standardize_virtual_spaces!(state, flip_xs, flip_ys; inv) + return _flip_virtual_spaces!(state, flip_xs, flip_ys; inv) end diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index e60e126de..6af200401 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -1,4 +1,4 @@ -struct LocalApprox <: ApproxAlgorithm +struct LocalApprox <: ApproximateAlgorithm trunc::TruncationStrategy end From 5acf0f99036491823fce237fba9cabcf2476cce4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 19:52:29 +0800 Subject: [PATCH 06/32] Add `approx` to test groups --- .github/workflows/Tests.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 621e19c4b..7049aeff3 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -34,6 +34,7 @@ jobs: - utility - bondenv - timeevol + - approx - toolbox os: - ubuntu-latest From 90b9286a7d239d0fa4c315d9f36bb7ea1c8d3554 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 13 Mar 2026 13:46:13 +0800 Subject: [PATCH 07/32] Use svd_trunc! to directly get truncation error --- src/algorithms/approximate/local_approx.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 6af200401..3e850b840 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -76,8 +76,7 @@ function localapprox_projector( ) R1 = qr_twolayer(A1, A2) R2 = lq_twolayer(B1, B2) - u, s, vh = svd_compact!(R1 * R2) - u, s, vh, ϵ = _truncate_compact((u, s, vh), trunc) + u, s, vh, ϵ = svd_trunc!(R1 * R2; trunc) sinv_sqrt = sdiag_pow(s, -0.5) P1 = R2 * vh' * sinv_sqrt P2 = sinv_sqrt * u' * R1 From 42c12476680a8682bb22afdf8ea22325f94c3e79 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 13 Mar 2026 21:00:54 +0800 Subject: [PATCH 08/32] Use eigh instead of svd --- src/algorithms/approximate/local_approx.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 3e850b840..574c20506 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -16,16 +16,15 @@ Calculate the QR decomposition of 2-layer PEPO tensor Only `R` is calculated and returned. """ function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) - @assert isdual(space(A1, 4)) && isdual(space(A2, 4)) + @assert isdual(virtualspace(A1, EAST)) && isdual(virtualspace(A2, EAST)) A2′ = twistdual(A2, [2, 3, 5, 6]) A1′ = twistdual(A1, [1, 3, 5, 6]) @tensoropt MdagM[x2 z z′ x2′] := conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] @tensoropt MdagM[x1 x2; x1′ x2′] := MdagM[x2 z z′ x2′] * conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] - # TODO: switch to eigh - _, s, R = svd_compact!(MdagM) - R = sdiag_pow(s, 0.5) * R + D, R = eigh_full!(MdagM) + R = sdiag_pow(D, 0.5) * R' return R end @@ -43,16 +42,15 @@ Calculate the LQ decomposition of 2-layer PEPO tensor Only `L` is calculated and returned. """ function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) - @assert !isdual(space(A1, 6)) && !isdual(space(A2, 6)) + @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) A2′ = twistdual(A2, [2, 3, 4, 5]) A1′ = twistdual(A1, [1, 3, 4, 5]) @tensoropt MMdag[x2 z z′ x2′] := A2[z z2; Y2 X2 y2 x2] * conj(A2′[z′ z2; Y2 X2 y2 x2′]) @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * A1[z1 z; Y1 X1 y1 x1] * conj(A1′[z1 z′; Y1 X1 y1 x1′]) - # TODO: switch to eigh - L, s, _ = svd_compact!(MMdag) - L = L * sdiag_pow(s, 0.5) + D, L = eigh_full!(MMdag) + L = L * sdiag_pow(D, 0.5) return L end From 98457078acc3fbac81c0a5dd29352e405fb16c6d Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 14 Mar 2026 20:55:43 +0800 Subject: [PATCH 09/32] Use eigh_trunc! to avoid small negative eigenvalues --- src/algorithms/approximate/local_approx.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 574c20506..0ef1d2739 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -23,7 +23,8 @@ function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] @tensoropt MdagM[x1 x2; x1′ x2′] := MdagM[x2 z z′ x2′] * conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] - D, R = eigh_full!(MdagM) + # avoid small negative eigenvalues due to numerical errors + D, R = eigh_trunc!(MdagM; trunc = trunctol(; atol = 1.0e-16)) R = sdiag_pow(D, 0.5) * R' return R end @@ -49,7 +50,8 @@ function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) A2[z z2; Y2 X2 y2 x2] * conj(A2′[z′ z2; Y2 X2 y2 x2′]) @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * A1[z1 z; Y1 X1 y1 x1] * conj(A1′[z1 z′; Y1 X1 y1 x1′]) - D, L = eigh_full!(MMdag) + # avoid small negative eigenvalues due to numerical errors + D, L = eigh_trunc!(MMdag; trunc = trunctol(; atol = 1.0e-16)) L = L * sdiag_pow(D, 0.5) return L end From 4e16e253bf9326b61cd8189d56beb2c752ab762b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 15 Mar 2026 08:29:38 +0800 Subject: [PATCH 10/32] Fix twists for fermions --- src/algorithms/approximate/local_approx.jl | 4 ++-- src/utility/util.jl | 16 ++++++++++++++++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 0ef1d2739..75e9b3e18 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -44,8 +44,8 @@ Only `L` is calculated and returned. """ function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) - A2′ = twistdual(A2, [2, 3, 4, 5]) - A1′ = twistdual(A1, [1, 3, 4, 5]) + A2′ = twistnondual(A2, [2, 3, 4, 5]) + A1′ = twistnondual(A1, [1, 3, 4, 5]) @tensoropt MMdag[x2 z z′ x2′] := A2[z z2; Y2 X2 y2 x2] * conj(A2′[z′ z2; Y2 X2 y2 x2′]) @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * diff --git a/src/utility/util.jl b/src/utility/util.jl index 3fa7e0cc1..d563e13f4 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -118,6 +118,22 @@ function twistdual!(t::AbstractTensorMap, is) end twistdual(t::AbstractTensorMap, is) = twistdual!(copy(t), is) +""" + twistnondual(t::AbstractTensorMap, i) + twistnondual!(t::AbstractTensorMap, i) + +Twist the i-th leg of a tensor `t` if it represents a non-dual space. +""" +function twistnondual!(t::AbstractTensorMap, i::Int) + !isdual(space(t, i)) || return t + return twist!(t, i) +end +function twistnondual!(t::AbstractTensorMap, is) + is′ = filter(i -> !isdual(space(t, i)), is) + return twist!(t, is′) +end +twistnondual(t::AbstractTensorMap, is) = twistnondual!(copy(t), is) + """ str(t) From 1860b2e40806722943a4e3fbe13b75f82a26c190 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 17 Mar 2026 09:34:57 +0800 Subject: [PATCH 11/32] Return truncation error --- src/algorithms/approximate/local_approx.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 75e9b3e18..c708dbb87 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -99,7 +99,7 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; trunc = alg.trunc ) - return (P1, P2) + return (P1, P2, ϵ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] Pys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) @@ -109,7 +109,7 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr rotr90(ρ1[_prev(r, Nr), c]), rotr90(ρ2[_prev(r, Nr), c]); trunc = alg.trunc ) - return (P1, P2) + return (P1, P2, ϵ) end # apply projectors As = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) @@ -120,5 +120,8 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end - return InfinitePEPO(cat(As; dims = 3)) + ϵx = map(Base.Fix2(getindex, 3), Pxs) + ϵy = map(Base.Fix2(getindex, 3), Pys) + info = (; ϵ = max(ϵx, ϵy)) + return InfinitePEPO(cat(As; dims = 3)), info end From f18a590a316c58c945a8526cdbec51e986ba180c Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 17 Mar 2026 11:37:47 +0800 Subject: [PATCH 12/32] Fix tests --- src/algorithms/approximate/local_approx.jl | 4 ++-- test/timeevol/tf_ising_finiteT.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index c708dbb87..88bc99e5c 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -120,8 +120,8 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end - ϵx = map(Base.Fix2(getindex, 3), Pxs) - ϵy = map(Base.Fix2(getindex, 3), Pys) + ϵx = maximum(map(Base.Fix2(getindex, 3), Pxs)) + ϵy = maximum(map(Base.Fix2(getindex, 3), Pys)) info = (; ϵ = max(ϵx, ϵy)) return InfinitePEPO(cat(As; dims = 3)), info end diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index c3726f231..f1c8b75d4 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -64,7 +64,7 @@ dt, nstep = 1.0e-3, 400 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) # use `approximate` to reach 2β - pepo2 = approximate(pepo, pepo, LocalApprox(trunc_pepo)) + pepo2, = approximate(pepo, pepo, LocalApprox(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) From 7437803908bf349e66d829dea00dc6235b81b25b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 31 Mar 2026 21:12:02 +0800 Subject: [PATCH 13/32] Return projectors in info --- src/algorithms/approximate/local_approx.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 88bc99e5c..79884fa00 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -122,6 +122,6 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr end ϵx = maximum(map(Base.Fix2(getindex, 3), Pxs)) ϵy = maximum(map(Base.Fix2(getindex, 3), Pys)) - info = (; ϵ = max(ϵx, ϵy)) + info = (; Pxs, Pys, ϵ = max(ϵx, ϵy)) return InfinitePEPO(cat(As; dims = 3)), info end From 93d31f1cdfee4c088894e472d1b1d2f06654fa9e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 1 Apr 2026 17:27:36 +0800 Subject: [PATCH 14/32] Increase test coverage --- src/algorithms/approximate/local_approx.jl | 6 +++--- test/approximate/local_approx.jl | 14 ++++++++++++-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 79884fa00..2a792673e 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -80,7 +80,7 @@ function localapprox_projector( sinv_sqrt = sdiag_pow(s, -0.5) P1 = R2 * vh' * sinv_sqrt P2 = sinv_sqrt * u' * R1 - return P1, s, P2, ϵ + return P1, P2, s, ϵ end """ @@ -95,7 +95,7 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr ρ2 = standardize_virtual_spaces(ρ2) # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] Pxs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - P1, s, P2, ϵ = localapprox_projector( + P1, P2, s, ϵ = localapprox_projector( ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; trunc = alg.trunc ) @@ -104,7 +104,7 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] Pys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) # TODO: reduce repeated rotations - P1, s, P2, ϵ = localapprox_projector( + P1, P2, s, ϵ = localapprox_projector( rotr90(ρ1[r, c]), rotr90(ρ2[r, c]), rotr90(ρ1[_prev(r, Nr), c]), rotr90(ρ2[_prev(r, Nr), c]); trunc = alg.trunc diff --git a/test/approximate/local_approx.jl b/test/approximate/local_approx.jl index 4052871ed..5106e1f54 100644 --- a/test/approximate/local_approx.jl +++ b/test/approximate/local_approx.jl @@ -34,10 +34,10 @@ end B1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) B2 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) - P1, s, P2, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = notrunc()) + P1, P2, s, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = notrunc()) @test P1 * P2 ≈ TensorKit.id(domain(P2)) - P1, sc, P2, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = truncrank(8)) + P1, P2, sc, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = truncrank(8)) A1 = removeunit(removeunit(removeunit(A1, 6), 5), 3) A2 = removeunit(removeunit(removeunit(A2, 6), 5), 3) B1 = removeunit(removeunit(removeunit(B1, 5), 4), 3) @@ -45,3 +45,13 @@ end @info "Truncation error = $(ϵ)." @test ϵ ≈ localapprox_cost(A1, A2, B1, B2, P1, P2) end + +@testset "Virtual space matching" begin + Vps = ComplexSpace.([2 2; 2 2]) + Vns = ComplexSpace.([2 4; 5 3]) + Ves = ComplexSpace.([3 5; 4 2]) + ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) + alg = LocalApprox(truncrank(2)) + ρ2, = approximate(ρ, ρ, alg) + @test ρ2 isa InfinitePEPO +end From 557eae001e63d70d2988cc23d79fef5a3c15cc2f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 10 Apr 2026 11:35:51 +0800 Subject: [PATCH 15/32] Update src/algorithms/approximate/approx_tools.jl Co-authored-by: Lander Burgelman <39218680+leburgel@users.noreply.github.com> --- src/algorithms/approximate/approx_tools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl index f34ebc199..06c397099 100644 --- a/src/algorithms/approximate/approx_tools.jl +++ b/src/algorithms/approximate/approx_tools.jl @@ -33,7 +33,7 @@ function _flip_virtual_spaces!( flip_ys[_next(r, Nr), c], flip_xs[r, _prev(c, Nc)], ] ) .+ (isa(state, InfinitePEPS) ? 1 : 2) - state.A[r, c] = flip(state.A[r, c], inds; inv) + state[r, c] = flip(state[r, c], inds; inv) end return state end From 39028faa1cc787c24261cf5ef32fd86acaebbdd3 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 10 Apr 2026 22:19:01 +0800 Subject: [PATCH 16/32] Apply suggestions from @leburgel --- src/PEPSKit.jl | 2 +- src/algorithms/approximate/approx_tools.jl | 63 ++++------------------ src/algorithms/approximate/local_approx.jl | 49 +++++++++++------ test/approximate/local_approx.jl | 15 ++---- test/timeevol/tf_ising_finiteT.jl | 12 +---- 5 files changed, 51 insertions(+), 90 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index cd3654b1a..ae7769fc7 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -144,7 +144,7 @@ export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint -export LocalApprox +export LocalApproximation export approximate export absorb_weight diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl index 06c397099..2c618be17 100644 --- a/src/algorithms/approximate/approx_tools.jl +++ b/src/algorithms/approximate/approx_tools.jl @@ -1,62 +1,21 @@ """ $(TYPEDEF) -Abstract super type for the truncation algorithms of -virtual spaces in InfinitePEPO or InfinitePEPS. +Abstract super type for approximation algorithms +for two-dimensional tensor networks """ abstract type ApproximateAlgorithm end +""" +Checks if the east and the north virtual space of `state` +follow the standard convention, i.e. are dual spaces. +""" function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) - Nr, Nc = size(state) - if isa(state, InfinitePEPO) - @assert size(state, 3) == 1 + isdual_easts = map(CartesianIndices(state.A)) do idx + return isdual(virtualspace(state[idx], EAST)) end - flip_xs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - return !isdual(virtualspace(state[r, c], EAST)) + isdual_norths = map(CartesianIndices(state.A)) do idx + return isdual(virtualspace(state[idx], NORTH)) end - flip_ys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - return !isdual(virtualspace(state[r, c], NORTH)) - end - return flip_xs, flip_ys -end - -function _flip_virtual_spaces!( - state::Union{InfinitePEPS, InfinitePEPO}, - flip_xs::AbstractMatrix{Bool}, flip_ys::AbstractMatrix{Bool}; - inv::Bool = false - ) - Nr, Nc = size(flip_xs) - for r in 1:Nr, c in 1:Nc - inds = findall( - [ - flip_ys[r, c], flip_xs[r, c], - flip_ys[_next(r, Nr), c], flip_xs[r, _prev(c, Nc)], - ] - ) .+ (isa(state, InfinitePEPS) ? 1 : 2) - state[r, c] = flip(state[r, c], inds; inv) - end - return state -end - -""" -Flip all north and east virtual spaces of `state` to non-dual spaces. -A new state is constructed. -""" -function standardize_virtual_spaces( - state::Union{InfinitePEPS, InfinitePEPO}; inv::Bool = false - ) - flip_xs, flip_ys = _check_virtual_dualness(state) - !all(flip_xs) && !all(flip_ys) && (return state) - return _flip_virtual_spaces!(deepcopy(state), flip_xs, flip_ys; inv) -end -""" -Flip all north and east virtual spaces of `state` to non-dual spaces. -Changes are in place. -""" -function standardize_virtual_spaces!( - state::Union{InfinitePEPS, InfinitePEPO}; inv::Bool = false - ) - flip_xs, flip_ys = _check_virtual_dualness(state) - !all(flip_xs) && !all(flip_ys) && (return state) - return _flip_virtual_spaces!(state, flip_xs, flip_ys; inv) + return isdual_easts, isdual_norths end diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index 2a792673e..ecfc4104c 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -1,9 +1,27 @@ -struct LocalApprox <: ApproximateAlgorithm +""" +$(TYPEDEF) + +Algorithm to approximate a two-layer network by truncating each pair +of virtual spaces with projectors that minimizes the local cost function, +which is the 2-norm of (e.g. for two layers of iPEPO) +``` + ↓ ╱ ↓ ╱ ↓ ╱ ↓ ╱ + ----A2---←---B2--- ----A2-←-|╲ ╱|--←-B2--- + ╱ | ╱ | ╱ | | ╲ ╱ | ╱ | + ↓ ↓ - ↓ |P1├-←-┤P2| ↓ + | ╱ | ╱ | ╱ | ╱ ╲ | | ╱ + ----A1---←---B1--- ----A1-←-|╱ ╲|--←-B1--- + ╱ ↓ ╱ ↓ ╱ ↓ ╱ ↓ +``` +on each bond of the network. +""" +struct LocalApproximation <: ApproximateAlgorithm trunc::TruncationStrategy end """ Calculate the QR decomposition of 2-layer PEPO tensor +with the east virtual legs transferred to the R tensor ``` ↓ ╱ ----A2-←- ┌-←- @@ -31,6 +49,7 @@ end """ Calculate the LQ decomposition of 2-layer PEPO tensor +with the west virtual legs transferred to the L tensor ``` ↓ ╱ --←-A2--- -←-┐ @@ -85,16 +104,18 @@ end """ Compute an approximation to the product of two 1-layer InfinitePEPOs `ρ1`, `ρ2` -with virtual bond truncated with `LocalApprox`. +with virtual bond truncated with `LocalApproximation`. """ -function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalApprox) - @assert size(ρ1) == size(ρ2) - @assert size(ρ1, 3) == size(ρ2, 3) == 1 - Nr, Nc, = size(ρ1) - ρ1 = standardize_virtual_spaces(ρ1) - ρ2 = standardize_virtual_spaces(ρ2) +function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalApproximation) + # sanity checks + (size(ρ1) == size(ρ2)) || error("Input PEPOs have different unit cell sizes.") + (size(ρ1, 3) == 1) || error("ρ1 should have only one layer.") + (size(ρ2, 3) == 1) || error("ρ2 should have only one layer.") + all(all.(_check_virtual_dualness(ρ1))) || error("East and north virtual spaces in ρ1 should be dual spaces.") + all(all.(_check_virtual_dualness(ρ2))) || error("East and north virtual spaces in ρ2 should be dual spaces.") # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] - Pxs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + Nr, Nc, = size(ρ1) + Px_errs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) P1, P2, s, ϵ = localapprox_projector( ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; trunc = alg.trunc @@ -102,7 +123,7 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr return (P1, P2, ϵ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] - Pys = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + Py_errs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) # TODO: reduce repeated rotations P1, P2, s, ϵ = localapprox_projector( rotr90(ρ1[r, c]), rotr90(ρ2[r, c]), @@ -113,15 +134,13 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr end # apply projectors As = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - Pw, Pe = Pxs[r, _prev(c, Nc)][2], Pxs[r, c][1] - Pn, Ps = Pys[r, c][1], Pys[_next(r, Nr), c][2] + Pw, Pe = Px_errs[r, _prev(c, Nc)][2], Px_errs[r, c][1] + Pn, Ps = Py_errs[r, c][1], Py_errs[_next(r, Nr), c][2] @tensoropt A[p1 p2; n e s w] := (ρ1[r, c])[p1 p; n1 e1 s1 w1] * (ρ2[r, c])[p p2; n2 e2 s2 w2] * Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end - ϵx = maximum(map(Base.Fix2(getindex, 3), Pxs)) - ϵy = maximum(map(Base.Fix2(getindex, 3), Pys)) - info = (; Pxs, Pys, ϵ = max(ϵx, ϵy)) + info = (; Px_errs, Py_errs) return InfinitePEPO(cat(As; dims = 3)), info end diff --git a/test/approximate/local_approx.jl b/test/approximate/local_approx.jl index 5106e1f54..4f6bc51ad 100644 --- a/test/approximate/local_approx.jl +++ b/test/approximate/local_approx.jl @@ -6,16 +6,7 @@ using PEPSKit using PEPSKit: localapprox_projector """ -Cost function of LocalApprox -``` - ↓ ╱ ↓ ╱ ↓ ╱ ↓ ╱ - ----A2---←---B2--- ----A2-←-|╲ ╱|--←-B2--- - ╱ | ╱ | ╱ | | ╲ ╱ | ╱ | - ↓ ↓ - ↓ |P1├-←-┤P2| ↓ - | ╱ | ╱ | ╱ | ╱ ╲ | | ╱ - ----A1---←---B1--- ----A1-←-|╱ ╲|--←-B1--- - ╱ ↓ ╱ ↓ ╱ ↓ ╱ ↓ -``` +Cost function of LocalApproximation. For test convenience, open virtual indices are made trivial and removed. """ function localapprox_cost(A1, A2, B1, B2, P1, P2) @@ -26,7 +17,7 @@ function localapprox_cost(A1, A2, B1, B2, P1, P2) return norm(net1 - net2) end -@testset "Cost function of LocalApprox" begin +@testset "Cost function of LocalApproximation" begin Random.seed!(0) Vaux, Vphy, V = ℂ^1, ℂ^10, ℂ^4 A1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ V ⊗ Vaux' ⊗ Vaux'), Inf) @@ -51,7 +42,7 @@ end Vns = ComplexSpace.([2 4; 5 3]) Ves = ComplexSpace.([3 5; 4 2]) ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) - alg = LocalApprox(truncrank(2)) + alg = LocalApproximation(truncrank(2)) ρ2, = approximate(ρ, ρ, alg) @test ρ2 isa InfinitePEPO end diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index 3409b98d4..de57dcdee 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -64,22 +64,14 @@ dt, nstep = 1.0e-3, 400 @test β ≈ info.t @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) - # use `approximate` to reach 2β - pepo2, = approximate(pepo, pepo, LocalApprox(trunc_pepo)) + # use `approximate` to reach 2β, or T = 1.25 + pepo2, = approximate(pepo, pepo, LocalApproximation(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) @info "tr(σ(x,z)ρ) at T = $(1 / (2β)): $(result_2β)." @test isapprox(abs.(result_2β), bm_2β, rtol = 5.0e-3) - # continue to get results at 2β, or T = 1.25 - pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0 = β, symmetrize_gates) - env = converge_env(InfinitePartitionFunction(pepo), 16) - result_2β = measure_mag(pepo, env) - @info "tr(σ(x,z)ρ) at T = $(1 / (2β)): $(result_2β)." - @test 2 * β ≈ info.t - @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4) - # Purification approach: results at 2β, or T = 1.25 alg = SimpleUpdate(; trunc = trunc_pepo, purified = true, bipartite, force_mpo) pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0; symmetrize_gates) From 6b53718384fa500157e6aa49ec619dd7ee7ba800 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 11 Apr 2026 06:55:32 +0800 Subject: [PATCH 17/32] Add hermitian projection in qr/lq_twolayer; test twists --- src/algorithms/approximate/local_approx.jl | 30 ++++++++++++++-------- test/approximate/local_approx.jl | 13 ++++++++++ 2 files changed, 33 insertions(+), 10 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index ecfc4104c..f32c31902 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -34,6 +34,13 @@ with the east virtual legs transferred to the R tensor Only `R` is calculated and returned. """ function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) + MdagM = _get_MdagM(A1, A2) + D, R = eigh_full!(MdagM) + # TODO: set small negative eigenvalues due to numerical errors to 0 + R = sdiag_pow(D, 0.5) * R' + return R +end +function _get_MdagM(A1::PEPOTensor, A2::PEPOTensor) @assert isdual(virtualspace(A1, EAST)) && isdual(virtualspace(A2, EAST)) A2′ = twistdual(A2, [2, 3, 5, 6]) A1′ = twistdual(A1, [1, 3, 5, 6]) @@ -41,10 +48,8 @@ function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] @tensoropt MdagM[x1 x2; x1′ x2′] := MdagM[x2 z z′ x2′] * conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] - # avoid small negative eigenvalues due to numerical errors - D, R = eigh_trunc!(MdagM; trunc = trunctol(; atol = 1.0e-16)) - R = sdiag_pow(D, 0.5) * R' - return R + project_hermitian!(MdagM) + return MdagM end """ @@ -62,17 +67,22 @@ with the west virtual legs transferred to the L tensor Only `L` is calculated and returned. """ function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) + MMdag = _get_MMdag(A1, A2) + D, L = eigh_full!(MMdag) + # TODO: set small negative eigenvalues due to numerical errors to 0 + L = L * sdiag_pow(D, 0.5) + return L +end +function _get_MMdag(A1::PEPOTensor, A2::PEPOTensor) @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) A2′ = twistnondual(A2, [2, 3, 4, 5]) A1′ = twistnondual(A1, [1, 3, 4, 5]) @tensoropt MMdag[x2 z z′ x2′] := - A2[z z2; Y2 X2 y2 x2] * conj(A2′[z′ z2; Y2 X2 y2 x2′]) + A2′[z z2; Y2 X2 y2 x2] * conj(A2[z′ z2; Y2 X2 y2 x2′]) @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * - A1[z1 z; Y1 X1 y1 x1] * conj(A1′[z1 z′; Y1 X1 y1 x1′]) - # avoid small negative eigenvalues due to numerical errors - D, L = eigh_trunc!(MMdag; trunc = trunctol(; atol = 1.0e-16)) - L = L * sdiag_pow(D, 0.5) - return L + A1′[z1 z; Y1 X1 y1 x1] * conj(A1[z1 z′; Y1 X1 y1 x1′]) + project_hermitian!(MMdag) + return MMdag end """ diff --git a/test/approximate/local_approx.jl b/test/approximate/local_approx.jl index 4f6bc51ad..b41b72fed 100644 --- a/test/approximate/local_approx.jl +++ b/test/approximate/local_approx.jl @@ -37,6 +37,19 @@ end @test ϵ ≈ localapprox_cost(A1, A2, B1, B2, P1, P2) end +@testset "Fermionic twists for qr/lq_twolayer" begin + Vphy = Vect[FermionParity](0 => 2, 1 => 2) + Vvir = Vect[FermionParity](0 => 2, 1 => 2) + for _ in 1:4 # multiple trials without setting seed + Aspace = (Vphy ⊗ Vphy' ← Vvir ⊗ Vvir ⊗ Vvir' ⊗ Vvir') + A1 = randn(ComplexF64, Aspace) + A2 = randn(ComplexF64, Aspace) + for MM in [PEPSKit._get_MMdag(A1, A2), PEPSKit._get_MdagM(A1, A2)] + @test isposdef(MM) + end + end +end + @testset "Virtual space matching" begin Vps = ComplexSpace.([2 2; 2 2]) Vns = ComplexSpace.([2 4; 5 3]) From ba55a23a72cadfa60adee69386e921bea0be34cb Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 11 Apr 2026 14:13:52 +0800 Subject: [PATCH 18/32] Add project_psd! to remove small negative eigenvalues --- src/algorithms/approximate/local_approx.jl | 4 ++-- src/utility/util.jl | 10 ++++++++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl index f32c31902..37140f463 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/approximate/local_approx.jl @@ -36,7 +36,7 @@ Only `R` is calculated and returned. function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) MdagM = _get_MdagM(A1, A2) D, R = eigh_full!(MdagM) - # TODO: set small negative eigenvalues due to numerical errors to 0 + project_psd!(D) R = sdiag_pow(D, 0.5) * R' return R end @@ -69,7 +69,7 @@ Only `L` is calculated and returned. function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) MMdag = _get_MMdag(A1, A2) D, L = eigh_full!(MMdag) - # TODO: set small negative eigenvalues due to numerical errors to 0 + project_psd!(D) L = L * sdiag_pow(D, 0.5) return L end diff --git a/src/utility/util.jl b/src/utility/util.jl index d563e13f4..9ac564c44 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -86,6 +86,16 @@ function ChainRulesCore.rrule( return spow, sdiag_pow_pullback end +""" +Given the eigenvalue spectrum of `D` expected from a +positive-semi-definite (PSD) map, sets small negative +eigenvalues that may arise from numerical errors to 0. +""" +function project_psd!(D::DiagonalTensorMap{T, S}) where {T, S} + D.data .= max.(zero(T), D.data) + return D +end + """ absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorMap) From 3b7f05ffd222cff5c7af9ff5a783e613410e2247 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 19 Apr 2026 18:27:14 +0800 Subject: [PATCH 19/32] Overload MPSKit.changebonds instead --- .github/workflows/Tests.yml | 2 +- src/PEPSKit.jl | 7 ++-- src/algorithms/approximate/approx_tools.jl | 21 ---------- .../local_approx.jl => changebonds/local.jl} | 41 +++++++++++-------- src/operators/infinitepepo.jl | 14 +++++++ .../local_approx.jl => changebonds/local.jl} | 20 ++++----- test/runtests.jl | 6 +-- test/timeevol/tf_ising_finiteT.jl | 4 +- 8 files changed, 57 insertions(+), 58 deletions(-) delete mode 100644 src/algorithms/approximate/approx_tools.jl rename src/algorithms/{approximate/local_approx.jl => changebonds/local.jl} (83%) rename test/{approximate/local_approx.jl => changebonds/local.jl} (78%) diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 4c0c86e89..651e42aff 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -35,7 +35,7 @@ jobs: - utility - bondenv - timeevol - - approx + - changebonds - toolbox os: - ubuntu-latest diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index ae7769fc7..588b1ff1a 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -105,8 +105,7 @@ include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_truncation.jl") -include("algorithms/approximate/approx_tools.jl") -include("algorithms/approximate/local_approx.jl") +include("algorithms/changebonds/local.jl") include("algorithms/time_evolution/trotter_gate.jl") include("algorithms/time_evolution/apply_gate.jl") @@ -144,8 +143,8 @@ export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint -export LocalApproximation -export approximate +export LocalTruncation +export changebonds export absorb_weight export ALSTruncation, FullEnvTruncation diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl deleted file mode 100644 index 2c618be17..000000000 --- a/src/algorithms/approximate/approx_tools.jl +++ /dev/null @@ -1,21 +0,0 @@ -""" -$(TYPEDEF) - -Abstract super type for approximation algorithms -for two-dimensional tensor networks -""" -abstract type ApproximateAlgorithm end - -""" -Checks if the east and the north virtual space of `state` -follow the standard convention, i.e. are dual spaces. -""" -function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) - isdual_easts = map(CartesianIndices(state.A)) do idx - return isdual(virtualspace(state[idx], EAST)) - end - isdual_norths = map(CartesianIndices(state.A)) do idx - return isdual(virtualspace(state[idx], NORTH)) - end - return isdual_easts, isdual_norths -end diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/changebonds/local.jl similarity index 83% rename from src/algorithms/approximate/local_approx.jl rename to src/algorithms/changebonds/local.jl index 37140f463..7df2f8419 100644 --- a/src/algorithms/approximate/local_approx.jl +++ b/src/algorithms/changebonds/local.jl @@ -1,8 +1,16 @@ """ $(TYPEDEF) -Algorithm to approximate a two-layer network by truncating each pair -of virtual spaces with projectors that minimizes the local cost function, +Abstract super type for algorithms to change virtual bonds +in two-dimensional tensor networks +""" +abstract type BondChangeAlgorithm end + +""" +$(TYPEDEF) + +Algorithm to truncate virtual bonds of a two-layer network +with projectors that minimizes the local cost function, which is the 2-norm of (e.g. for two layers of iPEPO) ``` ↓ ╱ ↓ ╱ ↓ ╱ ↓ ╱ @@ -15,7 +23,7 @@ which is the 2-norm of (e.g. for two layers of iPEPO) ``` on each bond of the network. """ -struct LocalApproximation <: ApproximateAlgorithm +struct LocalTruncation <: BondChangeAlgorithm trunc::TruncationStrategy end @@ -99,7 +107,7 @@ following truncation of two layers of InfinitePEPO ``` Reference: Physical Review B 100, 035449 (2019) """ -function localapprox_projector( +function virtual_projector( A1::PEPOTensor, A2::PEPOTensor, B1::PEPOTensor, B2::PEPOTensor; trunc::TruncationStrategy ) @@ -109,14 +117,14 @@ function localapprox_projector( sinv_sqrt = sdiag_pow(s, -0.5) P1 = R2 * vh' * sinv_sqrt P2 = sinv_sqrt * u' * R1 - return P1, P2, s, ϵ + return P1, P2, (; s, ϵ) end """ -Compute an approximation to the product of two 1-layer InfinitePEPOs `ρ1`, `ρ2` -with virtual bond truncated with `LocalApproximation`. +Truncate virtual bonds of the product of two 1-layer +InfinitePEPOs `ρ1`, `ρ2` with `LocalTruncation`. """ -function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalApproximation) +function MPSKit.changebonds(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) # sanity checks (size(ρ1) == size(ρ2)) || error("Input PEPOs have different unit cell sizes.") (size(ρ1, 3) == 1) || error("ρ1 should have only one layer.") @@ -125,32 +133,31 @@ function MPSKit.approximate(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalAppr all(all.(_check_virtual_dualness(ρ2))) || error("East and north virtual spaces in ρ2 should be dual spaces.") # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] Nr, Nc, = size(ρ1) - Px_errs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - P1, P2, s, ϵ = localapprox_projector( + Pxs_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + # TODO: support SiteDependentTruncation + return virtual_projector( ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; trunc = alg.trunc ) - return (P1, P2, ϵ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] - Py_errs = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) + Pys_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) # TODO: reduce repeated rotations - P1, P2, s, ϵ = localapprox_projector( + return virtual_projector( rotr90(ρ1[r, c]), rotr90(ρ2[r, c]), rotr90(ρ1[_prev(r, Nr), c]), rotr90(ρ2[_prev(r, Nr), c]); trunc = alg.trunc ) - return (P1, P2, ϵ) end # apply projectors As = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - Pw, Pe = Px_errs[r, _prev(c, Nc)][2], Px_errs[r, c][1] - Pn, Ps = Py_errs[r, c][1], Py_errs[_next(r, Nr), c][2] + Pw, Pe = Pxs_info[r, _prev(c, Nc)][2], Pxs_info[r, c][1] + Pn, Ps = Pys_info[r, c][1], Pys_info[_next(r, Nr), c][2] @tensoropt A[p1 p2; n e s w] := (ρ1[r, c])[p1 p; n1 e1 s1 w1] * (ρ2[r, c])[p p2; n2 e2 s2 w2] * Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end - info = (; Px_errs, Py_errs) + info = (; Pxs_info, Pys_info) return InfinitePEPO(cat(As; dims = 3)), info end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 172c56b05..1a65809e9 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -254,3 +254,17 @@ function _stack_tuples(A::Matrix{NTuple{N, T}}) where {N, T} end return out end + +""" +Checks if `state` follow the standard virtual arrow convention, +i.e. north and east domains are non-dual spaces. +""" +function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) + isdual_easts = map(CartesianIndices(state.A)) do idx + return isdual(virtualspace(state[idx], EAST)) + end + isdual_norths = map(CartesianIndices(state.A)) do idx + return isdual(virtualspace(state[idx], NORTH)) + end + return isdual_easts, isdual_norths +end diff --git a/test/approximate/local_approx.jl b/test/changebonds/local.jl similarity index 78% rename from test/approximate/local_approx.jl rename to test/changebonds/local.jl index b41b72fed..a1853b7fd 100644 --- a/test/approximate/local_approx.jl +++ b/test/changebonds/local.jl @@ -3,13 +3,13 @@ using Random using LinearAlgebra using TensorKit using PEPSKit -using PEPSKit: localapprox_projector +using PEPSKit: virtual_projector """ -Cost function of LocalApproximation. +Cost function of LocalTruncation. For test convenience, open virtual indices are made trivial and removed. """ -function localapprox_cost(A1, A2, B1, B2, P1, P2) +function localcompress_cost(A1, A2, B1, B2, P1, P2) @tensor net1[pa1 pb1; pa2′ pb2′] := A1[pa1 pa; D1] * A2[pa pa2′; D2] * B1[pb1 pb; D1] * B2[pb pb2′; D2] @tensor net2[pa1 pb1; pa2′ pb2′] := P1[Da1 Da2; D] * P2[D; Db1 Db2] * @@ -17,7 +17,7 @@ function localapprox_cost(A1, A2, B1, B2, P1, P2) return norm(net1 - net2) end -@testset "Cost function of LocalApproximation" begin +@testset "Cost function of LocalTruncation" begin Random.seed!(0) Vaux, Vphy, V = ℂ^1, ℂ^10, ℂ^4 A1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ V ⊗ Vaux' ⊗ Vaux'), Inf) @@ -25,16 +25,16 @@ end B1 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) B2 = normalize(randn(Vphy ⊗ Vphy' ← Vaux ⊗ Vaux ⊗ Vaux' ⊗ V'), Inf) - P1, P2, s, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = notrunc()) + P1, P2, info = virtual_projector(A1, A2, B1, B2; trunc = notrunc()) @test P1 * P2 ≈ TensorKit.id(domain(P2)) - P1, P2, sc, ϵ = localapprox_projector(A1, A2, B1, B2; trunc = truncrank(8)) + P1, P2, info = virtual_projector(A1, A2, B1, B2; trunc = truncrank(8)) A1 = removeunit(removeunit(removeunit(A1, 6), 5), 3) A2 = removeunit(removeunit(removeunit(A2, 6), 5), 3) B1 = removeunit(removeunit(removeunit(B1, 5), 4), 3) B2 = removeunit(removeunit(removeunit(B2, 5), 4), 3) - @info "Truncation error = $(ϵ)." - @test ϵ ≈ localapprox_cost(A1, A2, B1, B2, P1, P2) + @info "Truncation error = $(info.ϵ)." + @test info.ϵ ≈ localcompress_cost(A1, A2, B1, B2, P1, P2) end @testset "Fermionic twists for qr/lq_twolayer" begin @@ -55,7 +55,7 @@ end Vns = ComplexSpace.([2 4; 5 3]) Ves = ComplexSpace.([3 5; 4 2]) ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) - alg = LocalApproximation(truncrank(2)) - ρ2, = approximate(ρ, ρ, alg) + alg = LocalTruncation(truncrank(2)) + ρ2, = changebonds(ρ, ρ, alg) @test ρ2 isa InfinitePEPO end diff --git a/test/runtests.jl b/test/runtests.jl index 7a9066e68..8ff3cd426 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -114,9 +114,9 @@ end include("timeevol/j1j2_finiteT.jl") end end - if GROUP == "ALL" || GROUP == "APPROX" - @time @safetestset "Local approximation" begin - include("approximate/local_approx.jl") + if GROUP == "ALL" || GROUP == "CHANGEBONDS" + @time @safetestset "Local truncation" begin + include("changebonds/local.jl") end end if GROUP == "ALL" || GROUP == "TOOLBOX" diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index de57dcdee..52bd74c72 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -64,8 +64,8 @@ dt, nstep = 1.0e-3, 400 @test β ≈ info.t @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) - # use `approximate` to reach 2β, or T = 1.25 - pepo2, = approximate(pepo, pepo, LocalApproximation(trunc_pepo)) + # use `changebonds` to reach 2β, or T = 1.25 + pepo2, = changebonds(pepo, pepo, LocalTruncation(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) From 831a50ccc78163ae393f40d6367288e122593847 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 4 May 2026 10:11:17 +0800 Subject: [PATCH 20/32] Fix merge hiccups --- src/PEPSKit.jl | 2 -- src/utility/util.jl | 17 ----------------- 2 files changed, 19 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 5529375b8..552450026 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -107,8 +107,6 @@ include("algorithms/truncation/bond_truncation.jl") include("algorithms/changebonds/local.jl") -include("algorithms/changebonds/local.jl") - include("algorithms/time_evolution/apply_gate.jl") include("algorithms/time_evolution/apply_mpo.jl") include("algorithms/time_evolution/trotter_gate.jl") diff --git a/src/utility/util.jl b/src/utility/util.jl index 47d3dde88..9ac564c44 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -128,23 +128,6 @@ function twistdual!(t::AbstractTensorMap, is) end twistdual(t::AbstractTensorMap, is) = twistdual!(copy(t), is) -# lifted from #311, to be removed once that's merged -""" - twistnondual(t::AbstractTensorMap, i) - twistnondual!(t::AbstractTensorMap, i) - -Twist the i-th leg of a tensor `t` if it represents a non-dual space. -""" -function twistnondual!(t::AbstractTensorMap, i::Int) - !isdual(space(t, i)) || return t - return twist!(t, i) -end -function twistnondual!(t::AbstractTensorMap, is) - is′ = filter(i -> !isdual(space(t, i)), is) - return twist!(t, is′) -end -twistnondual(t::AbstractTensorMap, is) = twistnondual!(copy(t), is) - """ twistnondual(t::AbstractTensorMap, i) twistnondual!(t::AbstractTensorMap, i) From 2ef5720290eb223fc846bd6ca960395ab5b2a1f8 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 19 May 2026 09:55:01 +0800 Subject: [PATCH 21/32] Rename to `compress` --- src/PEPSKit.jl | 4 ++-- .../{changebonds => compress}/local.jl | 20 +++++++++---------- test/{changebonds => compress}/local.jl | 2 +- test/timeevol/tf_ising_finiteT.jl | 4 ++-- 4 files changed, 14 insertions(+), 16 deletions(-) rename src/algorithms/{changebonds => compress}/local.jl (91%) rename test/{changebonds => compress}/local.jl (98%) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 1d30aebd1..ddef8397a 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -109,7 +109,7 @@ include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_tensor.jl") include("algorithms/truncation/bond_truncation.jl") -include("algorithms/changebonds/local.jl") +include("algorithms/compress/local.jl") include("algorithms/time_evolution/apply_gate.jl") include("algorithms/time_evolution/apply_mpo.jl") @@ -150,7 +150,7 @@ export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint export LocalTruncation -export changebonds +export compress export absorb_weight export ALSTruncation, FullEnvTruncation diff --git a/src/algorithms/changebonds/local.jl b/src/algorithms/compress/local.jl similarity index 91% rename from src/algorithms/changebonds/local.jl rename to src/algorithms/compress/local.jl index 7df2f8419..261d39dd6 100644 --- a/src/algorithms/changebonds/local.jl +++ b/src/algorithms/compress/local.jl @@ -1,10 +1,10 @@ """ $(TYPEDEF) -Abstract super type for algorithms to change virtual bonds -in two-dimensional tensor networks +Abstract super type for virtual space truncation algorithms when +compressing two or more layers of two-dimensional tensor networks """ -abstract type BondChangeAlgorithm end +abstract type CompressAlgorithm end """ $(TYPEDEF) @@ -23,7 +23,7 @@ which is the 2-norm of (e.g. for two layers of iPEPO) ``` on each bond of the network. """ -struct LocalTruncation <: BondChangeAlgorithm +struct LocalTruncation <: CompressAlgorithm trunc::TruncationStrategy end @@ -52,9 +52,8 @@ function _get_MdagM(A1::PEPOTensor, A2::PEPOTensor) @assert isdual(virtualspace(A1, EAST)) && isdual(virtualspace(A2, EAST)) A2′ = twistdual(A2, [2, 3, 5, 6]) A1′ = twistdual(A1, [1, 3, 5, 6]) - @tensoropt MdagM[x2 z z′ x2′] := - conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] - @tensoropt MdagM[x1 x2; x1′ x2′] := MdagM[x2 z z′ x2′] * + @tensoropt MdagM[x1 x2; x1′ x2′] := + conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] * conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] project_hermitian!(MdagM) return MdagM @@ -85,9 +84,8 @@ function _get_MMdag(A1::PEPOTensor, A2::PEPOTensor) @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) A2′ = twistnondual(A2, [2, 3, 4, 5]) A1′ = twistnondual(A1, [1, 3, 4, 5]) - @tensoropt MMdag[x2 z z′ x2′] := - A2′[z z2; Y2 X2 y2 x2] * conj(A2[z′ z2; Y2 X2 y2 x2′]) - @tensoropt MMdag[x1 x2; x1′ x2′] := MMdag[x2 z z′ x2′] * + @tensoropt MMdag[x1 x2; x1′ x2′] := + A2′[z z2; Y2 X2 y2 x2] * conj(A2[z′ z2; Y2 X2 y2 x2′]) * A1′[z1 z; Y1 X1 y1 x1] * conj(A1[z1 z′; Y1 X1 y1 x1′]) project_hermitian!(MMdag) return MMdag @@ -124,7 +122,7 @@ end Truncate virtual bonds of the product of two 1-layer InfinitePEPOs `ρ1`, `ρ2` with `LocalTruncation`. """ -function MPSKit.changebonds(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) +function compress(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) # sanity checks (size(ρ1) == size(ρ2)) || error("Input PEPOs have different unit cell sizes.") (size(ρ1, 3) == 1) || error("ρ1 should have only one layer.") diff --git a/test/changebonds/local.jl b/test/compress/local.jl similarity index 98% rename from test/changebonds/local.jl rename to test/compress/local.jl index a1853b7fd..5c14b1421 100644 --- a/test/changebonds/local.jl +++ b/test/compress/local.jl @@ -56,6 +56,6 @@ end Ves = ComplexSpace.([3 5; 4 2]) ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) alg = LocalTruncation(truncrank(2)) - ρ2, = changebonds(ρ, ρ, alg) + ρ2, = compress(ρ, ρ, alg) @test ρ2 isa InfinitePEPO end diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index 187cd9440..cc23c78bf 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -64,8 +64,8 @@ dt, nstep = 1.0e-3, 400 @test β ≈ info.t @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) - # use `changebonds` to reach 2β, or T = 1.25 - pepo2, = changebonds(pepo, pepo, LocalTruncation(trunc_pepo)) + # use `compress` to reach 2β, or T = 1.25 + pepo2, = compress(pepo, pepo, LocalTruncation(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) From fae86111fb63c288566ce5cdccbd3087ffc4c0e5 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 19 May 2026 10:20:29 +0800 Subject: [PATCH 22/32] Improve PSD projection + square root --- src/algorithms/compress/local.jl | 14 ++++++++------ src/utility/util.jl | 10 ---------- 2 files changed, 8 insertions(+), 16 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index 261d39dd6..141373c40 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -44,9 +44,10 @@ Only `R` is calculated and returned. function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) MdagM = _get_MdagM(A1, A2) D, R = eigh_full!(MdagM) - project_psd!(D) - R = sdiag_pow(D, 0.5) * R' - return R + # remove small negative eigenvalues due to numerical noises + T = eltype(D.data) + D.data .= sqrt.(max.(zero(T), D.data)) + return lmul!(D, R') end function _get_MdagM(A1::PEPOTensor, A2::PEPOTensor) @assert isdual(virtualspace(A1, EAST)) && isdual(virtualspace(A2, EAST)) @@ -76,9 +77,10 @@ Only `L` is calculated and returned. function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) MMdag = _get_MMdag(A1, A2) D, L = eigh_full!(MMdag) - project_psd!(D) - L = L * sdiag_pow(D, 0.5) - return L + # remove small negative eigenvalues due to numerical noises + T = eltype(D.data) + D.data .= sqrt.(max.(zero(T), D.data)) + return rmul!(L, D) end function _get_MMdag(A1::PEPOTensor, A2::PEPOTensor) @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) diff --git a/src/utility/util.jl b/src/utility/util.jl index cb7c1073d..87e3a65f4 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -48,16 +48,6 @@ function ChainRulesCore.rrule( return spow, sdiag_pow_pullback end -""" -Given the eigenvalue spectrum of `D` expected from a -positive-semi-definite (PSD) map, sets small negative -eigenvalues that may arise from numerical errors to 0. -""" -function project_psd!(D::DiagonalTensorMap{T, S}) where {T, S} - D.data .= max.(zero(T), D.data) - return D -end - """ absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorMap) From 38a8b04c93dfac181acd8f1494e0433d4c8d9065 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 19 May 2026 10:20:52 +0800 Subject: [PATCH 23/32] Reorder tests --- test/compress/local.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/test/compress/local.jl b/test/compress/local.jl index 5c14b1421..de7c6cbdf 100644 --- a/test/compress/local.jl +++ b/test/compress/local.jl @@ -17,6 +17,19 @@ function localcompress_cost(A1, A2, B1, B2, P1, P2) return norm(net1 - net2) end +@testset "Fermionic twists" begin + Vphy = Vect[FermionParity](0 => 2, 1 => 2) + Vvir = Vect[FermionParity](0 => 2, 1 => 2) + for _ in 1:4 # multiple trials without setting seed + Aspace = (Vphy ⊗ Vphy' ← Vvir ⊗ Vvir ⊗ Vvir' ⊗ Vvir') + A1 = randn(ComplexF64, Aspace) + A2 = randn(ComplexF64, Aspace) + for MM in [PEPSKit._get_MMdag(A1, A2), PEPSKit._get_MdagM(A1, A2)] + @test isposdef(MM) + end + end +end + @testset "Cost function of LocalTruncation" begin Random.seed!(0) Vaux, Vphy, V = ℂ^1, ℂ^10, ℂ^4 @@ -37,19 +50,6 @@ end @test info.ϵ ≈ localcompress_cost(A1, A2, B1, B2, P1, P2) end -@testset "Fermionic twists for qr/lq_twolayer" begin - Vphy = Vect[FermionParity](0 => 2, 1 => 2) - Vvir = Vect[FermionParity](0 => 2, 1 => 2) - for _ in 1:4 # multiple trials without setting seed - Aspace = (Vphy ⊗ Vphy' ← Vvir ⊗ Vvir ⊗ Vvir' ⊗ Vvir') - A1 = randn(ComplexF64, Aspace) - A2 = randn(ComplexF64, Aspace) - for MM in [PEPSKit._get_MMdag(A1, A2), PEPSKit._get_MdagM(A1, A2)] - @test isposdef(MM) - end - end -end - @testset "Virtual space matching" begin Vps = ComplexSpace.([2 2; 2 2]) Vns = ComplexSpace.([2 4; 5 3]) From 44aba113dafcc850bf94cffec7906ab2ae6d5baf Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 19 May 2026 10:23:46 +0800 Subject: [PATCH 24/32] Rename qr/lq_twolayer --- src/algorithms/compress/local.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index 141373c40..f635c2d45 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -28,7 +28,7 @@ struct LocalTruncation <: CompressAlgorithm end """ -Calculate the QR decomposition of 2-layer PEPO tensor +Calculate the `left_orth` of 2-layer PEPO tensor with the east virtual legs transferred to the R tensor ``` ↓ ╱ @@ -41,7 +41,7 @@ with the east virtual legs transferred to the R tensor ``` Only `R` is calculated and returned. """ -function qr_twolayer(A1::PEPOTensor, A2::PEPOTensor) +function left_orth_twolayer(A1::PEPOTensor, A2::PEPOTensor) MdagM = _get_MdagM(A1, A2) D, R = eigh_full!(MdagM) # remove small negative eigenvalues due to numerical noises @@ -61,7 +61,7 @@ function _get_MdagM(A1::PEPOTensor, A2::PEPOTensor) end """ -Calculate the LQ decomposition of 2-layer PEPO tensor +Calculate the `right_orth` of 2-layer PEPO tensor with the west virtual legs transferred to the L tensor ``` ↓ ╱ @@ -74,7 +74,7 @@ with the west virtual legs transferred to the L tensor ``` Only `L` is calculated and returned. """ -function lq_twolayer(A1::PEPOTensor, A2::PEPOTensor) +function right_orth_twolayer(A1::PEPOTensor, A2::PEPOTensor) MMdag = _get_MMdag(A1, A2) D, L = eigh_full!(MMdag) # remove small negative eigenvalues due to numerical noises @@ -111,8 +111,8 @@ function virtual_projector( A1::PEPOTensor, A2::PEPOTensor, B1::PEPOTensor, B2::PEPOTensor; trunc::TruncationStrategy ) - R1 = qr_twolayer(A1, A2) - R2 = lq_twolayer(B1, B2) + R1 = left_orth_twolayer(A1, A2) + R2 = right_orth_twolayer(B1, B2) u, s, vh, ϵ = svd_trunc!(R1 * R2; trunc) sinv_sqrt = sdiag_pow(s, -0.5) P1 = R2 * vh' * sinv_sqrt From 4d802135613f2a5a750fc97d2b13cd3306f8e790 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 20 May 2026 08:05:03 +0800 Subject: [PATCH 25/32] Remove CompressAlgorithm abstract type --- src/algorithms/compress/local.jl | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index f635c2d45..fa8099bc8 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -1,14 +1,6 @@ """ $(TYPEDEF) -Abstract super type for virtual space truncation algorithms when -compressing two or more layers of two-dimensional tensor networks -""" -abstract type CompressAlgorithm end - -""" -$(TYPEDEF) - Algorithm to truncate virtual bonds of a two-layer network with projectors that minimizes the local cost function, which is the 2-norm of (e.g. for two layers of iPEPO) @@ -23,7 +15,7 @@ which is the 2-norm of (e.g. for two layers of iPEPO) ``` on each bond of the network. """ -struct LocalTruncation <: CompressAlgorithm +struct LocalTruncation trunc::TruncationStrategy end From da824f41744df08b5ff7281e9be68fd9a95577a2 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 20 May 2026 08:09:19 +0800 Subject: [PATCH 26/32] Use autoopt for MdagM and MMdag --- src/algorithms/compress/local.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index fa8099bc8..be74d07e9 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -45,9 +45,9 @@ function _get_MdagM(A1::PEPOTensor, A2::PEPOTensor) @assert isdual(virtualspace(A1, EAST)) && isdual(virtualspace(A2, EAST)) A2′ = twistdual(A2, [2, 3, 5, 6]) A1′ = twistdual(A1, [1, 3, 5, 6]) - @tensoropt MdagM[x1 x2; x1′ x2′] := - conj(A2[z z2; Y2 x2 y2 X2]) * A2′[z′ z2; Y2 x2′ y2 X2] * - conj(A1[z1 z; Y1 x1 y1 X1]) * A1′[z1 z′; Y1 x1′ y1 X1] + @autoopt @tensor MdagM[Dx1 Dx2; Dx1′ Dx2′] := + conj(A2[dz dz2; DY2 Dx2 Dy2 DX2]) * A2′[dz′ dz2; DY2 Dx2′ Dy2 DX2] * + conj(A1[dz1 dz; DY1 Dx1 Dy1 DX1]) * A1′[dz1 dz′; DY1 Dx1′ Dy1 DX1] project_hermitian!(MdagM) return MdagM end @@ -78,9 +78,9 @@ function _get_MMdag(A1::PEPOTensor, A2::PEPOTensor) @assert !isdual(virtualspace(A1, WEST)) && !isdual(virtualspace(A2, WEST)) A2′ = twistnondual(A2, [2, 3, 4, 5]) A1′ = twistnondual(A1, [1, 3, 4, 5]) - @tensoropt MMdag[x1 x2; x1′ x2′] := - A2′[z z2; Y2 X2 y2 x2] * conj(A2[z′ z2; Y2 X2 y2 x2′]) * - A1′[z1 z; Y1 X1 y1 x1] * conj(A1[z1 z′; Y1 X1 y1 x1′]) + @autoopt @tensor MMdag[Dx1 Dx2; Dx1′ Dx2′] := + A2′[dz dz2; DY2 DX2 Dy2 Dx2] * conj(A2[dz′ dz2; DY2 DX2 Dy2 Dx2′]) * + A1′[dz1 dz; DY1 DX1 Dy1 Dx1] * conj(A1[dz1 dz′; DY1 DX1 Dy1 Dx1′]) project_hermitian!(MMdag) return MMdag end From 5ea481f8ed2975af63dcbd42f35cf519ab1df18b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 2 Jul 2026 11:47:04 +0800 Subject: [PATCH 27/32] Cartesian `eachindex` for InfinitePEPS / PEPO / PartitionFunction --- src/operators/infinitepepo.jl | 7 ++++--- src/states/infinitepartitionfunction.jl | 1 + src/states/infinitepeps.jl | 1 + test/bp/rotation.jl | 2 +- test/types/infinitepartitionfunction.jl | 1 + test/types/infinitesquarenetwork.jl | 3 +++ test/types/suweight.jl | 8 ++++---- 7 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index b15b70d11..d7ea5b1a9 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -143,6 +143,7 @@ Base.@propagate_inbounds Base.getindex(A::InfinitePEPO, I...) = Base.@propagate_inbounds Base.setindex!(A::InfinitePEPO, v, I...) = periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePEPO, args...) = axes(unitcell(A), args...) +Base.eachindex(A::InfinitePEPO) = eachindex(IndexCartesian(), unitcell(A)) eachcoordinate(A::InfinitePEPO) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePEPO, dirs) return collect(Iterators.product(dirs, axes(A, 1), axes(A, 2))) @@ -263,11 +264,11 @@ end Checks if `state` follow the standard virtual arrow convention, i.e. north and east domains are non-dual spaces. """ -function _check_virtual_dualness(state::Union{InfinitePEPS, InfinitePEPO}) - isdual_easts = map(CartesianIndices(state.A)) do idx +function _check_virtual_dualness(state::InfiniteState) + isdual_easts = map(eachindex(state)) do idx return isdual(virtualspace(state[idx], EAST)) end - isdual_norths = map(CartesianIndices(state.A)) do idx + isdual_norths = map(eachindex(state)) do idx return isdual(virtualspace(state[idx], NORTH)) end return isdual_easts, isdual_norths diff --git a/src/states/infinitepartitionfunction.jl b/src/states/infinitepartitionfunction.jl index 8adb2517f..b685bdf90 100644 --- a/src/states/infinitepartitionfunction.jl +++ b/src/states/infinitepartitionfunction.jl @@ -131,6 +131,7 @@ Base.@propagate_inbounds Base.getindex(A::InfinitePartitionFunction, I...) = Base.@propagate_inbounds Base.setindex!(A::InfinitePartitionFunction, v, I...) = periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePartitionFunction, args...) = axes(unitcell(A), args...) +Base.eachindex(A::InfinitePartitionFunction) = eachindex(IndexCartesian(), unitcell(A)) eachcoordinate(A::InfinitePartitionFunction) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePartitionFunction, dirs) return collect(Iterators.product(dirs, axes(A, 1), axes(A, 2))) diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index f4c147549..d25b19e7a 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -136,6 +136,7 @@ Base.@propagate_inbounds Base.getindex(A::InfinitePEPS, I...) = Base.@propagate_inbounds Base.setindex!(A::InfinitePEPS, v, I...) = periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePEPS, args...) = axes(unitcell(A), args...) +Base.eachindex(A::InfinitePEPS) = eachindex(IndexCartesian(), unitcell(A)) eachcoordinate(A::InfinitePEPS) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePEPS, dirs) return collect(Iterators.product(dirs, axes(A, 1), axes(A, 2))) diff --git a/test/bp/rotation.jl b/test/bp/rotation.jl index 8b6fc7428..6effb1898 100644 --- a/test/bp/rotation.jl +++ b/test/bp/rotation.jl @@ -20,7 +20,7 @@ function meas_sites( op::O, ψ::InfinitePEPS, env::Union{BPEnv, CTMRGEnv} ) where {O <: AbstractTensorMap{<:Any, <:Any, 1, 1}} lattice = physicalspace(ψ) - return map(CartesianIndices(ψ.A)) do site1 + return map(eachindex(ψ)) do site1 lo = LocalOperator(lattice, (site1,) => op) return expectation_value(ψ, lo, env) end diff --git a/test/types/infinitepartitionfunction.jl b/test/types/infinitepartitionfunction.jl index 19be8e02b..0e653fe2e 100644 --- a/test/types/infinitepartitionfunction.jl +++ b/test/types/infinitepartitionfunction.jl @@ -24,4 +24,5 @@ sizes = [(1, 1), (3, 3)] @test (rotl90 ∘ rotl90)(pf) ≈ rot180(pf) @test (rotr90 ∘ rotr90 ∘ rotr90)(pf) ≈ rotl90(pf) @test length(pf) == prod(sz) + @test eachindex(pf) == CartesianIndices(size(pf)) end diff --git a/test/types/infinitesquarenetwork.jl b/test/types/infinitesquarenetwork.jl index ff6b57b53..9b6e35d4f 100644 --- a/test/types/infinitesquarenetwork.jl +++ b/test/types/infinitesquarenetwork.jl @@ -25,6 +25,9 @@ sizes = [(1, 1), (3, 3)] pf = InfinitePartitionFunction(pf_tensor; unitcell = sz) pepo = InfinitePEPO(pepo_tensor; unitcell = (sz..., 2)) + @test eachindex(peps) == CartesianIndices(size(peps)) + @test eachindex(pepo) == CartesianIndices(size(pepo)) + peps_n = InfiniteSquareNetwork(peps) pf_n = InfiniteSquareNetwork(pf) pepo_n = InfiniteSquareNetwork(peps, pepo) diff --git a/test/types/suweight.jl b/test/types/suweight.jl index cc7440920..b65b8e705 100644 --- a/test/types/suweight.jl +++ b/test/types/suweight.jl @@ -21,16 +21,16 @@ end function test_rotation(state::Union{InfinitePEPS, InfinitePEPO}, wts::SUWeight) InfiniteState = (state isa InfinitePEPS) ? InfinitePEPS : InfinitePEPO A1 = InfiniteState( - map(CartesianIndices(state.A)) do idx - return absorb_weight(state.A[idx], wts, idx[1], idx[2], Tuple(1:4)) + map(eachindex(state)) do idx + return absorb_weight(state[idx], wts, idx[1], idx[2], Tuple(1:4)) end ) for n in 1:4 rot = compose_n(rotl90, n) state2, wts2 = rot(state), rot(wts) A2 = InfiniteState( - map(CartesianIndices(state2.A)) do idx - return absorb_weight(state2.A[idx], wts2, idx[1], idx[2], Tuple(1:4)) + map(eachindex(state2)) do idx + return absorb_weight(state2[idx], wts2, idx[1], idx[2], Tuple(1:4)) end ) @test A2 ≈ rot(A1) From 95cd0dfd7debb331eae0d6c9f70ed40f2a30df9c Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 2 Jul 2026 11:58:27 +0800 Subject: [PATCH 28/32] Reduce repeated rotations --- src/algorithms/compress/local.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index be74d07e9..767f77ddc 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -133,11 +133,10 @@ function compress(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] + ρ1′, ρ2′ = rotr90.(unitcell(ρ1)), rotr90.(unitcell(ρ2)) Pys_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) - # TODO: reduce repeated rotations return virtual_projector( - rotr90(ρ1[r, c]), rotr90(ρ2[r, c]), - rotr90(ρ1[_prev(r, Nr), c]), rotr90(ρ2[_prev(r, Nr), c]); + ρ1′[r, c], ρ2′[r, c], ρ1′[_prev(r, Nr), c], ρ2′[_prev(r, Nr), c]; trunc = alg.trunc ) end From 0505fddbc74fb1f376511c9c4c209d48328fa58d Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 2 Jul 2026 12:34:41 +0800 Subject: [PATCH 29/32] Lazy multiplication of iPEPOs --- src/algorithms/compress/local.jl | 23 +++++++++++------------ src/operators/infinitepepo.jl | 13 +++++++++++++ test/compress/local.jl | 2 +- test/timeevol/tf_ising_finiteT.jl | 2 +- test/types/infinitesquarenetwork.jl | 10 ++++++++++ 5 files changed, 36 insertions(+), 14 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index 767f77ddc..5e55c6cd9 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -113,27 +113,26 @@ function virtual_projector( end """ -Truncate virtual bonds of the product of two 1-layer -InfinitePEPOs `ρ1`, `ρ2` with `LocalTruncation`. + compress(ρ::InfinitePEPO, alg::LocalTruncation) + +Compress a 2-layer iPEPO into a 1-layer iPEPO by +truncating the virtual bonds with `LocalTruncation`. """ -function compress(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) +function compress(ρ::InfinitePEPO, alg::LocalTruncation) # sanity checks - (size(ρ1) == size(ρ2)) || error("Input PEPOs have different unit cell sizes.") - (size(ρ1, 3) == 1) || error("ρ1 should have only one layer.") - (size(ρ2, 3) == 1) || error("ρ2 should have only one layer.") - all(all.(_check_virtual_dualness(ρ1))) || error("East and north virtual spaces in ρ1 should be dual spaces.") - all(all.(_check_virtual_dualness(ρ2))) || error("East and north virtual spaces in ρ2 should be dual spaces.") + size(ρ, 3) == 2 || error("Input PEPO should have two layers.") + all(all.(_check_virtual_dualness(ρ))) || error("East and north virtual spaces in ρ should be dual spaces.") # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] - Nr, Nc, = size(ρ1) + Nr, Nc, = size(ρ) Pxs_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) # TODO: support SiteDependentTruncation return virtual_projector( - ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; + ρ[r, c, 2], ρ[r, c, 1], ρ[r, _next(c, Nc), 2], ρ[r, _next(c, Nc), 1]; trunc = alg.trunc ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] - ρ1′, ρ2′ = rotr90.(unitcell(ρ1)), rotr90.(unitcell(ρ2)) + ρ1′, ρ2′ = rotr90.(unitcell(ρ)[:, :, 2]), rotr90.(unitcell(ρ)[:, :, 1]) Pys_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) return virtual_projector( ρ1′[r, c], ρ2′[r, c], ρ1′[_prev(r, Nr), c], ρ2′[_prev(r, Nr), c]; @@ -145,7 +144,7 @@ function compress(ρ1::InfinitePEPO, ρ2::InfinitePEPO, alg::LocalTruncation) Pw, Pe = Pxs_info[r, _prev(c, Nc)][2], Pxs_info[r, c][1] Pn, Ps = Pys_info[r, c][1], Pys_info[_next(r, Nr), c][2] @tensoropt A[p1 p2; n e s w] := - (ρ1[r, c])[p1 p; n1 e1 s1 w1] * (ρ2[r, c])[p p2; n2 e2 s2 w2] * + (ρ[r, c, 2])[p1 p; n1 e1 s1 w1] * (ρ[r, c, 1])[p p2; n2 e2 s2 w2] * Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index d7ea5b1a9..e073e46a6 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -213,6 +213,19 @@ end VI.scalartype(::Type{NT}) where {NT <: InfinitePEPO} = scalartype(eltype(NT)) VI.zerovector(A::InfinitePEPO) = InfinitePEPO(zerovector(unitcell(A))) +## Math + +""" + ρ1 * ρ2 + +Lazily multiply two infinite PEPOs by stacking the layers of `ρ1` under those of `ρ2`. +""" +function Base.:*(ρ1::InfinitePEPO, ρ2::InfinitePEPO) + size(ρ1)[1:2] == size(ρ2)[1:2] || + throw(ArgumentError("Input PEPOs have different unit cell sizes.")) + return InfinitePEPO(cat(unitcell(ρ2), unitcell(ρ1); dims = 3)) +end + ## (Approximate) equality function Base.:(==)(A₁::InfinitePEPO, A₂::InfinitePEPO) return all(zip(unitcell(A₁), unitcell(A₂))) do (p₁, p₂) diff --git a/test/compress/local.jl b/test/compress/local.jl index de7c6cbdf..d83ab7885 100644 --- a/test/compress/local.jl +++ b/test/compress/local.jl @@ -56,6 +56,6 @@ end Ves = ComplexSpace.([3 5; 4 2]) ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) alg = LocalTruncation(truncrank(2)) - ρ2, = compress(ρ, ρ, alg) + ρ2, = compress(ρ * ρ, alg) @test ρ2 isa InfinitePEPO end diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index 52accf8ed..b713fc8f1 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -65,7 +65,7 @@ dt, nstep = 1.0e-3, 400 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) # use `compress` to reach 2β, or T = 1.25 - pepo2, = compress(pepo, pepo, LocalTruncation(trunc_pepo)) + pepo2, = compress(pepo * pepo, LocalTruncation(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) diff --git a/test/types/infinitesquarenetwork.jl b/test/types/infinitesquarenetwork.jl index 9b6e35d4f..51283a7b5 100644 --- a/test/types/infinitesquarenetwork.jl +++ b/test/types/infinitesquarenetwork.jl @@ -20,10 +20,12 @@ sizes = [(1, 1), (3, 3)] peps_tensor = PEPSTensor(randn, T, P, V) pf_tensor = PFTensor(randn, T, V) pepo_tensor = randn(T, P ⊗ P' ← V ⊗ V ⊗ V' ⊗ V') + pepo_tensor2 = randn(T, P ⊗ P' ← V ⊗ V ⊗ V' ⊗ V') peps = InfinitePEPS(peps_tensor; unitcell = sz) pf = InfinitePartitionFunction(pf_tensor; unitcell = sz) pepo = InfinitePEPO(pepo_tensor; unitcell = (sz..., 2)) + pepo2 = InfinitePEPO(pepo_tensor2; unitcell = (sz..., 1)) @test eachindex(peps) == CartesianIndices(size(peps)) @test eachindex(pepo) == CartesianIndices(size(pepo)) @@ -54,5 +56,13 @@ sizes = [(1, 1), (3, 3)] @test pf_n + pf_n ≈ 2 * pf_n @test (rotr90 ∘ rotr90)(pf_n) ≈ rot180(pf_n) + pepo_product = pepo2 * pepo + @test size(pepo_product) == (sz..., 3) + @test all(pepo_product[I] == pepo[I] for I in CartesianIndices(size(pepo))) + @test all( + pepo_product[r, c, size(pepo, 3) + h] == pepo2[r, c, h] for + r in axes(pepo2, 1), c in axes(pepo2, 2), h in axes(pepo2, 3) + ) + @test length(peps_n) == prod(sz) end From 4aa2833733d8efa434ac295bc338ea4e3d9e3909 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 7 Jul 2026 17:18:04 +0800 Subject: [PATCH 30/32] Use iPEPO tuple for `compress` --- src/algorithms/compress/local.jl | 25 +++++++++++++++---------- test/compress/local.jl | 2 +- test/timeevol/tf_ising_finiteT.jl | 2 +- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index 5e55c6cd9..a66a83155 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -113,26 +113,31 @@ function virtual_projector( end """ - compress(ρ::InfinitePEPO, alg::LocalTruncation) + compress((ρ1, ρ2), alg::LocalTruncation) -Compress a 2-layer iPEPO into a 1-layer iPEPO by -truncating the virtual bonds with `LocalTruncation`. +Compress two 1-layer iPEPOs into a 1-layer iPEPO by truncating the virtual bonds +with `LocalTruncation`. In the tuple `(ρ1, ρ2)`, `ρ1` is the lower layer and +`ρ2` is the upper layer. """ -function compress(ρ::InfinitePEPO, alg::LocalTruncation) +function compress(ρs::Tuple{InfinitePEPO, InfinitePEPO}, alg::LocalTruncation) + ρ1, ρ2 = ρs # sanity checks - size(ρ, 3) == 2 || error("Input PEPO should have two layers.") - all(all.(_check_virtual_dualness(ρ))) || error("East and north virtual spaces in ρ should be dual spaces.") + size(ρ1) == size(ρ2) || error("Input PEPOs have different unit cell sizes.") + size(ρ1, 3) == 1 || error("ρ1 should have only one layer.") + size(ρ2, 3) == 1 || error("ρ2 should have only one layer.") + all(all.(_check_virtual_dualness(ρ1))) || error("East and north virtual spaces in ρ1 should be dual spaces.") + all(all.(_check_virtual_dualness(ρ2))) || error("East and north virtual spaces in ρ2 should be dual spaces.") # x-bond projectors: [r, c] on bond [r, c]--[r, c+1] - Nr, Nc, = size(ρ) + Nr, Nc, = size(ρ1) Pxs_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) # TODO: support SiteDependentTruncation return virtual_projector( - ρ[r, c, 2], ρ[r, c, 1], ρ[r, _next(c, Nc), 2], ρ[r, _next(c, Nc), 1]; + ρ1[r, c], ρ2[r, c], ρ1[r, _next(c, Nc)], ρ2[r, _next(c, Nc)]; trunc = alg.trunc ) end # y-bond projectors: [r, c] on bond [r, c]--[r-1, c] - ρ1′, ρ2′ = rotr90.(unitcell(ρ)[:, :, 2]), rotr90.(unitcell(ρ)[:, :, 1]) + ρ1′, ρ2′ = rotr90.(unitcell(ρ1)), rotr90.(unitcell(ρ2)) Pys_info = map(Iterators.product(1:Nr, 1:Nc)) do (r, c) return virtual_projector( ρ1′[r, c], ρ2′[r, c], ρ1′[_prev(r, Nr), c], ρ2′[_prev(r, Nr), c]; @@ -144,7 +149,7 @@ function compress(ρ::InfinitePEPO, alg::LocalTruncation) Pw, Pe = Pxs_info[r, _prev(c, Nc)][2], Pxs_info[r, c][1] Pn, Ps = Pys_info[r, c][1], Pys_info[_next(r, Nr), c][2] @tensoropt A[p1 p2; n e s w] := - (ρ[r, c, 2])[p1 p; n1 e1 s1 w1] * (ρ[r, c, 1])[p p2; n2 e2 s2 w2] * + (ρ1[r, c])[p1 p; n1 e1 s1 w1] * (ρ2[r, c])[p p2; n2 e2 s2 w2] * Pn[n1 n2; n] * Pe[e1 e2; e] * Ps[s; s1 s2] * Pw[w; w1 w2] return A end diff --git a/test/compress/local.jl b/test/compress/local.jl index d83ab7885..0817378cd 100644 --- a/test/compress/local.jl +++ b/test/compress/local.jl @@ -56,6 +56,6 @@ end Ves = ComplexSpace.([3 5; 4 2]) ρ = InfinitePEPO(randn, ComplexF64, Vps, Vns, Ves) alg = LocalTruncation(truncrank(2)) - ρ2, = compress(ρ * ρ, alg) + ρ2, = compress((ρ, ρ), alg) @test ρ2 isa InfinitePEPO end diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index b713fc8f1..573b84e90 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/test/timeevol/tf_ising_finiteT.jl @@ -65,7 +65,7 @@ dt, nstep = 1.0e-3, 400 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2) # use `compress` to reach 2β, or T = 1.25 - pepo2, = compress(pepo * pepo, LocalTruncation(trunc_pepo)) + pepo2, = compress((pepo, pepo), LocalTruncation(trunc_pepo)) normalize!.(pepo2.A) env2 = converge_env(InfinitePartitionFunction(pepo2), 16) result_2β = measure_mag(pepo2, env2) From 4598411f96c7b4cbe79a5157463eddd543041a20 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Jul 2026 11:25:01 +0800 Subject: [PATCH 31/32] Fix type annotation --- src/algorithms/compress/local.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl index a66a83155..3403b6e0d 100644 --- a/src/algorithms/compress/local.jl +++ b/src/algorithms/compress/local.jl @@ -119,7 +119,7 @@ Compress two 1-layer iPEPOs into a 1-layer iPEPO by truncating the virtual bonds with `LocalTruncation`. In the tuple `(ρ1, ρ2)`, `ρ1` is the lower layer and `ρ2` is the upper layer. """ -function compress(ρs::Tuple{InfinitePEPO, InfinitePEPO}, alg::LocalTruncation) +function compress(ρs::Tuple{<:InfinitePEPO, <:InfinitePEPO}, alg::LocalTruncation) ρ1, ρ2 = ρs # sanity checks size(ρ1) == size(ρ2) || error("Input PEPOs have different unit cell sizes.") From 3fd8da3832372e33ebfde1f70e266bb5b9567a72 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Jul 2026 11:26:01 +0800 Subject: [PATCH 32/32] Reorder functions --- src/operators/infinitepepo.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index e073e46a6..5148ffe07 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -208,6 +208,20 @@ function _is_bipartite(psi::InfiniteState) return true end +""" +Checks if `state` follow the standard virtual arrow convention, +i.e. north and east domains are non-dual spaces. +""" +function _check_virtual_dualness(state::InfiniteState) + isdual_easts = map(eachindex(state)) do idx + return isdual(virtualspace(state[idx], EAST)) + end + isdual_norths = map(eachindex(state)) do idx + return isdual(virtualspace(state[idx], NORTH)) + end + return isdual_easts, isdual_norths +end + ## Vector interface VI.scalartype(::Type{NT}) where {NT <: InfinitePEPO} = scalartype(eltype(NT)) @@ -272,17 +286,3 @@ function _stack_tuples(A::Matrix{NTuple{N, T}}) where {N, T} end return out end - -""" -Checks if `state` follow the standard virtual arrow convention, -i.e. north and east domains are non-dual spaces. -""" -function _check_virtual_dualness(state::InfiniteState) - isdual_easts = map(eachindex(state)) do idx - return isdual(virtualspace(state[idx], EAST)) - end - isdual_norths = map(eachindex(state)) do idx - return isdual(virtualspace(state[idx], NORTH)) - end - return isdual_easts, isdual_norths -end