From 5204349ed999483403452971c73373a859966ffa Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 5 Dec 2025 17:27:00 +0800 Subject: [PATCH 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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