[WIP] Local truncation for two layers of InfinitePEPO

.github/workflows/Tests.yml

@@ -34,6 +34,7 @@ jobs:
           - utility
           - bondenv
           - timeevol
+          - approx
           - toolbox
         os:
           - ubuntu-latest

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")
@@ -103,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

src/algorithms/approximate/approx_tools.jl

@@ -0,0 +1,62 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for the truncation algorithms of
+virtual spaces in InfinitePEPO or InfinitePEPS.
+"""
+abstract type ApproximateAlgorithm 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 _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.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 _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)
+end

src/algorithms/approximate/local_approx.jl

@@ -0,0 +1,125 @@
+struct LocalApprox <: ApproximateAlgorithm
+    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(cat(As; dims = 3))
+end

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

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)

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")