diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index a6e199c54..cd0974636 100644 --- a/src/PEPSKit.jl +++ b/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 diff --git a/src/algorithms/compress/local.jl b/src/algorithms/compress/local.jl new file mode 100644 index 000000000..3403b6e0d --- /dev/null +++ b/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 diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index a180a953c..5148ffe07 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))) @@ -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₂) 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/src/utility/util.jl b/src/utility/util.jl index 2c58f8801..87e3a65f4 100644 --- a/src/utility/util.jl +++ b/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) 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/compress/local.jl b/test/compress/local.jl new file mode 100644 index 000000000..0817378cd --- /dev/null +++ b/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 diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl index e1137cc04..573b84e90 100644 --- a/test/timeevol/tf_ising_finiteT.jl +++ b/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) 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..51283a7b5 100644 --- a/test/types/infinitesquarenetwork.jl +++ b/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 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)