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 diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4a39d6bd3..5ef6bbf32 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") @@ -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 diff --git a/src/algorithms/approximate/approx_tools.jl b/src/algorithms/approximate/approx_tools.jl new file mode 100644 index 000000000..f34ebc199 --- /dev/null +++ b/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 diff --git a/src/algorithms/approximate/local_approx.jl b/src/algorithms/approximate/local_approx.jl new file mode 100644 index 000000000..6af200401 --- /dev/null +++ b/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 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/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) 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")