Compress 2-layer InfinitePEPO with local truncation

src/PEPSKit.jl

@@ -112,6 +112,8 @@ include("algorithms/truncation/fullenv_truncation.jl")
 include("algorithms/truncation/bond_tensor.jl")
 include("algorithms/truncation/bond_truncation.jl")
 
+include("algorithms/compress/local.jl")
+
 include("algorithms/time_evolution/apply_gate.jl")
 include("algorithms/time_evolution/apply_mpo.jl")
 include("algorithms/time_evolution/get_cluster.jl")
@@ -152,6 +154,9 @@ export leading_boundary
 export PEPSOptimize, FixedPointGradient, GeomSum, ManualIter
 export fixedpoint
 
+export LocalTruncation
+export compress
+
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
 export SimpleUpdate

src/algorithms/compress/local.jl

@@ -0,0 +1,158 @@
+"""
+$(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)
+```
+        ↓ ╱      ↓ ╱            ↓ ╱                 ↓ ╱
+    ----A2---←---B2---      ----A2-←-|╲       ╱|--←-B2---
+      ╱ |      ╱ |            ╱ |    | ╲     ╱ |  ╱ |
+        ↓        ↓       -      ↓    |P1├-←-┤P2|    ↓
+        | ╱      | ╱            | ╱  | ╱     ╲ |    | ╱
+    ----A1---←---B1---      ----A1-←-|╱       ╲|--←-B1---
+      ╱ ↓      ╱ ↓            ╱ ↓                 ╱ ↓
+```
+on each bond of the network.
+"""
+struct LocalTruncation
+    trunc::TruncationStrategy
+end
+
+"""
+Calculate the `left_orth` of 2-layer PEPO tensor
+with the east virtual legs transferred to the R tensor
+```
+        ↓ ╱
+    ----A2-←-           ┌-←-
+      ╱ |               |
+        ↓     =  (Q)-←--R
+        | ╱             |
+    ----A1-←-           └-←-
+      ╱ ↓
+```
+Only `R` is calculated and returned.
+"""
+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
+    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))
+    A2′ = twistdual(A2, [2, 3, 5, 6])
+    A1′ = twistdual(A1, [1, 3, 5, 6])
+    @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
+
+"""
+Calculate the `right_orth` of 2-layer PEPO tensor
+with the west virtual legs transferred to the L tensor
+```
+        ↓ ╱  
+    --←-A2---    -←-┐
+      ╱ |           |
+        ↓      =    L--←-(Q)
+        | ╱         |
+    --←-A1---    -←-┘
+      ╱ ↓    
+```
+Only `L` is calculated and returned.
+"""
+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
+    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))
+    A2′ = twistnondual(A2, [2, 3, 4, 5])
+    A1′ = twistnondual(A1, [1, 3, 4, 5])
+    @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
+
+"""
+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 virtual_projector(
+        A1::PEPOTensor, A2::PEPOTensor, B1::PEPOTensor, B2::PEPOTensor;
+        trunc::TruncationStrategy
+    )
+    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
+    P2 = sinv_sqrt * u' * R1
+    return P1, P2, (; s, ϵ)
+end
+
+"""
+    compress((ρ1, ρ2), alg::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(ρs::Tuple{<:InfinitePEPO, <:InfinitePEPO}, alg::LocalTruncation)
+    ρ1, ρ2 = ρs
+    # 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]
+    Nr, Nc, = size(ρ1)
+    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
+        )
+    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)
+        return virtual_projector(
+            ρ1′[r, c], ρ2′[r, c], ρ1′[_prev(r, Nr), c], ρ2′[_prev(r, Nr), c];
+            trunc = alg.trunc
+        )
+    end
+    # apply projectors
+    As = map(Iterators.product(1:Nr, 1:Nc)) do (r, c)
+        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 = (; Pxs_info, Pys_info)
+    return InfinitePEPO(cat(As; dims = 3)), info
+end

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)))
@@ -207,11 +208,38 @@ 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))
 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₂)

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

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

src/utility/util.jl

@@ -80,7 +80,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)

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

test/compress/local.jl

@@ -0,0 +1,61 @@
+using Test
+using Random
+using LinearAlgebra
+using TensorKit
+using PEPSKit
+using PEPSKit: virtual_projector
+
+"""
+Cost function of LocalTruncation.
+For test convenience, open virtual indices are made trivial and removed.
+"""
+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] *
+        A1[pa1 pa; Da1] * A2[pa pa2′; Da2] * B1[pb1 pb; Db1] * B2[pb pb2′; Db2]
+    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
+    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, P2, info = virtual_projector(A1, A2, B1, B2; trunc = notrunc())
+    @test P1 * P2 ≈ TensorKit.id(domain(P2))
+
+    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 = $(info.ϵ)."
+    @test info.ϵ ≈ localcompress_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 = LocalTruncation(truncrank(2))
+    ρ2, = compress((ρ, ρ), alg)
+    @test ρ2 isa InfinitePEPO
+end

test/timeevol/tf_ising_finiteT.jl

@@ -64,13 +64,13 @@ dt, nstep = 1.0e-3, 400
     @test β ≈ info.t
     @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
-    # 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)
+    # 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)
     @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)
+    @test isapprox(abs.(result_2β), bm_2β, rtol = 5.0e-3)
 
     # Purification approach: results at 2β, or T = 1.25
     alg = SimpleUpdate(; trunc = trunc_pepo, purified = true, bipartite, force_mpo)

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

test/types/infinitesquarenetwork.jl

@@ -20,10 +20,15 @@ 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))
 
     peps_n = InfiniteSquareNetwork(peps)
     pf_n = InfiniteSquareNetwork(pf)
@@ -51,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

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)