From c101c796c78f0200db5d758a996eb33a0be6518e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 6 Apr 2026 17:36:13 +0800 Subject: [PATCH 01/27] Change ALS tensor axis order --- .../contractions/bondenv/als_solve.jl | 113 +++++++----------- .../contractions/bondenv/gaugefix.jl | 6 +- src/algorithms/time_evolution/apply_mpo.jl | 55 +++++---- src/algorithms/truncation/bond_truncation.jl | 19 +-- .../truncation/fullenv_truncation.jl | 21 +--- test/bondenv/bond_truncate.jl | 6 +- 6 files changed, 91 insertions(+), 129 deletions(-) diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index c244c7370..af083f618 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -19,11 +19,9 @@ Construct the tensor └-----------------------------------┘ ``` """ -function _tensor_Ra( - benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} +function _tensor_Ra(benv::BondEnv, b::MPSTensor) return @autoopt @tensor Ra[DX1 Db1; DX0 Db0] := ( - benv[DX1 DY1; DX0 DY0] * b[Db0 DY0; db] * conj(b[Db1 DY1; db]) + benv[DX1 DY1; DX0 DY0] * b[Db0 db; DY0] * conj(b[Db1 db; DY1]) ) end @@ -44,10 +42,10 @@ Construct the tensor ``` """ function _tensor_Sa( - benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, b::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sa[DX1 Db1; da] := ( - benv[DX1 DY1; DX0 DY0] * conj(b[Db1 DY1; db]) * a2b2[DX0 DY0; da db] + return @autoopt @tensor Sa[DX1 da; Db1] := ( + benv[DX1 DY1; DX0 DY0] * conj(b[Db1 db; DY1]) * a2b2[DX0 DY0; da db] ) end @@ -67,11 +65,9 @@ Construct the tensor └-----------------------------------┘ ``` """ -function _tensor_Rb( - benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} +function _tensor_Rb(benv::BondEnv, a::MPSTensor) return @autoopt @tensor Rb[Da1 DY1; Da0 DY0] := ( - benv[DX1 DY1; DX0 DY0] * a[DX0 Da0; da] * conj(a[DX1 Da1; da]) + benv[DX1 DY1; DX0 DY0] * a[DX0 da; Da0] * conj(a[DX1 da; Da1]) ) end @@ -92,10 +88,10 @@ Construct the tensor ``` """ function _tensor_Sb( - benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, a::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sb[Da1 DY1; db] := ( - benv[DX1 DY1; DX0 DY0] * conj(a[DX1 Da1; da]) * a2b2[DX0 DY0; da db] + return @autoopt @tensor Sb[Da1 db; DY1] := ( + benv[DX1 DY1; DX0 DY0] * conj(a[DX1 da; Da1]) * a2b2[DX0 DY0; da db] ) end @@ -116,30 +112,10 @@ Calculate the inner product ``` """ function inner_prod( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * - conj(a1b1[DX1 DY1; da db]) * - a2b2[DX0 DY0; da db] -end - -""" -$(SIGNATURES) - -Calculate the fidelity between two evolution steps -``` - |⟨a1,b1|a2,b2⟩|^2 - -------------------------- - ⟨a1,b1|a1,b1⟩⟨a2,b2|a2,b2⟩ -``` -""" -function fidelity( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - b12 = inner_prod(benv, a1b1, a2b2) - b11 = inner_prod(benv, a1b1, a1b1) - b22 = inner_prod(benv, a2b2, a2b2) - return abs2(b12) / abs(b11 * b22) + conj(a1b1[DX1 DY1; da db]) * a2b2[DX0 DY0; da db] end """ @@ -153,14 +129,12 @@ Contract the axis between `a` and `b` tensors ``` """ function _combine_ab( - a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 1, 2} + a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2} ) where {T <: Number, S <: ElementarySpace} return @tensor ab[DX DY; da db] := a[DX da; D] * b[D; db DY] end -function _combine_ab( - a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} - return @tensor ab[DX DY; da db] := a[DX D; da] * b[D DY; db] +function _combine_ab(a::MPSTensor, b::MPSTensor) + return @tensor ab[DX DY; da db] := a[DX da; D] * b[D db; DY] end """ @@ -168,41 +142,40 @@ $(SIGNATURES) Calculate the cost function ``` - f(a,b) = ‖ |a1,b1⟩ - |a2,b2⟩ ‖^2 - = ⟨a1,b1|a1,b1⟩ - 2 Re⟨a1,b1|a2,b2⟩ + ⟨a2,b2|a2,b2⟩ + f(a,b) = ‖ |ψ1⟩ - |ψ2⟩ ‖^2 + = ⟨ψ1|benv|ψ1⟩ - 2 Re⟨ψ1|benv|ψ2⟩ + ⟨ψ2|benv|ψ2⟩ +``` +and the fidelity +``` + |⟨ψ1|benv|ψ2⟩|² + ------------------------ + ⟨ψ1|benv|ψ1⟩⟨ψ2|benv|ψ2⟩ ``` """ -function cost_function_als( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - t1 = inner_prod(benv, a1b1, a1b1) - t2 = inner_prod(benv, a2b2, a2b2) - t3 = inner_prod(benv, a1b1, a2b2) - return real(t1) + real(t2) - 2 * real(t3) +function cost_function_als(benv, ψ1, ψ2) + b12 = inner_prod(benv, ψ1, ψ2) + b11 = inner_prod(benv, ψ1, ψ1) + b22 = inner_prod(benv, ψ2, ψ2) + cost = real(b11) + real(b22) - 2 * real(b12) + fid = abs2(b12) / abs(b11 * b22) + return cost, fid end """ $(SIGNATURES) -Solve the equations `Rx x = Sx` (x = a, b) with initial guess `x0` -``` - ┌---------------------------┐ - | ┌----┐ | - └---| |--- 1 -- x -- 2 --┘ - | | ↓ - | Rx | -3 - | | - ┌---| |--- -1 -2 --┐ - | └----┘ | - └---------------------------┘ -``` -""" -function _solve_ab( - Rx::AbstractTensorMap{T, S, 2, 2}, - Sx::AbstractTensorMap{T, S, 2, 1}, - x0::AbstractTensorMap{T, S, 2, 1}, - ) where {T <: Number, S <: ElementarySpace} - f(x) = (@tensor Sx2[-1 -2; -3] := Rx[-1 -2; 1 2] * x[1 2; -3]) - x1, info = linsolve(f, Sx, x0, 0, 1) +Solve the equations `Rx x = Sx` with initial guess `x0`. +""" +function _solve_als( + Rx::AbstractTensorMap{T, S, N, N}, + Sx::GenericMPSTensor{S, N}, + x0::GenericMPSTensor{S, N}; kwargs... + ) where {T, S, N} + @assert N >= 2 + pR = (codomainind(Rx), domainind(Rx)) + pX = ((1, (3:(N + 1))...), (2,)) + pRX = ((1, N + 1, (2:(N - 1))...), (N,)) + f(x) = tensorcontract(Rx, pR, false, x, pX, false, pRX) + x1, info = linsolve(f, Sx, x0, 0, 1; kwargs...) return x1, info end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index bfd0f7f01..0572d4852 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -1,5 +1,5 @@ """ -Replace bond environment `benv` by its positive approximant `Z† Z` +Find positive approximant `Z† Z` of a norm tensor `benv` (returns the "half environment" `Z`) ``` ┌-----------------┐ ┌---------------┐ @@ -11,9 +11,9 @@ Replace bond environment `benv` by its positive approximant `Z† Z` └-----------------┘ └---------------┘ ``` """ -function positive_approx(benv::BondEnv) +function positive_approx(benv::AbstractTensorMap{T, S, N, N}) where {T, S, N} # eigen-decomposition: benv = U D U' - D, U = eigh_full((benv + benv') / 2) + D, U = eigh_full!(project_hermitian(benv)) # determine if `env` is (mostly) positive or negative sgn = sign(sum(D.data)) # When optimizing the truncation of a bond, diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl index b60818530..b7ae6f42f 100644 --- a/src/algorithms/time_evolution/apply_mpo.jl +++ b/src/algorithms/time_evolution/apply_mpo.jl @@ -112,7 +112,7 @@ Then the fidelity is just ``` =# """ -Perform QR decomposition through a PEPS tensor +Perform QR decomposition through a `GenericMPSTensor` ``` ╱ ╱ -←-R0-←-M-←- => ---Q-←-R1-←- @@ -120,19 +120,22 @@ Perform QR decomposition through a PEPS tensor ``` """ function qr_through( - R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} + R0::MPSBondTensor, M::GenericMPSTensor{S, N}; normalize::Bool = true + ) where {S, N} @assert !isdual(codomain(R0, 1)) @assert !isdual(domain(M, 1)) && !isdual(codomain(M, 1)) - @tensor A[-1 -2 -3 -4; -5] := R0[-1; 1] * M[1 -2 -3 -4; -5] + pR = (codomainind(R0), domainind(R0)) + pM = ((1,), Tuple(2:(N + 1))) + pRM = (codomainind(M), domainind(M)) + A = tensorcontract(R0, pR, false, M, pM, false, pRM) _, r = left_orth!(A; positive = true) normalize && normalize!(r, Inf) return r end # for `M` at the left end of the MPS function qr_through( - ::Nothing, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} + ::Nothing, M::GenericMPSTensor{S, N}; normalize::Bool = true + ) where {S, N} @assert !isdual(domain(M, 1)) _, r = left_orth(M; positive = true) normalize && normalize!(r, Inf) @@ -140,7 +143,7 @@ function qr_through( end """ -Perform LQ decomposition through a tensor +Perform LQ decomposition through a `GenericMPSTensor` ``` ╱ ╱ -←-L0-←-Q-←- <= -←-M-←-L1-←- @@ -148,21 +151,24 @@ Perform LQ decomposition through a tensor ``` """ function lq_through( - M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true - ) where {S <: ElementarySpace} + M::GenericMPSTensor{S, N}, L1::MPSBondTensor; normalize::Bool = true + ) where {S, N} @assert !isdual(domain(L1, 1)) @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1)) - @tensor A[-1; -2 -3 -4 -5] := M[-1 -2 -3 -4; 1] * L1[1; -5] + pM = (codomainind(M), domainind(M)) + pL = (codomainind(L1), domainind(L1)) + pML = ((1,), Tuple(2:(N + 1))) + A = tensorcontract(M, pM, false, L1, pL, false, pML) l, _ = right_orth!(A; positive = true) normalize && normalize!(l, Inf) return l end # for `M` at the right end of the MPS function lq_through( - M::GenericMPSTensor{S, 4}, ::Nothing; normalize::Bool = true - ) where {S <: ElementarySpace} + M::GenericMPSTensor{S, N}, ::Nothing; normalize::Bool = true + ) where {S, N} @assert !isdual(codomain(M, 1)) - A = permute(M, ((1,), (2, 3, 4, 5))) + A = permute(M, ((1,), Tuple(2:(N + 1))); copy = true) l, _ = right_orth!(A; positive = true) normalize && normalize!(l, Inf) return l @@ -171,7 +177,7 @@ end """ Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond """ -function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}} +function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor} # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3 N = length(Ms) # get the first R and the last L @@ -214,11 +220,14 @@ function _proj_from_RL( end """ -Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices, -find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds. +Given a cluster `Ms`, find all projectors `Pa`, `Pb` +and Schmidt weights `wts` on internal bonds. """ -function _get_allprojs(Ms, Rs, Ls, truncs::Vector{E}) where {E <: TruncationStrategy} +function _get_allprojs( + Ms::Vector{T}, truncs::Vector{E} + ) where {T <: GenericMPSTensor, E <: TruncationStrategy} N = length(Ms) + Rs, Ls = _get_allRLs(Ms) @assert length(truncs) == N - 1 projs_errs = map(1:(N - 1)) do i trunc = if isa(truncs[i], FixedSpaceTruncation) @@ -258,14 +267,16 @@ Find projectors to truncate internal bonds of the cluster `Ms`. """ function _cluster_truncate!( Ms::Vector{T}, truncs::Vector{E} - ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationStrategy} - Rs, Ls = _get_allRLs(Ms) - Pas, Pbs, wts, ϵs = _get_allprojs(Ms, Rs, Ls, truncs) + ) where {T <: GenericMPSTensor, E <: TruncationStrategy} + Pas, Pbs, wts, ϵs = _get_allprojs(Ms, truncs) # apply projectors # M1 -- (Pa1,wt1,Pb1) -- M2 -- (Pa2,wt2,Pb2) -- M3 for (i, (Pa, Pb)) in enumerate(zip(Pas, Pbs)) - @tensor (Ms[i])[-1 -2 -3 -4; -5] := (Ms[i])[-1 -2 -3 -4; 1] * Pa[1; -5] - @tensor (Ms[i + 1])[-1 -2 -3 -4; -5] := Pb[-1; 1] * (Ms[i + 1])[1 -2 -3 -4; -5] + Ms[i] = Ms[i] * twistdual(Pa, 1) + pP = ((1,), (2,)) + pM = ((1,), Tuple(2:numind(Ms[i + 1]))) + pPM = (codomainind(Ms[i + 1]), domainind(Ms[i + 1])) + Ms[i + 1] = tensorcontract(Pb, pP, false, Ms[i + 1], pM, false, pPM) end return wts, ϵs, Pas, Pbs end diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index d8feb8cb3..57fd68cc7 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -74,17 +74,11 @@ function bond_truncate( perm_ab = ((1, 3), (4, 2)) a, s0, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) a, b = absorb_s(a, s0, b) - #= temporarily reorder axes of a and b to - 1 -a/b- 2 - ↓ - 3 - =# - perm = ((1, 3), (2,)) - a, b = permute(a, perm), permute(b, perm) + # put b in MPS axis order + b = permute(b, ((1, 2), (3,))) ab = _combine_ab(a, b) # cost function will be normalized by initial value - cost00 = cost_function_als(benv, ab, a2b2) - fid = fidelity(benv, ab, a2b2) + cost00, fid = cost_function_als(benv, ab, a2b2) cost0, fid0, Δcost, Δfid, Δs = cost00, fid, NaN, NaN, NaN verbose && @info "ALS init" * _als_message(0, cost0, fid, Δcost, Δfid, Δs, 0.0) for iter in 1:(alg.maxiter) @@ -99,15 +93,14 @@ function bond_truncate( =# Ra = _tensor_Ra(benv, b) Sa = _tensor_Sa(benv, b, a2b2) - a, info_a = _solve_ab(Ra, Sa, a) + a, info_a = _solve_als(Ra, Sa, a) # Fixing `a`, solve for `b` from `Rb b = Sb` Rb = _tensor_Rb(benv, a) Sb = _tensor_Sb(benv, a, a2b2) - b, info_b = _solve_ab(Rb, Sb, b) + b, info_b = _solve_als(Rb, Sb, b) @debug "Bond truncation info" info_a info_b ab = _combine_ab(a, b) - cost = cost_function_als(benv, ab, a2b2) - fid = fidelity(benv, ab, a2b2) + cost, fid = cost_function_als(benv, ab, a2b2) # TODO: replace with truncated svdvals (without calculating u, vh) _, s, _ = svd_trunc!(permute(ab, perm_ab); trunc = alg.trunc) # fidelity, cost and normalized bond-s change diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl index ea339c8d8..fd9c685bd 100644 --- a/src/algorithms/truncation/fullenv_truncation.jl +++ b/src/algorithms/truncation/fullenv_truncation.jl @@ -49,24 +49,7 @@ between two states specified by the bond matrices `b1`, `b2` function inner_prod( benv::BondEnv{T, S}, b1::AbstractTensorMap{T, S, 1, 1}, b2::AbstractTensorMap{T, S, 1, 1} ) where {T <: Number, S <: ElementarySpace} - val = @tensor conj(b1[1; 2]) * benv[1 2; 3 4] * b2[3; 4] - return val -end - -""" -$(SIGNATURES) - -Given the bond environment `benv`, calculate the fidelity -between two states specified by the bond matrices `b1`, `b2` -``` - F(b1, b2) = (⟨b1|b2⟩ ⟨b2|b1⟩) / (⟨b1|b1⟩ ⟨b2|b2⟩) -``` -""" -function fidelity( - benv::BondEnv{T, S}, b1::AbstractTensorMap{T, S, 1, 1}, b2::AbstractTensorMap{T, S, 1, 1} - ) where {T <: Number, S <: ElementarySpace} - return abs2(inner_prod(benv, b1, b2)) / - real(inner_prod(benv, b1, b1) * inner_prod(benv, b2, b2)) + return @tensor conj(b1[1; 2]) * benv[1 2; 3 4] * b2[3; 4] end """ @@ -252,7 +235,7 @@ function fullenv_truncate( l, info_l = linsolve(Base.Fix1(*, B), p, l, 0, 1) @debug "Bond truncation info" info_l info_r @tensor b1[-1; -2] = l[-1 1] * vh[1; -2] - fid = fidelity(benv, b0, b1) + _, fid = cost_function_als(benv, b0, b1) u, s, vh = svd_trunc!(b1; trunc = alg.trunc) # determine convergence s_nrm = norm(s0, Inf) diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 58b0a4158..25f56c200 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -5,6 +5,7 @@ using TensorKit using PEPSKit using LinearAlgebra using KrylovKit +using PEPSKit: cost_function_als Random.seed!(0) maxiter = 600 @@ -19,6 +20,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') # random positive-definite environment Z = randn(Float64, Vext ← Vbond) benv = Z' * Z + normalize!(benv, Inf) # untruncated bond tensor a2b2 = randn(Float64, Vbondl ⊗ Vbondr ← Vphy' ⊗ Vphy') a2, s, b2 = svd_compact(permute(a2b2, perm_ab)) @@ -26,7 +28,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') # bond tensor (truncated SVD initialization) a0, s, b0 = svd_trunc(permute(a2b2, perm_ab); trunc = trunc) a0, b0 = PEPSKit.absorb_s(a0, s, b0) - fid0 = PEPSKit.fidelity(benv, PEPSKit._combine_ab(a0, b0), a2b2) + fid0 = cost_function_als(benv, PEPSKit._combine_ab(a0, b0), a2b2)[2] @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() for (label, alg) in ( @@ -36,7 +38,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') a1, ss[label], b1, info = PEPSKit.bond_truncate(a2, b2, benv, alg) @info "$label improved fidelity = $(info.fid)." # display(ss[label]) - @test info.fid ≈ PEPSKit.fidelity(benv, PEPSKit._combine_ab(a1, b1), a2b2) + @test info.fid ≈ cost_function_als(benv, PEPSKit._combine_ab(a1, b1), a2b2)[2] @test info.fid > fid0 end @test isapprox(ss["ALS"], ss["FET"], atol = 1.0e-3) From 2fdd026f4544b214c1de1c815dd7c9192051d904 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 6 Apr 2026 17:59:50 +0800 Subject: [PATCH 02/27] bondenv_ctm for PEPO --- .../contractions/bondenv/benv_ctm.jl | 44 ++++++----- .../contractions/bondenv/benv_tools.jl | 10 +++ test/bondenv/benv_ctm.jl | 73 ++++++++++++------- 3 files changed, 80 insertions(+), 47 deletions(-) diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 9b97c1576..45262cc67 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -1,14 +1,16 @@ """ Construct the environment (norm) tensor ``` - C1---T1---------T1---C2 r-1 - | ‖ ‖ | - T4===XX== ==YY===T2 r - | ‖ ‖ | - C4---T3---------T3---C3 r+1 - c-1 c c+1 c+2 + -1 C1---E1---------E1---C2 + | ‖ ‖ | + 0 E4===XX== ==YY===E2 + | ‖ ‖ | + 1 C4---E3---------E3---C3 + -1 0 1 2 ``` -where `XX = X' X` and `YY = Y' Y` (stacked together). +with `X` at position `[row, col]`. +`X, Y` are unitary tensors produced when finding the reduced site tensors +with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together). Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` @@ -21,7 +23,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in └---------------------┘ ``` """ -function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) +function bondenv_ctm( + row::Int, col::Int, X::T, Y::T, env::CTMRGEnv + ) where {T} Nr, Nc = size(env.corners)[[2, 3]] cm1 = _prev(col, Nc) cp1 = _next(col, Nc) @@ -32,29 +36,29 @@ function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv c2 = env.corners[2, rm1, cp2] c3 = env.corners[3, rp1, cp2] c4 = env.corners[4, rp1, cm1] - t1X, t1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1] - t2 = env.edges[2, row, cp2] - t3X, t3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1] - t4 = env.edges[4, row, cm1] + e1X, e1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1] + e2 = env.edges[2, row, cp2] + e3X, e3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1] + e4 = env.edges[4, row, cm1] #= index labels - C1--χ4--T1X---------χ6---------T1Y--χ8---C2 r-1 + C1--χ4--E1X---------χ6---------E1Y--χ8---C2 r-1 | ‖ ‖ | χ2 DNX DNY χ10 | ‖ ‖ | - T4==DWX==XX===DX== ==DY===YY==DEY==T2 r + E4==DWX==XX===DX== ==DY===YY==DEY==E2 r | ‖ ‖ | χ1 DSX DSY χ9 | ‖ ‖ | - C4--χ3--T3X---------χ5---------T3Y--χ7---C3 r+1 + C4--χ3--E3X---------χ5---------E3Y--χ7---C3 r+1 c-1 c c+1 c+2 =# + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) @autoopt @tensor benv[DX1 DY1; DX0 DY0] := - c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] * - X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] * - c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] * - Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] - + c4[χ3 χ1] * e4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * e3X[χ5 DSX0 DSX1 χ3] * + twistdual(X, 1)[pX DNX0 DX0 DSX0 DWX0] * conj(X[pX DNX1 DX1 DSX1 DWX1]) * e1X[χ4 DNX0 DNX1 χ6] * + c3[χ9 χ7] * e2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * e3Y[χ7 DSY0 DSY1 χ5] * + twistdual(Y, 1)[pY DNY0 DEY0 DSY0 DY0] * conj(Y[pY DNY1 DEY1 DSY1 DY1]) * e1Y[χ6 DNY0 DNY1 χ8] normalize!(benv, Inf) return benv end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 91ad5958e..38400b919 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,3 +1,13 @@ const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace} +const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace} +const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} +# Orthogonal tensors obtained PEPSTensor/PEPOTensor +# with one physical leg factored out by `_qr_bond` const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace} const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace} + +"Convert tensor `t` connected by the bond to be truncated to a `PEPSTensor`." +_prepare_site_tensor(t::PEPSTensor) = t +_prepare_site_tensor(t::PEPOTensor) = first(fuse_physicalspaces(t)) +_prepare_site_tensor(t::PEPSOrth) = permute(insertleftunit(t, 1), ((1,), (2, 3, 4, 5))) +_prepare_site_tensor(t::PEPOOrth) = permute(t, ((1,), (2, 3, 4, 5))) diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index c6406ba43..7676c5139 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -7,8 +7,12 @@ using Random Random.seed!(100) Nr, Nc = 2, 2 +Envspace = Vect[FermionParity ⊠ U1Irrep]( + (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 +) +ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) # create Hubbard iPEPS using simple update -function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) +function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0) H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1) peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc)) @@ -20,30 +24,45 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) return peps end -peps = get_hubbard_state() -# calculate CTMRG environment -Envspace = Vect[FermionParity ⊠ U1Irrep]( - (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 -) -ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) -env, = leading_boundary(CTMRGEnv(rand, ComplexF64, peps, Envspace), peps, ctm_alg) -for row in 1:Nr, col in 1:Nc - cp1 = PEPSKit._next(col, Nc) - A, B = peps.A[row, col], peps.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) - benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) - @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] - Z = PEPSKit.positive_approx(benv) - # verify that gauge fixing can greatly reduce - # condition number for physical state bond envs - cond1 = cond(Z' * Z) - Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) - benv2 = Z2' * Z2 - cond2 = cond(benv2) - @test 1 <= cond2 < cond1 - @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" - # verify gauge fixing is done correctly - @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] - @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] - @test half ≈ half2 +function get_hubbard_pepo(t::Float64 = 1.0, U::Float64 = 8.0) + H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) + pepo = PEPSKit.infinite_temperature_density_matrix(H) + wts = SUWeight(pepo) + alg = SimpleUpdate(; + trunc = truncerror(; atol = 1.0e-10) & truncrank(4), bipartite = false + ) + pepo, = time_evolve(pepo, H, 2.0e-3, 500, alg, wts; check_interval = 100) + normalize!.(pepo.A, Inf) + return pepo +end + +function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) + network = isa(state, InfinitePEPS) ? state : InfinitePEPS(state) + env, = leading_boundary(CTMRGEnv(rand, ComplexF64, network, Envspace), network, ctm_alg) + for row in 1:Nr, col in 1:Nc + cp1 = PEPSKit._next(col, Nc) + A, B = state.A[row, col], state.A[row, cp1] + X, a, b, Y = PEPSKit._qr_bond(A, B) + benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) + Z = PEPSKit.positive_approx(benv) + # verify that gauge fixing can greatly reduce + # condition number for physical state bond envs + cond1 = cond(Z' * Z) + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + cond2 = cond(benv2) + @test 1 <= cond2 < cond1 + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" + # verify gauge fixing is done correctly + @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @test half ≈ half2 + end + return +end + +peps = get_hubbard_peps() +pepo = get_hubbard_pepo() +for state in (peps, pepo) + test_benv_ctm(state) end From 1a13c39d8beb03713caae94b1a8d2b3441aa83c4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 6 Apr 2026 20:25:37 +0800 Subject: [PATCH 03/27] Rotation of LocalCircuit --- src/PEPSKit.jl | 1 + src/algorithms/time_evolution/trotter_gate.jl | 113 --------------- src/operators/localcircuit.jl | 137 ++++++++++++++++++ test/runtests.jl | 3 + test/types/localcircuit.jl | 21 +++ 5 files changed, 162 insertions(+), 113 deletions(-) create mode 100644 src/operators/localcircuit.jl create mode 100644 test/types/localcircuit.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index a0317a6b1..f1c537ee6 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -61,6 +61,7 @@ include("utility/symmetrization.jl") include("operators/infinitepepo.jl") include("operators/transfermatrix.jl") include("operators/localoperator.jl") +include("operators/localcircuit.jl") include("operators/lattices/squarelattice.jl") include("operators/models.jl") diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl index 84a4882b5..24248ad9c 100644 --- a/src/algorithms/time_evolution/trotter_gate.jl +++ b/src/algorithms/time_evolution/trotter_gate.jl @@ -1,116 +1,3 @@ -""" -$(TYPEDEF) - -Circuit consisting of local gates and MPOs. - -## Fields - -$(TYPEDFIELDS) -""" -struct LocalCircuit{O, S} - "lattice of physical spaces on which the gates act" - lattice::Matrix{S} - - "list of `sites => gate` pairs that make up the circuit" - gates::Vector{Pair{Vector{CartesianIndex{2}}, O}} - - LocalCircuit{O, S}(lattice::Matrix{S}) where {O, S} = - new{O, S}(lattice, Vector{Pair{Vector{CartesianIndex{2}}, O}}()) -end - -LocalCircuit{O}(lattice::Matrix{<:ElementarySpace}) where {O} = - LocalCircuit{O, eltype(lattice)}(lattice) -LocalCircuit{O}(lattice, gates::Pair...) where {O} = LocalCircuit{O}(lattice, gates) - -function LocalCircuit{O}(lattice, terms) where {O} - operator = LocalCircuit{O}(lattice) - for (inds, term) in terms - add_factor!(operator, inds, term) - end - return operator -end - -# Default to Any for eltype: needs to be abstract anyways so not that much to gain -LocalCircuit(lattice, terms) = LocalCircuit{Any}(lattice, terms) -LocalCircuit(lattice, terms::Pair...) = LocalCircuit(lattice, terms) - -add_factor!(operator::LocalCircuit, inds::Tuple, term::AbstractTensorMap) = add_factor!(operator, collect(inds), term) -add_factor!(operator::LocalCircuit, inds::Vector, term::AbstractTensorMap) = add_factor!(operator, map(CartesianIndex{2}, inds), term) -function add_factor!(operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap) - # input checks - length(inds) == numin(term) == numout(term) || throw(ArgumentError("Incompatible number of indices and tensor legs")) - for (i, ind) in enumerate(inds) - ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) - physicalspace(operator, ind_translated) == domain(term)[i] == codomain(term)[i] || - throw(SpaceMismatch("Incompatible physical spaces at $(ind).")) - end - - # permute input - if !issorted(inds) - I = sortperm(inds) - inds = inds[I] - term = permute(term, (Tuple(I), Tuple(I) .+ numout(term))) - end - - # translate coordinates - I1 = first(inds) - I1_mod = CartesianIndex(mod1.(Tuple(I1), size(operator))) - inds .-= (I1 - I1_mod) - - push!(operator.gates, inds => term) - - return operator -end -# for MPO term -# TODO: consider directly use MPSKit.FiniteMPO -function add_factor!( - operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::Vector{M} - ) where {M <: AbstractTensorMap} - # input checks - length(inds) >= 2 || throw(ArgumentError("Gate MPO must act on 2 or more sites.")) - length(inds) == length(term) || throw(ArgumentError("Incompatible number of indices and length of gate MPO.")) - allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) - for (i, (ind, t)) in enumerate(zip(inds, term)) - ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) - out_ax = (i == 1) ? 1 : 2 - in_ax = (i == 1) ? 2 : 3 - physicalspace(operator, ind_translated) == space(t, out_ax) == space(t, in_ax)' || - throw(SpaceMismatch("Incompatible physical spaces at $(ind).")) - if i >= 2 - ind_prev = inds[i - 1] - sum(Tuple(ind - ind_prev) .^ 2) == 1 || throw(ArgumentError("Two consecutive sites in `inds` must be nearest neighbours for MPO terms.")) - end - end - # for MPO term, `inds` should not be sorted - push!(operator.gates, inds => term) - return operator -end - -function checklattice(::Type{Bool}, H1::LocalCircuit, H2::LocalCircuit) - return physicalspace(H1) == physicalspace(H2) -end -function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalCircuit) - return physicalspace(peps) == physicalspace(O) -end -function checklattice(::Type{Bool}, H::LocalCircuit, peps::InfinitePEPS) - return checklattice(Bool, peps, H) -end -function checklattice(::Type{Bool}, pepo::InfinitePEPO, O::LocalCircuit) - return size(pepo, 3) == 1 && physicalspace(pepo) == physicalspace(O) -end -function checklattice(::Type{Bool}, O::LocalCircuit, pepo::InfinitePEPO) - return checklattice(Bool, pepo, O) -end - -""" - physicalspace(gates::LocalCircuit) - -Return lattice of physical spaces on which the `LocalCircuit` is defined. -""" -physicalspace(gates::LocalCircuit) = gates.lattice -physicalspace(gates::LocalCircuit, args...) = physicalspace(gates)[args...] -Base.size(gates::LocalCircuit) = size(physicalspace(gates)) - const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2} """ diff --git a/src/operators/localcircuit.jl b/src/operators/localcircuit.jl new file mode 100644 index 000000000..f60a17ad2 --- /dev/null +++ b/src/operators/localcircuit.jl @@ -0,0 +1,137 @@ +""" +$(TYPEDEF) + +Circuit consisting of local gates and MPOs. + +## Fields + +$(TYPEDFIELDS) +""" +struct LocalCircuit{O, S} + "lattice of physical spaces on which the gates act" + lattice::Matrix{S} + + "list of `sites => gate` pairs that make up the circuit" + gates::Vector{Pair{Vector{CartesianIndex{2}}, O}} + + LocalCircuit{O, S}(lattice::Matrix{S}) where {O, S} = + new{O, S}(lattice, Vector{Pair{Vector{CartesianIndex{2}}, O}}()) +end + +LocalCircuit{O}(lattice::Matrix{<:ElementarySpace}) where {O} = + LocalCircuit{O, eltype(lattice)}(lattice) +LocalCircuit{O}(lattice, gates::Pair...) where {O} = LocalCircuit{O}(lattice, gates) + +function LocalCircuit{O}(lattice, terms) where {O} + operator = LocalCircuit{O}(lattice) + for (inds, term) in terms + add_factor!(operator, inds, term) + end + return operator +end + +# Default to Any for eltype: needs to be abstract anyways so not that much to gain +LocalCircuit(lattice, terms) = LocalCircuit{Any}(lattice, terms) +LocalCircuit(lattice, terms::Pair...) = LocalCircuit(lattice, terms) + +add_factor!(operator::LocalCircuit, inds::Tuple, term::AbstractTensorMap) = add_factor!(operator, collect(inds), term) +add_factor!(operator::LocalCircuit, inds::Vector, term::AbstractTensorMap) = add_factor!(operator, map(CartesianIndex{2}, inds), term) +# for AbstractTensorMap term +function add_factor!(operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap) + # input checks + length(inds) == numin(term) == numout(term) || throw(ArgumentError("Incompatible number of indices and tensor legs")) + for (i, ind) in enumerate(inds) + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) + physicalspace(operator, ind_translated) == domain(term)[i] == codomain(term)[i] || + throw(SpaceMismatch("Incompatible physical spaces at $(ind).")) + end + # permute input + if !issorted(inds) + I = sortperm(inds) + inds = inds[I] + term = permute(term, (Tuple(I), Tuple(I) .+ numout(term))) + end + # translate coordinates + I1 = first(inds) + I1_mod = CartesianIndex(mod1.(Tuple(I1), size(operator))) + inds .-= (I1 - I1_mod) + push!(operator.gates, inds => term) + return operator +end +# for MPO term +# TODO: consider directly use MPSKit.FiniteMPO +function add_factor!( + operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::Vector{M} + ) where {M <: AbstractTensorMap} + # input checks + length(inds) >= 2 || throw(ArgumentError("Gate MPO must act on 2 or more sites.")) + length(inds) == length(term) || throw(ArgumentError("Incompatible number of indices and length of gate MPO.")) + allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + for (i, (ind, t)) in enumerate(zip(inds, term)) + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) + out_ax = (i == 1) ? 1 : 2 + in_ax = (i == 1) ? 2 : 3 + physicalspace(operator, ind_translated) == space(t, out_ax) == space(t, in_ax)' || + throw(SpaceMismatch("Incompatible physical spaces at $(ind).")) + if i >= 2 + ind_prev = inds[i - 1] + sum(Tuple(ind - ind_prev) .^ 2) == 1 || throw(ArgumentError("Two consecutive sites in `inds` must be nearest neighbours for MPO terms.")) + end + end + # for MPO term, `inds` should not be sorted + push!(operator.gates, inds => term) + return operator +end + +function checklattice(::Type{Bool}, H1::LocalCircuit, H2::LocalCircuit) + return physicalspace(H1) == physicalspace(H2) +end +function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalCircuit) + return physicalspace(peps) == physicalspace(O) +end +function checklattice(::Type{Bool}, H::LocalCircuit, peps::InfinitePEPS) + return checklattice(Bool, peps, H) +end +function checklattice(::Type{Bool}, pepo::InfinitePEPO, O::LocalCircuit) + return size(pepo, 3) == 1 && physicalspace(pepo) == physicalspace(O) +end +function checklattice(::Type{Bool}, O::LocalCircuit, pepo::InfinitePEPO) + return checklattice(Bool, pepo, O) +end + +""" + physicalspace(gates::LocalCircuit) + +Return lattice of physical spaces on which the `LocalCircuit` is defined. +""" +physicalspace(gates::LocalCircuit) = gates.lattice +physicalspace(gates::LocalCircuit, args...) = physicalspace(gates)[args...] +Base.size(gates::LocalCircuit) = size(physicalspace(gates)) + +# Equality +# ----------- + +Base.:(==)(O₁::LocalCircuit, O₂::LocalCircuit) = + physicalspace(O₁) == physicalspace(O₂) && O₁.gates == O₂.gates + +# Rotation +# ---------------------- + +function Base.rotr90(H::LocalCircuit) + Hsize = size(H) + lattice2 = rotr90(physicalspace(H)) + terms2 = (siterotr90.(inds, Ref(Hsize)) => term for (inds, term) in H.gates) + return LocalCircuit(lattice2, terms2) +end +function Base.rotl90(H::LocalCircuit) + Hsize = size(H) + lattice2 = rotl90(physicalspace(H)) + terms2 = (siterotl90.(inds, Ref(Hsize)) => term for (inds, term) in H.gates) + return LocalCircuit(lattice2, terms2) +end +function Base.rot180(H::LocalCircuit) + Hsize = size(H) + lattice2 = rot180(physicalspace(H)) + terms2 = (siterot180.(inds, Ref(Hsize)) => term for (inds, term) in H.gates) + return LocalCircuit(lattice2, terms2) +end diff --git a/test/runtests.jl b/test/runtests.jl index f59691b93..24bc65bb9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,6 +23,9 @@ end @time @safetestset "LocalOperator" begin include("types/localoperator.jl") end + @time @safetestset "LocalCircuit" begin + include("types/localcircuit.jl") + end end if GROUP == "ALL" || GROUP == "CTMRG" @time @safetestset "Gauge Fixing" begin diff --git a/test/types/localcircuit.jl b/test/types/localcircuit.jl new file mode 100644 index 000000000..16d96597c --- /dev/null +++ b/test/types/localcircuit.jl @@ -0,0 +1,21 @@ +using TensorKit +using PEPSKit +using PEPSKit: LocalCircuit +using Test + +@testset "Rotation" begin + op = LocalCircuit( + [ℂ^1 ℂ^2 ℂ^3; ℂ^4 ℂ^5 ℂ^6], + ( + ((1, 1), (1, 2)) => randn(ℂ^1, ℂ^1) ⊗ randn(ℂ^2, ℂ^2), + ((2, 1), (1, 1)) => randn(ℂ^4, ℂ^4) ⊗ randn(ℂ^1, ℂ^1), + ((1, 2), (2, 3)) => randn(ℂ^2, ℂ^2) ⊗ randn(ℂ^6, ℂ^6), + ((1, 3), (2, 2)) => randn(ℂ^3, ℂ^3) ⊗ randn(ℂ^5, ℂ^5), + )... + ) + @test rot180(rot180(op)) == op + @test rotl90(rotl90(op)) == rot180(op) == rotr90(rotr90(op)) + @test physicalspace(rotl90(op)) == rotl90(physicalspace(op)) + @test physicalspace(rotr90(op)) == rotr90(physicalspace(op)) + @test physicalspace(rot180(op)) == rot180(physicalspace(op)) +end From bfb826aca6d464ac0572441a99063a3beff351f0 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 7 Apr 2026 13:33:02 +0800 Subject: [PATCH 04/27] Use vector of sites in `nearest_neighbours` etc --- src/operators/lattices/squarelattice.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/operators/lattices/squarelattice.jl b/src/operators/lattices/squarelattice.jl index 005c0bdad..102546594 100644 --- a/src/operators/lattices/squarelattice.jl +++ b/src/operators/lattices/squarelattice.jl @@ -29,19 +29,19 @@ function vertices(lattice::InfiniteSquare) end function nearest_neighbours(lattice::InfiniteSquare) - neighbors = Tuple{CartesianIndex, CartesianIndex}[] + neighbors = Vector{CartesianIndex{2}}[] for idx in vertices(lattice) - push!(neighbors, (idx, idx + CartesianIndex(0, 1))) - push!(neighbors, (idx, idx + CartesianIndex(1, 0))) + push!(neighbors, [idx, idx + CartesianIndex(0, 1)]) + push!(neighbors, [idx, idx + CartesianIndex(1, 0)]) end return neighbors end function next_nearest_neighbours(lattice::InfiniteSquare) - neighbors = Tuple{CartesianIndex, CartesianIndex}[] + neighbors = Vector{CartesianIndex{2}}[] for idx in vertices(lattice) - push!(neighbors, (idx, idx + CartesianIndex(1, 1))) - push!(neighbors, (idx + CartesianIndex(0, 1), idx + CartesianIndex(1, 0))) + push!(neighbors, [idx, idx + CartesianIndex(1, 1)]) + push!(neighbors, [idx + CartesianIndex(0, 1), idx + CartesianIndex(1, 0)]) end return neighbors end From f6554f386ab695370510d9eadaa6b15d83e0cabb Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Apr 2026 15:22:03 +0800 Subject: [PATCH 05/27] Neighbourhood tensor update --- src/PEPSKit.jl | 7 +- .../contractions/bondenv/benv_ntu.jl | 239 ++++++++++++++++++ .../contractions/bondenv/benv_tools.jl | 201 +++++++++++++++ .../contractions/bondenv/gaugefix.jl | 31 ++- src/algorithms/time_evolution/ntupdate.jl | 222 ++++++++++++++++ .../time_evolution/ntupdate3site.jl | 89 +++++++ 6 files changed, 778 insertions(+), 11 deletions(-) create mode 100644 src/algorithms/contractions/bondenv/benv_ntu.jl create mode 100644 src/algorithms/time_evolution/ntupdate.jl create mode 100644 src/algorithms/time_evolution/ntupdate3site.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index f1c537ee6..ce16b89b5 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -89,6 +89,7 @@ include("algorithms/contractions/bp_contractions.jl") include("algorithms/contractions/bondenv/benv_tools.jl") include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") +include("algorithms/contractions/bondenv/benv_ntu.jl") include("algorithms/contractions/bondenv/benv_ctm.jl") include("algorithms/contractions/correlator/peps.jl") include("algorithms/contractions/correlator/pepo_1layer.jl") @@ -111,6 +112,8 @@ include("algorithms/time_evolution/apply_mpo.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") +include("algorithms/time_evolution/ntupdate.jl") +include("algorithms/time_evolution/ntupdate3site.jl") include("algorithms/time_evolution/gaugefix_su.jl") include("algorithms/bp/beliefpropagation.jl") @@ -143,7 +146,8 @@ export fixedpoint export absorb_weight export ALSTruncation, FullEnvTruncation -export SimpleUpdate +export NNEnv, NNpEnv, NNNEnv +export SimpleUpdate, NeighbourUpdate export TimeEvolver, timestep, time_evolve export InfiniteSquareNetwork @@ -151,6 +155,7 @@ export InfinitePartitionFunction export InfinitePEPS, InfiniteTransferPEPS export SUWeight export InfinitePEPO, InfiniteTransferPEPO +export infinite_temperature_density_matrix export BPEnv, BeliefPropagation export BPGauge, SUGauge diff --git a/src/algorithms/contractions/bondenv/benv_ntu.jl b/src/algorithms/contractions/bondenv/benv_ntu.jl new file mode 100644 index 000000000..2c88c1cfe --- /dev/null +++ b/src/algorithms/contractions/bondenv/benv_ntu.jl @@ -0,0 +1,239 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + +""" +Algorithms to construct bond environment for Neighborhood Tensor Update (NTU). +""" +abstract type NeighbourEnv end + +""" +SVD with `truncrank(1)`. +""" +function _svd_cut!(t::AbstractTensorMap) + t1, s, t2 = svd_trunc!(t; trunc = truncrank(1)) + return absorb_s(t1, s, t2) +end + +""" +Algorithm struct for "NTU-NN" bond environment. +""" +struct NNEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN" approximation. +``` + -1 ●=======● + ║ ║ + 0 ●===X== ==Y===● + ║ ║ + 1 ●=======● + -1 0 1 2 +``` +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNEnv + ) where {T, S <: InfiniteState} + neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)] + m = collect_neighbors(state, row, col, neighbors) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + # southwest half + @autoopt @tensor benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] := + cor_se(m[1, 1])[Dse1 Dse0 Ds1 Ds0] * + cor_sw(m[1, 0])[Dsw1 Dsw0 Ds1 Ds0] * + edge_w(X, hair_w(m[0, -1]))[Dnw1 Dnw0 Dw1 Dw0 Dsw1 Dsw0] + normalize!(benv_sw, Inf) + # northeast half + @autoopt @tensor benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] := + cor_nw(m[-1, 0])[Dn1 Dn0 Dnw1 Dnw0] * + cor_ne(m[-1, 1])[Dne1 Dne0 Dn1 Dn0] * + edge_e(Y, hair_e(m[0, 2]))[Dne1 Dne0 Dse1 Dse0 De1 De0] + normalize!(benv_ne, Inf) + @tensor benv[Dw1 De1; Dw0 De0] := + benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] * benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] + normalize!(benv, Inf) + return benv +end + +""" +Algorithm struct for "NTU-NN+" bond environment. +""" +struct NNpEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN+" approximation. +``` + -2 ●.......● + ║ ║ + -1 ○===●=======●===○ + ║ ║ ║ ║ + 0 ●===●===X== ==Y===●===● + ║ ║ ║ ║ + 1 ○===●=======●===○ + ║ ║ + 2 ●.......● + -2 -1 0 1 2 3 +``` +Dotted lines and ○ are splitted using SVD with `truncrank(1)`. +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNpEnv + ) where {T, S <: InfiniteState} + neighbors = [ + (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (1, 2), (0, 2), (-1, 2), + (-1, 1), (-1, 0), (0, -2), (2, 0), (2, 1), (0, 3), (-2, 1), (-2, 0), + ] + ms = collect_neighbors(state, row, col, neighbors) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + + # ---- hairs (size D^2) with a 1D auxiliary leg ---- + + @tensor top[-1 -2; -3 -4] := cor_nw(ms[-2, 0])[1 2 -1 -2] * cor_ne(ms[-2, 1])[-3 -4 1 2] + tl, tr = _svd_cut!(top) + + @tensor bot[-1 -2; -3 -4] := cor_sw(ms[2, 0])[-1 -2 1 2] * cor_se(ms[2, 1])[-3 -4 1 2] + bl, br = _svd_cut!(bot) + + nw = permute(cor_nw(ms[-1, -1]), ((3, 4), (1, 2))) + nw1, nw2 = _svd_cut!(nw) + + ne = permute(cor_ne(ms[-1, 2]), ((3, 4), (1, 2))) + ne1, ne2 = _svd_cut!(ne) + + sw = permute(cor_sw(ms[1, -1]), ((1, 2), (3, 4))) + sw1, sw2 = _svd_cut!(sw) + + se = permute(cor_se(ms[1, 2]), ((3, 4), (1, 2))) + se1, se2 = _svd_cut!(se) + + m = ms[0, -1] + @tensoropt hW[anw DXw1 DXw0 asw] := + hair_w(ms[0, -2])[Dw21 Dw20] * + nw1[Dnw11 Dnw10 anw] * sw1[Dsw11 Dsw10 asw] * + twistdual(m, 1)[p Dnw10 DXw0 Dsw10 Dw20] * + conj(m[p Dnw11 DXw1 Dsw11 Dw21]) + + m = ms[0, 2] + @tensoropt hE[ane DYe1 DYe0 ase] := + hair_e(ms[0, 3])[De21 De20] * + ne2[ane Dne21 Dne20] * se2[ase Dse21 Dse20] * + twistdual(m, 1)[p Dne20 De20 Dse20 DYe0] * + conj(m[p Dne21 De21 Dse21 DYe1]) + + # ---- corners (size D^4) with a 1D auxiliary leg ---- + + m = ms[-1, 0] + @tensoropt NW[at Dn1 Dn0 DXn1 DXn0 anw] := + tl[Dtl1 Dtl0 at] * nw2[anw Dnw21 Dnw20] * + twistdual(m, 1)[p Dtl0 Dn0 DXn0 Dnw20] * + conj(m[p Dtl1 Dn1 DXn1 Dnw21]) + + m = ms[-1, 1] + @tensoropt NE[at Dn1 Dn0 DYn1 DYn0 ane] := + tr[at Dtr1 Dtr0] * ne1[Dne11 Dne10 ane] * + twistdual(m, 1)[p Dtr0 Dne10 DYn0 Dn0] * + conj(m[p Dtr1 Dne11 DYn1 Dn1]) + + m = ms[1, 0] + @tensoropt SW[asw DXs1 DXs0 Ds1 Ds0 ab] := + bl[Dbl1 Dbl0 ab] * sw2[asw Dsw21 Dsw20] * + twistdual(m, 1)[p DXs0 Ds0 Dbl0 Dsw20] * + conj(m[p DXs1 Ds1 Dbl1 Dsw21]) + + m = ms[1, 1] + @tensoropt SE[ase DYs1 DYs0 Ds1 Ds0 ab] := + br[ab Dbr1 Dbr0] * se1[Dse11 Dse10 ase] * + twistdual(m, 1)[p DYs0 Dse10 Dbr0 Ds0] * + conj(m[p DYs1 Dse11 Dbr1 Ds1]) + + # ---- left half ---- + @tensoropt benvL[at Dn1 Dn0 Dw1 Dw0 Ds1 Ds0 ab] := + hW[anw DXw1 DXw0 asw] * + NW[at Dn1 Dn0 DXn1 DXn0 anw] * + SW[asw DXs1 DXs0 Ds1 Ds0 ab] * + twistdual(X, 1)[p DXn0 Dw0 DXs0 DXw0] * + conj(X[p DXn1 Dw1 DXs1 DXw1]) + normalize!(benvL, Inf) + + # ---- right half ---- + + @tensoropt benvR[at Dn1 Dn0 De1 De0 Ds1 Ds0 ab] := + hE[ane DYe1 DYe0 ase] * + NE[at Dn1 Dn0 DYn1 DYn0 ane] * + SE[ase DYs1 DYs0 Ds1 Ds0 ab] * + twistdual(Y, 1)[p DYn0 DYe0 DYs0 De0] * + conj(Y[p DYn1 DYe1 DYs1 De1]) + normalize!(benvR, Inf) + + # ---- the full NN+ environment ---- + + @tensor benv[Dw1, De1; Dw0, De0] := + benvL[at Dn1 Dn0 Dw1 Dw0 Ds1 Ds0 ab] * + benvR[at Dn1 Dn0 De1 De0 Ds1 Ds0 ab] + normalize!(benv, Inf) + return benv +end + +""" +Algorithm struct for "NTU-NNN" bond environment. +""" +struct NNNEnv <: NeighbourEnv end +""" +Calculates the bond environment within "NTU-NNN" approximation. +``` + -1 ●===●=======●===● + ║ ║ ║ ║ + 0 ●===X== ==Y===● + ║ ║ ║ ║ + 1 ●===●=======●===● + -1 0 1 2 +``` +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNNEnv + ) where {T, S <: InfiniteState} + neighbors = [ + (-1, -1), (0, -1), (1, -1), + (1, 0), (1, 1), (1, 2), (0, 2), + (-1, 2), (-1, 1), (-1, 0), + ] + m = collect_neighbors(state, row, col, neighbors) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + #= left half + -1 ●======●=== -1/-2 + ║ ║ + 0 ●======X=== -3/-4 + ║ ║ + ....D1.....D2... + ║ ║ + 1 ●==D3==●=== -5/-6 + -1 0 + =# + vecl = enlarge_corner_nw(cor_nw(m[-1, -1]), edge_n(m[-1, 0]), edge_w(m[0, -1]), X) + @tensor vecl[:] := + cor_sw(m[1, -1])[D11 D10 D31 D30] * + edge_s(m[1, 0])[D21 D20 -5 -6 D31 D30] * + vecl[D11 D10 D21 D20 -1 -2 -3 -4] + normalize!(vecl, Inf) + #= right half + -1 -1/-2===●==D1==● + ║ ║ + ........D2.....D3... + ║ ║ + 0 -3/-4===Y======● + ║ ║ + 1 -5/-6===●======● + 0 1 + =# + vecr = enlarge_corner_se(cor_se(m[1, 2]), edge_s(m[1, 1]), edge_e(m[0, 2]), Y) + @tensor vecr[:] := + edge_n(m[-1, 1])[D11 D10 D21 D20 -1 -2] * + cor_ne(m[-1, 2])[D31 D30 D11 D10] * + vecr[D21 D20 D31 D30 -3 -4 -5 -6] + normalize!(vecr, Inf) + # combine left and right part + @tensor benv[-1 -2; -3 -4] := vecl[1 2 -1 -3 3 4] * vecr[1 2 -2 -4 3 4] + normalize!(benv, Inf) + return benv +end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 38400b919..3c8217e06 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,3 +1,10 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace} const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace} const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} @@ -11,3 +18,197 @@ _prepare_site_tensor(t::PEPSTensor) = t _prepare_site_tensor(t::PEPOTensor) = first(fuse_physicalspaces(t)) _prepare_site_tensor(t::PEPSOrth) = permute(insertleftunit(t, 1), ((1,), (2, 3, 4, 5))) _prepare_site_tensor(t::PEPOOrth) = permute(t, ((1,), (2, 3, 4, 5))) + +""" +Extract tensors in an InfinitePEPS or 1-layer InfinitePEPO +at positions `neighbors` relative to `(row, col)` +""" +function collect_neighbors( + state::InfiniteState, row::Int, col::Int, neighbors::Vector{Tuple{Int, Int}} + ) + Nr, Nc = size(state) + return Dict( + nb => _prepare_site_tensor(state[mod1(row + nb[1], Nr), mod1(col + nb[2], Nc)]) + for nb in neighbors + ) +end + +""" + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}) + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}, hairs::Vector{H}) where {T, H <: Union{Nothing, Hair}} + +Contract the physical axes and the virtual axes of `ket` with `bra` to obtain the tensor on the boundary of the bond environment. +Virtual axes specified by `open_axs` (in ascending order) are not contracted. + +# Examples + +- West "hair" tensor (`open_axs = [EAST]`) + ``` + ╱| + ┌-----ket----- 2 + | ╱ | | + | | | | + | | | ╱ + └---|-bra----- 1 + |╱ + ``` +- Northwest corner tensor (`open_axs = [EAST, SOUTH]`, `hairs = [h, nothing]`) + ``` + ╱| + ┌-----ket----- 2 + | ╱ | h + | 4 | | + | | | + | | ╱ + └-----bra----- 1 + ╱ + 3 + ``` +- West edge tensor (`open_axs = [1, 2, 3]`) + ``` + 2 + ╱ + ┌-----ket----- 4 + | ╱ | + | 6 | 1 + | | ╱ + └-----bra----- 3 + ╱ + 5 + ``` +""" +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int} + ) + # no hair tensors to be attached to virtual legs + return benv_tensor(ket, bra, open_axs, fill(nothing, 4 - length(open_axs))) +end +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}, hairs::Vector{H} + ) where {H <: Union{Nothing, Hair}} + @assert length(hairs) == 4 - length(open_axs) + ket2, nax = copy(ket), numind(ket) + axs, open_axs2 = (2:5), open_axs .+ 1 + # contract with hair tensors + hair_axs = Tuple(ax for ax in axs if ax ∉ open_axs2) + for (h, ax) in zip(hairs, hair_axs) + if h === nothing + twistdual!(ket2, ax) + continue + end + # ensure the hair doesn't change the virtual spaces + @assert space(h, 1) == space(h, 2)' + ket_indices = collect(-1:-1:-nax) + ket_indices[ax] = 1 + ket2 = ncon([h, ket2], [[-ax, 1], ket_indices]) + end + # combine bra and ket + indexlist = [-collect(1:2:(2 * nax)), -collect(2:2:(2 * nax))] + for ax in 1:nax + if ax ∉ open_axs2 + indexlist[1][ax] = indexlist[2][ax] = ax + end + end + return ncon([bra, ket2], indexlist, [true, false]) +end + +#= Free axes of different boundary tensors +(C/E/H mean corner/edge/hair) + + H_n + | + + C_nw - - E_n - - C_ne + | | | + + | | | + H_w - E_w - - ket - - E_e - H_e + | | | + + | | | + C_sw - - E_s - - C_se + + | + H_s +=# +const open_axs_hair = Dict(:n => [SOUTH], :e => [WEST], :s => [NORTH], :w => [EAST]) +const open_axs_cor = Dict( + :nw => [EAST, SOUTH], :ne => [SOUTH, WEST], :se => [NORTH, WEST], :sw => [NORTH, EAST] +) +const open_axs_edge = Dict( + :n => [EAST, SOUTH, WEST], + :e => [NORTH, SOUTH, WEST], + :s => [NORTH, EAST, WEST], + :w => [NORTH, EAST, SOUTH], +) + +# construction of hairs +for (dir, open_axs) in open_axs_hair + fname = Symbol("hair_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h1, h2, h3) = benv_tensor(ket, ket, $open_axs, [h1, h2, h3]) + end +end + +# construction of corners +for (dir, open_axs) in open_axs_cor + fname = Symbol("cor_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h1, h2) = benv_tensor(ket, ket, $open_axs, [h1, h2]) + end +end + +# construction of edges +for (dir, open_axs) in open_axs_edge + fname = Symbol("edge_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h) = benv_tensor(ket, ket, $open_axs, [h]) + end +end + +""" +Enlarge the northwest corner +``` + ctl══ D1 ══ et ══ -5/-6 + ║ ║ + D2 D3 + ║ ║ + el ══ D4 ══ X ═══ -7/-8 + ║ ║ + -1/-2 -3/-4 +``` +""" +function enlarge_corner_nw( + ctl::AbstractTensor{E, S, 4}, + et::AbstractTensor{E, S, 6}, el::AbstractTensor{E, S, 6}, + ket::PEPSTensor, bra::PEPSTensor = ket + ) where {E, S} + return @tensoropt ctl2[:] := ctl[D11 D10 D21 D20] * + et[-5 -6 D31 D30 D11 D10] * el[D21 D20 D41 D40 -1 -2] * + conj(bra[d D31 -7 -3 D41]) * twistdual(ket, 1)[d D30 -8 -4 D40] +end + +""" +Enlarge the southeast corner +``` + -1/-2 -3/-4 + ║ ║ + -5/-6 ═════ Y ══ D1 ═══ er + ║ ║ + D2 D3 + ║ ║ + -7/-8 ═════ eb ═ D4 ══ cbr +``` +""" +function enlarge_corner_se( + cbr::AbstractTensor{E, S, 4}, + eb::AbstractTensor{E, S, 6}, er::AbstractTensor{E, S, 6}, + ket::PEPSTensor, bra::PEPSTensor = ket + ) where {E, S} + return @tensoropt cbr2[:] := cbr[D31 D30 D41 D40] * + eb[D21 D20 D41 D40 -7 -8] * er[-3 -4 D31 D30 D11 D10] * + conj(bra[d -1 D11 D21 -5]) * twistdual(ket, 1)[d -2 D10 D20 -6] +end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 0572d4852..9df4c1599 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -92,12 +92,13 @@ function fixgauge_benv( end """ -When the (half) bond environment `Z` consists of two `PEPSOrth` tensors `X`, `Y` as +When the (half) bond environment `Z` consists of +two `PEPSOrth` or `PEPOOrth` tensors `X`, `Y` as ``` - ┌---------------┐ ┌-------------------┐ - | | = | | , - └---Z-- --┘ └--Z0---X-- --Y--┘ - ↓ ↓ + ┌-----------------------┐ + | | + └---Z---(X)-- --(Y)---┘ + ↓ ``` apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: ``` @@ -106,15 +107,25 @@ apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: -4 - X - 1 - Rinv - -2 -4 - Linv - 1 - Y - -2 | | -3 -3 + + -2 -2 + | | + -5 - X - 1 - Rinv - -3 -5 - Linv - 1 - Y - -3 + | ╲ | ╲ + -4 -1 -4 -1 ``` """ function _fixgauge_benvXY( - X::PEPSOrth{T, S}, - Y::PEPSOrth{T, S}, - Linv::AbstractTensorMap{T, S, 1, 1}, - Rinv::AbstractTensorMap{T, S, 1, 1}, - ) where {T <: Number, S <: ElementarySpace} + X::PEPSOrth, Y::PEPSOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, + ) @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] return X, Y end +function _fixgauge_benvXY( + X::PEPOOrth, Y::PEPOOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, + ) + @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3] + @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5] + return X, Y +end diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl new file mode 100644 index 000000000..f08d9fbfd --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -0,0 +1,222 @@ +""" +$(TYPEDEF) + +Algorithm struct for neighbourhood tensor update (NTU) of InfinitePEPS or InfinitePEPO. + +## Fields + +$(TYPEDFIELDS) + +Reference: +- Physical Review B 104, 094411 (2021) +- Physical Review B 106, 195105 (2022) +""" +@kwdef struct NeighbourUpdate <: TimeEvolution + "Bond truncation algorithm after applying time evolution gate" + opt_alg::Union{ALSTruncation, FullEnvTruncation} = + ALSTruncation(; trunc = truncerror(; atol = 1.0e-10)) + "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" + imaginary_time::Bool = true + "Algorithm to construct NTU bond environment." + bondenv_alg::NeighbourEnv = NNEnv() + "When true, fix gauge of bond environment" + fixgauge::Bool = true + "When true, assume bipartite unit cell structure" + bipartite::Bool = false +end + +""" +Internal state of neighbourhood tensor update algorithm +""" +struct NTUState{S <: InfiniteState, N <: Number} + "number of performed iterations" + iter::Int + "evolved time" + t::N + "PEPS/PEPO" + psi::S +end + +""" + TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::NeighbourUpdate; t0::Number = 0.0, symmetrize_gates::Bool = false + ) + +Initialize a TimeEvolver with Hamiltonian `H` and neighbourhood tensor update `alg`, +starting from the initial state `psi0`. + +- The initial time is specified by `t0`. +- Use `symmetrize_gates = true` for second-order Trotter decomposition. +""" +function TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::NeighbourUpdate; t0::Number = 0.0, symmetrize_gates::Bool = false + ) + _timeevol_sanity_check(psi0, physicalspace(H), alg) + dt′ = _get_dt(psi0, dt, alg.imaginary_time) + gate = trotterize(H, dt′; symmetrize_gates) + state = NTUState(0, t0, psi0) + return TimeEvolver(alg, dt, nstep, gate, state) +end + +""" +Neighbourhood tensor update optimized for nearest neighbor gates +utilizing reduced bond tensors with the physical leg. +""" +function _ntu_iter( + state::InfiniteState, gate::NNGate, wts::SUWeight, + sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate + ) + @assert length(sites) == 2 + return _bond_truncate(state, wts, Tuple(sites), alg; gate) +end + +""" +One round of neighbourhood tensor update +""" +function ntu_iter( + state::InfiniteState, circuit::LocalCircuit, alg::NeighbourUpdate + ) + Nr, Nc, = size(state) + state2, wts = deepcopy(state), SUWeight(state) + info = (; fid = 1.0) + for (sites, gate) in circuit.gates + if length(sites) == 1 + # 1-site gate + # TODO: special treatment for bipartite state + site = sites[1] + r, c = mod1(site[1], Nr), mod1(site[2], Nc) + state2[r, c] = _apply_sitegate(state2[r, c], gate) + info′ = (; fid = 1.0) + elseif length(sites) == 2 + (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) + if alg.bipartite + length(sites) > 2 && error("Multi-site MPO gates are not compatible with bipartite states.") + r > 1 && continue + end + state2, wts, info′ = _ntu_iter(state2, gate, wts, sites, alg) + (!alg.bipartite) && continue + if d == 1 + rp1, cp1 = _next(r, Nr), _next(c, Nc) + state2[rp1, cp1] = deepcopy(state2[r, c]) + state2[rp1, c] = deepcopy(state2[r, cp1]) + wts[1, rp1, cp1] = deepcopy(wts[1, r, c]) + else + rm1, cm1 = _prev(r, Nr), _prev(c, Nc) + state2[rm1, cm1] = deepcopy(state2[r, c]) + state2[r, cm1] = deepcopy(state2[rm1, c]) + wts[2, rm1, cm1] = deepcopy(wts[2, r, c]) + end + else + # N-site MPO gate (N ≥ 2) + alg.bipartite && error("Multi-site MPO gates are not compatible with bipartite states.") + state2, wts, info′ = _ntu_iter(state2, gate, wts, sites, alg) + end + # record the worst fidelity + (info′.fid < info.fid) && (info = info′) + end + return state2, wts, info +end + +function Base.iterate(it::TimeEvolver{<:NeighbourUpdate}, state = it.state) + iter, t = state.iter, state.t + (iter == it.nstep) && return nothing + psi, wts, info = ntu_iter(state.psi, it.gate, it.alg) + iter, t = iter + 1, t + it.dt + # update internal state + it.state = NTUState(iter, t, psi) + info = (; (; t, wts)..., info...) + return (psi, info), it.state +end + +""" + timestep( + it::TimeEvolver{<:NeighbourUpdate}, psi::InfiniteState; + iter::Int = it.state.iter, t::Float64 = it.state.t + ) -> (psi, info) + +Given the TimeEvolver iterator `it`, perform one step of NTU time evolution +on the input state `psi`. + +- Using `iter` and `t` to reset the current iteration number and evolved time + respectively of the TimeEvolver `it`. +""" +function MPSKit.timestep( + it::TimeEvolver{<:NeighbourUpdate}, psi::InfiniteState; + iter::Int = it.state.iter, t::Float64 = it.state.t + ) + _timeevol_sanity_check(psi, physicalspace(it.state.psi), it.alg) + state = NTUState(iter, t, psi) + result = iterate(it, state) + if result === nothing + @warn "TimeEvolver `it` has already reached the end." + return nothing + else + return first(result) + end +end + +""" + time_evolve( + it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 500 + ) -> (psi, info) + +Perform time evolution to the end of `TimeEvolver` iterator `it`. + +- `check_interval` sets the number of iterations between outputs of information. +""" +function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 500) + time_start = time0 = time() + @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" + info0 = nothing + for (psi, info) in it + iter = it.state.iter + stop = (iter == it.nstep) + showinfo = (iter == 1) || (iter % check_interval == 0) || stop + time1 = time() + if showinfo + Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) + @info @sprintf( + "NTU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, it.state.t, Δλ, time1 - time0 + ) + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + end + info0, time0 = info, time() + end + return +end + +""" + time_evolve( + psi0::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, + dt::Number, nstep::Int, alg::NeighbourUpdate; + t0::Number = 0.0, symmetrize_gates::Bool = false, + check_interval::Int = 10 + ) -> (psi, info) + +Perform time evolution on the initial state `psi0` and initial environment `env0` +with Hamiltonian `H`, using `NeighbourUpdate` algorithm `alg`, time step `dt` for +`nstep` number of steps. + +- Convergence check for ground state search is not supported. +- Set `symmetrize_gates = true` for second-order Trotter decomposition. +- Use `t0` to specify the initial time of `psi0`. +- `check_interval` sets the interval to output information (and check convergence). + Output during the evolution can be turned off by setting `check_interval <= 0`. +- `info` is a NamedTuple containing information of the evolution, + including the time `info.t` evolved since `psi0`. +""" +function MPSKit.time_evolve( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::NeighbourUpdate; symmetrize_gates::Bool = false, + t0::Number = 0.0, check_interval::Int = 10 + ) + it = TimeEvolver(psi0, H, dt, nstep, alg; t0, symmetrize_gates) + return time_evolve(it; check_interval) +end diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl new file mode 100644 index 000000000..00055f54d --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -0,0 +1,89 @@ +""" +Neighbourhood tensor update with N-site MPO `gate` (N ≥ 2). +""" +function _ntu_iter( + state::InfiniteState, gate::Vector{T}, wts::SUWeight, + sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate + ) where {T <: AbstractTensorMap} + Nr, Nc = size(state) + + # apply gate MPO without truncation + Ms, _, invperms = _get_cluster(state, sites) + flips = [isdual(space(M, 1)) for M in Ms[2:end]] + _flip_virtuals!(Ms, flips) # flip virtual arrows in `Ms` to ← + _apply_gatempo!(Ms, gate) + _flip_virtuals!(Ms, flips) # restore virtual arrows in `Ms` + state2, wts2 = deepcopy(state), deepcopy(wts) + for (M, s, invperm) in zip(Ms, sites, invperms) + s′ = CartesianIndex(mod1(s[1], Nr), mod1(s[2], Nc)) + state2[s′] = permute(M, invperm) + end + + # truncate each bond sequentially along the path + info = (; fid = 1.0) + for bondsites in zip(sites, Iterators.drop(sites, 1)) + state2, wts2, info′ = _bond_truncate(state2, wts2, bondsites, alg) + # record the worst fidelity + (info′.fid < info.fid) && (info = info′) + end + return state2, wts2, info +end + +""" +Truncate a nearest neighbor bond between `site1` and `site2` +after rotating the bond to standard x direction `A ← B`. +""" +function _bond_truncate( + state::InfiniteState, wts::SUWeight, + (site1, site2)::NTuple{2, CartesianIndex{2}}, + alg::NeighbourUpdate; gate::Union{NNGate, Nothing} = nothing + ) + # rotate bond to standard x direction `A ← B` + ucell = size(state)[1:2] + bond, rev = _nn_bondrev(site1, site2, ucell) + state2 = _bond_rotation(state, bond[1], rev; inv = false) + wts2 = _bond_rotation(wts, bond[1], rev; inv = false) + + # rotated bond tensors + siteA = if bond[1] == 1 + rev ? siterot180(site1, ucell) : site1 + else + rev ? siterotl90(site1, ucell) : siterotr90(site1, ucell) + end + row, col = mod1.(Tuple(siteA), size(state2)[1:2]) + cp1 = _next(col, size(state2, 2)) + A, B = state2[row, col], state2[row, cp1] + + # create bond environment + X, a, b, Y = _qr_bond(A, B; trunc = trunctol(; rtol = 1.0e-12)) + benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg) + @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" + if alg.fixgauge + Z = positive_approx(benv) + Z, a, b, (Linv, Rinv) = fixgauge_benv(Z, a, b) + X, Y = _fixgauge_benvXY(X, Y, Linv, Rinv) + benv = Z' * Z + @debug "cond(L) = $(LinearAlgebra.cond(Linv)); cond(R): $(LinearAlgebra.cond(Rinv))" + @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))" + end + + # (optional) apply the NN gate + if !(gate === nothing) + a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) + else + a = permute(a, ((1, 2), (3,))) + b = permute(b, ((1,), (2, 3))) + end + + a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) + A, B = _qr_bond_undo(X, a, b, Y) + normalize!(A, Inf) + normalize!(B, Inf) + normalize!(s, Inf) + state2[row, col], state2[row, cp1], wts2[1, row, col] = A, B, s + + # rotate back tensors and bond weight + state2 = _bond_rotation(state2, bond[1], rev; inv = true) + wts2 = _bond_rotation(wts2, bond[1], rev; inv = true) + return state2, wts2, info +end From 2540e48b7716d150a26156fe304b8000b33babc9 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Apr 2026 16:07:59 +0800 Subject: [PATCH 06/27] Fix benv_tensor for dual physical space --- src/algorithms/contractions/bondenv/benv_tools.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 3c8217e06..c3ae13ded 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -109,6 +109,7 @@ function benv_tensor( indexlist[1][ax] = indexlist[2][ax] = ax end end + twistdual!(ket2, 1) return ncon([bra, ket2], indexlist, [true, false]) end From 24145b1d3ab3dee85a299bb0f0eeaf1700714865 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Apr 2026 16:29:45 +0800 Subject: [PATCH 07/27] Add NTU tests --- test/bondenv/benv_ntu.jl | 72 ++++++++++++++++++++++++ test/runtests.jl | 6 ++ test/timeevol/j1j2_finiteT.jl | 68 +++++++++++----------- test/timeevol/tf_ising_realtime.jl | 90 ++++++++++++++++++++++++++++++ 4 files changed, 201 insertions(+), 35 deletions(-) create mode 100644 test/bondenv/benv_ntu.jl create mode 100644 test/timeevol/tf_ising_realtime.jl diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl new file mode 100644 index 000000000..d3e2e4ffd --- /dev/null +++ b/test/bondenv/benv_ntu.jl @@ -0,0 +1,72 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random + +Nr, Nc = 2, 2 +Random.seed!(20) + +function test_ntu_env( + state::Union{InfinitePEPS, InfinitePEPO}, env_alg::Alg + ) where {Alg <: PEPSKit.NeighbourEnv} + @info "Testing $(typeof(env_alg))" + for row in 1:Nr, col in 1:Nc + cp1 = PEPSKit._next(col, Nc) + A, B = state.A[row, col], state.A[row, cp1] + X, a, b, Y = PEPSKit._qr_bond(A, B) + @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY] + benv = PEPSKit.bondenv_ntu(row, col, X, Y, state, env_alg) + # this is a result of `_qr_bond` + @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] + # NTU bond environments are exact and should be positive definite + @test benv' ≈ benv + benv = project_hermitian(benv) + nrm1 = PEPSKit.inner_prod(benv, ab, ab) + @test nrm1 ≈ real(nrm1) + D, U = eigh_full(benv) + @test minimum(D.data) >= 0 + # gauge fixing + Z = PEPSKit.sdiag_pow(D, 0.5) * U' + @assert benv ≈ Z' * Z + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + # gauge fixing should reduce condition number + cond = LinearAlgebra.cond(benv) + cond2 = LinearAlgebra.cond(benv2) + @test cond2 <= cond + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)" + # verify gauge transformation of X, Y + @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY] + nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2) + X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv) + benv3 = PEPSKit.bondenv_ntu(row, col, X2, Y2, state, env_alg) + benv3 *= norm(benv2, Inf) + nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2) + @test benv2 ≈ benv3 + @test nrm1 ≈ nrm2 ≈ nrm3 + end + return +end + +@testset "NTU bond environment ($(sym))" for sym in [U1Irrep, FermionParity] + Pspace = Vect[sym](0 => 1, 1 => 1) + V2 = Vect[sym](0 => 4, 1 => 1) + V3 = Vect[sym](0 => 3, 1 => 2) + V4 = Vect[sym](0 => 2, 1 => 2) + V5 = Vect[sym](0 => 2, 1 => 3) + Pspaces = fill(Pspace, (Nr, Nc)) + Nspaces = [V2 V2; V4 V4] + Espaces = [V3 V5; V5 V3] + + for state in [ + InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces), + InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces), + ] + normalize!.(state.A, Inf) + for env_alg in (NNEnv(), NNpEnv(), NNNEnv()) + test_ntu_env(state, env_alg) + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 24bc65bb9..66f8ac952 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -90,6 +90,9 @@ end @time @safetestset "Gauge fixing" begin include("bondenv/benv_gaugefix.jl") end + @time @safetestset "Exact NTU bond environments" begin + include("bondenv/benv_ntu.jl") + end @time @safetestset "Full update bond environment" begin include("bondenv/benv_ctm.jl") end @@ -110,6 +113,9 @@ end @time @safetestset "J1-J2 model at finite temperature" begin include("timeevol/j1j2_finiteT.jl") end + @time @safetestset "Transverse Field Ising model real-time evolution" begin + include("timeevol/tf_ising_realtime.jl") + end end if GROUP == "ALL" || GROUP == "TOOLBOX" @time @safetestset "Density matrix from double-layer PEPO" begin diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl index 63576e6a1..1217333ec 100644 --- a/test/timeevol/j1j2_finiteT.jl +++ b/test/timeevol/j1j2_finiteT.jl @@ -5,9 +5,8 @@ import MPSKitModels: σˣ, σᶻ using PEPSKit # Benchmark energy from high-temperature expansion -# at β = 0.3, 0.6 -# Physical Review B 86, 045139 (2012) Fig. 15-16 -bm = [-0.1235, -0.213] +const βs = [0.2, 0.4, 0.6] +const bm = [-0.08624893, -0.15688984, -0.21300888] function converge_env(state, χ::Int) trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12) @@ -23,37 +22,36 @@ ham = j1_j2_model( ) pepo0 = PEPSKit.infinite_temperature_density_matrix(ham) wts0 = SUWeight(pepo0) -# 7 = 1 (spin-0) + 2 x 3 (spin-1) -trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12) -check_interval = 100 -dt, nstep = 1.0e-3, 600 +dt, nstep, check_interval = 5.0e-3, 40, 40 -# PEPO approach -alg = SimpleUpdate(; trunc = trunc_pepo, purified = false) -evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0) -pepo, wts, info = time_evolve(evolver; check_interval) -env = converge_env(InfinitePartitionFunction(pepo), 16) -energy = expectation_value(pepo, ham, env) / (Nr * Nc) -@info "β = $(dt * nstep): tr(ρH) = $(energy)" -@test dt * nstep ≈ info.t -@test energy ≈ bm[2] atol = 5.0e-3 - -# PEPS (purified PEPO) approach -alg = SimpleUpdate(; trunc = trunc_pepo, purified = true) -evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0) -pepo, wts, info = time_evolve(evolver; check_interval) -env = converge_env(InfinitePartitionFunction(pepo), 16) -energy = expectation_value(pepo, ham, env) / (Nr * Nc) -@info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)" -@test energy ≈ bm[1] atol = 5.0e-3 - -# test BP gauge fixing for purified iPEPO -bp_alg = BeliefPropagation(; maxiter = 100, tol = 1.0e-9) -bp_env, = leading_boundary(BPEnv(ones, Float64, pepo), pepo, bp_alg) -pepo, = gauge_fix(pepo, BPGauge(), bp_env) +@testset "Simple update" begin + # 7 = 1 (spin-0) + 2 x 3 (spin-1) + trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12) + alg = SimpleUpdate(; trunc = trunc_pepo, purified = true) + pepo, wts = deepcopy(pepo0), deepcopy(wts0) + for (β, bme) in zip(βs, bm) + t0 = β - βs[1] + pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0, check_interval) + # measure energy + env = converge_env(InfinitePEPS(pepo), 16) + energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) + @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)" + @test energy ≈ bme atol = 5.0e-3 + end +end -env = converge_env(InfinitePEPS(pepo), 16) -energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) -@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)" -@test dt * nstep ≈ info.t -@test energy ≈ bm[2] atol = 5.0e-3 +@testset "Neighbourhood tensor update" begin + trunc_pepo = truncrank(4) & truncerror(; atol = 1.0e-12) + opt_alg = ALSTruncation(; trunc = trunc_pepo, tol = 1.0e-10) + alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv()) + pepo = deepcopy(pepo0) + for (β, bme) in zip(βs, bm) + t0 = β - βs[1] + pepo, info = time_evolve(pepo, ham, dt, nstep, alg; t0, check_interval) + # measure energy + env = converge_env(InfinitePEPS(pepo), 16) + energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) + @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)" + @test energy ≈ bme atol = 2.0e-2 + end +end diff --git a/test/timeevol/tf_ising_realtime.jl b/test/timeevol/tf_ising_realtime.jl new file mode 100644 index 000000000..977d5f283 --- /dev/null +++ b/test/timeevol/tf_ising_realtime.jl @@ -0,0 +1,90 @@ +using Test +using TensorKit +import MPSKitModels: S_zz, σˣ +using PEPSKit +using Printf +using Accessors: @set + +const hc = 3.044382 +const formatter = Printf.Format("t = %.2f, ⟨σˣ⟩ = %.7e + %.7e im.") +# real time evolution of ⟨σx⟩ +# benchmark data from Physical Review B 104, 094411 (2021) +# Figure 6(a) calculated with D = 8 and χ = 32 +const data = [ + # 0.01 9.9920027e-1 + 0.06 9.7274912e-1 + 0.11 9.1973182e-1 + 0.16 8.6230618e-1 + 0.21 8.1894325e-1 + 0.26 8.0003708e-1 + 0.31 8.0081082e-1 + 0.36 8.0979257e-1 + # 0.41 8.1559623e-1 + # 0.46 8.1541661e-1 + # 0.51 8.1274128e-1 +] + +# the fully polarized state +peps0 = InfinitePEPS(zeros, ComplexF64, ℂ^2, ℂ^1; unitcell = (2, 2)) +for t in peps0.A + t[1, 1, 1, 1, 1] = 1.0 + t[2, 1, 1, 1, 1] = 1.0 +end +lattice = collect(space(t, 1) for t in peps0.A) + +# Hamiltonian +op = LocalOperator(lattice, ((1, 1),) => σˣ()) +ham = transverse_field_ising(ComplexF64, Trivial, InfiniteSquare(2, 2); J = 1.0, g = hc) + +# truncation strategy +Dcut, chi = 4, 16 +trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dcut) +trunc_env = truncerror(; atol = 1.0e-10) & truncrank(chi) + +ctm_alg = SequentialCTMRG(; + tol = 1.0e-8, maxiter = 50, verbosity = 2, + trunc = trunc_env, projector_alg = :fullinfinite +) + +interval = 5 +ntu_alg = NeighbourUpdate(; + opt_alg = FullEnvTruncation(; trunc = trunc_peps, tol = 1.0e-10), + bondenv_alg = NNEnv(), imaginary_time = false +) + +# do one step of NTU to match benchmark data +peps0, = time_evolve(peps0, ham, 0.01, 6, ntu_alg) +@info "Space of `peps0[1, 1]` = $(space(peps0[1, 1]))." +env0 = CTMRGEnv(ones, ComplexF64, peps0, ℂ^1) +env0, = leading_boundary(env0, peps0, ctm_alg) +# measure magnetization +magx = expectation_value(peps0, op, env0) +@info Printf.format(formatter, 0.06, real(magx), imag(magx)) +@test isapprox(magx, data[1, 2]; atol = 0.005) + +@testset "Neigborhood tensor update" begin + peps, env = deepcopy(peps0), deepcopy(env0) + count = 2 + evolver = TimeEvolver(peps, ham, 0.01, 30, ntu_alg; t0 = 0.06) + spaces0 = collect(space(t) for t in peps.A) + for (peps, info) in evolver + !(evolver.state.iter % interval == 0) && continue + spaces = collect(space(t) for t in peps.A) + if spaces0 == spaces + env, = leading_boundary(env, peps, ctm_alg) + else + env = complex(CTMRGEnv(info.wts)) + env, = leading_boundary(env, peps, ctm_alg) + end + # monitor the growth of env dimension + corner = env.corners[1, 1, 1] + corner_dim = dim.(space(corner, ax) for ax in 1:numind(corner)) + @info "Dimension of env.corner[1, 1, 1] = $(corner_dim)." + # measure magnetization + magx = expectation_value(peps, op, env) + @info Printf.format(formatter, info.t, real(magx), imag(magx)) + @test isapprox(magx, data[count, 2]; atol = 0.005) + count += 1 + spaces0 = spaces + end +end From f9d029ae9ddbde0fa7940893b1624d4b36ad2837 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 8 Apr 2026 23:03:45 +0800 Subject: [PATCH 08/27] Remove unneeded bipartite restrictions --- src/algorithms/time_evolution/ntupdate.jl | 5 +---- src/algorithms/time_evolution/simpleupdate.jl | 5 +---- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index f08d9fbfd..b8b90b00f 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -91,10 +91,7 @@ function ntu_iter( info′ = (; fid = 1.0) elseif length(sites) == 2 (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) - if alg.bipartite - length(sites) > 2 && error("Multi-site MPO gates are not compatible with bipartite states.") - r > 1 && continue - end + alg.bipartite && r > 1 && continue state2, wts, info′ = _ntu_iter(state2, gate, wts, sites, alg) (!alg.bipartite) && continue if d == 1 diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 9d6842844..a28f86d8b 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -142,10 +142,7 @@ function su_iter( state2[r, c] = _apply_sitegate(state2[r, c], gate; alg.purified) elseif length(sites) == 2 (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) - if alg.bipartite - length(sites) > 2 && error("Multi-site MPO gates are not compatible with bipartite states.") - r > 1 && continue - end + alg.bipartite && r > 1 && continue ϵ′ = _su_iter!(state2, gate, env2, sites, alg) ϵ = max(ϵ, ϵ′) (!alg.bipartite) && continue From c51b54b6e5028404aebe16ffb96abb6597428809 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 12 Apr 2026 10:08:05 +0800 Subject: [PATCH 09/27] Test `timestep` for NTU --- test/timeevol/timestep.jl | 36 ++++++++++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl index 39153a898..b9c2ea9fb 100644 --- a/test/timeevol/timestep.jl +++ b/test/timeevol/timestep.jl @@ -3,14 +3,16 @@ using Random using TensorKit using PEPSKit +Nr, Nc = 2, 2 +H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1)) +Pspace, Vspace = ℂ^2, ℂ^4 +ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc)) +dt, nstep = 0.1, 20 +trunc = truncerror(; atol = 1.0e-10) & truncrank(4) + @testset "SimpleUpdate timestep" begin - Nr, Nc = 2, 2 - H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1)) - Pspace, Vspace = ℂ^2, ℂ^4 - ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc)) + alg = SimpleUpdate(; trunc) env0 = SUWeight(ψ0) - alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4)) - dt, nstep = 1.0e-2, 50 # manual timestep evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0) ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing @@ -31,3 +33,25 @@ using PEPSKit @test env1 == env2 == env3 @test info1 == info2 == info3 end + +@testset "NeighbourUpdate timestep" begin + alg = NeighbourUpdate(; opt_alg = ALSTruncation(; trunc)) + # manual timestep + evolver = TimeEvolver(ψ0, H, dt, nstep, alg) + ψ1, info1 = deepcopy(ψ0), nothing + for iter in 0:(nstep - 1) + ψ1, info1 = timestep(evolver, ψ1) + end + # time_evolve + ψ2, info2 = time_evolve(ψ0, H, dt, nstep, alg) + # for-loop syntax + ## manually reset internal state of evolver + evolver.state = PEPSKit.NTUState(0, 0.0, ψ0) + ψ3, info3 = nothing, nothing, nothing + for state in evolver + ψ3, info3 = state + end + # results should be *exactly* the same + @test ψ1 == ψ2 == ψ3 + @test info1 == info2 == info3 +end From e8f07eb723d0bb2d078c73e5164ba7f44df57328 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 12 Apr 2026 11:59:59 +0800 Subject: [PATCH 10/27] Convergence check for NTU --- src/PEPSKit.jl | 2 +- src/algorithms/time_evolution/ntupdate.jl | 104 +++++++++++++++++++--- 2 files changed, 91 insertions(+), 15 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index ce16b89b5..4700303a5 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -3,7 +3,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Random using Compat -using Accessors: @set, @reset +using Accessors: @set, @reset, @insert using VectorInterface import VectorInterface as VI diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index b8b90b00f..98fced498 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -79,7 +79,7 @@ function ntu_iter( state::InfiniteState, circuit::LocalCircuit, alg::NeighbourUpdate ) Nr, Nc, = size(state) - state2, wts = deepcopy(state), SUWeight(state) + state2, wts = copy(state), SUWeight(state) info = (; fid = 1.0) for (sites, gate) in circuit.gates if length(sites) == 1 @@ -96,14 +96,14 @@ function ntu_iter( (!alg.bipartite) && continue if d == 1 rp1, cp1 = _next(r, Nr), _next(c, Nc) - state2[rp1, cp1] = deepcopy(state2[r, c]) - state2[rp1, c] = deepcopy(state2[r, cp1]) - wts[1, rp1, cp1] = deepcopy(wts[1, r, c]) + state2[rp1, cp1] = copy(state2[r, c]) + state2[rp1, c] = copy(state2[r, cp1]) + wts[1, rp1, cp1] = copy(wts[1, r, c]) else rm1, cm1 = _prev(r, Nr), _prev(c, Nc) - state2[rm1, cm1] = deepcopy(state2[r, c]) - state2[r, cm1] = deepcopy(state2[rm1, c]) - wts[2, rm1, cm1] = deepcopy(wts[2, r, c]) + state2[rm1, cm1] = copy(state2[r, c]) + state2[r, cm1] = copy(state2[rm1, c]) + wts[2, rm1, cm1] = copy(wts[2, r, c]) end else # N-site MPO gate (N ≥ 2) @@ -156,14 +156,91 @@ end """ time_evolve( - it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 500 + it::TimeEvolver{<:NeighbourUpdate}, + [H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm]; + tol::Float64 = 1.0e-7, check_interval::Int = 10 ) -> (psi, info) -Perform time evolution to the end of `TimeEvolver` iterator `it`. +Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`, +or until convergence of energy set by a positive `tol`. -- `check_interval` sets the number of iterations between outputs of information. +To enable convergence check (for imaginary time evolution of InfinitePEPS only), +provide the Hamiltonian `H`, CTMRG environment `env`, CTMRG algorithm `ctm_alg` +and setting `tol > 0`. + +`check_interval` sets the number of iterations between energy checks +(for ground state search) and outputs of information. """ -function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 500) +function MPSKit.time_evolve( + it::TimeEvolver{<:NeighbourUpdate}, + H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm; + tol::Float64 = 1.0e-7, check_interval::Int = 10 + ) + @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" + time_start = time0 = time() + psi0 = copy(it.state.psi) + @assert (psi0 isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + # initial energy + env, = leading_boundary(env, psi0, ctm_alg) + energy = expectation_value(psi0, H, env) / prod(size(psi0)) + @info @sprintf("NTU iter 0: E = %.4e", energy) + info0 = (; energy, env) + # start evolving + energy0, ΔE = energy, 0.0 + iter0, t0 = it.state.iter, it.state.t + for (psi, info) in it + iter = it.state.iter + showinfo = (iter == 1) || (iter % check_interval == 0) || (iter == it.nstep) + !showinfo && continue + # bond weight change + Δλ = hasproperty(info0, :wts) ? compare_weights(info.wts, info0.wts) : NaN + # reconverge environment + if all(space(t) == space(t0) for (t, t0) in zip(psi.A, psi0.A)) + # recreate `env` from bond weights if psi virtual space changed + env = CTMRGEnv(info.wts) + end + env, = leading_boundary(env, psi, ctm_alg) + # measure energy + energy = expectation_value(psi, H, env) / prod(size(psi)) + ΔE = energy - energy0 + info = @insert info.energy = energy + info = @insert info.env = env + # show information + time1 = time() + @info @sprintf( + "NTU iter %-6d: E = %.5f, ΔE = %.3e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, energy, ΔE, Δλ, time1 - time0 + ) + # determine whether to stop evolution + stop = false + if (ΔE <= 0 && abs(ΔE) < tol) + stop = true + @info "NTU: energy has converged." + end + if ΔE > 0 + stop = true + @warn "NTU: energy has increased. Abort evolution and return results from last check." + psi, info, energy = psi0, info0, energy0 + it.state = NTUState(iter0, t0, psi0) + end + if iter == it.nstep + stop = true + @info "NTU: reached maximum iteration." + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + else + iter0, t0 = it.state.iter, it.state.t + psi0, energy0, info0 = psi, energy, info + end + time0 = time() + end + return +end + +function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 50) time_start = time0 = time() @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" info0 = nothing @@ -197,11 +274,10 @@ end check_interval::Int = 10 ) -> (psi, info) -Perform time evolution on the initial state `psi0` and initial environment `env0` -with Hamiltonian `H`, using `NeighbourUpdate` algorithm `alg`, time step `dt` for +Perform time evolution on the initial state `psi0` with Hamiltonian `H`, +using `NeighbourUpdate` algorithm `alg`, time step `dt` for `nstep` number of steps. -- Convergence check for ground state search is not supported. - Set `symmetrize_gates = true` for second-order Trotter decomposition. - Use `t0` to specify the initial time of `psi0`. - `check_interval` sets the interval to output information (and check convergence). From a4de5ff54a339d5e4cc612beb77403959fbf3a7a Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 15 Apr 2026 13:30:47 +0800 Subject: [PATCH 11/27] Streamline convergence checks --- src/algorithms/time_evolution/ntupdate.jl | 110 ++++++++-------------- 1 file changed, 39 insertions(+), 71 deletions(-) diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index 98fced498..cda69b436 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -155,34 +155,56 @@ function MPSKit.timestep( end """ - time_evolve( - it::TimeEvolver{<:NeighbourUpdate}, - [H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm]; - tol::Float64 = 1.0e-7, check_interval::Int = 10 - ) -> (psi, info) +$(SIGNATURES) -Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`, -or until convergence of energy set by a positive `tol`. +Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`. +`check_interval` sets the number of iterations between outputs of information. +""" +function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 50) + time_start = time0 = time() + @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" + info0 = nothing + for (psi, info) in it + iter = it.state.iter + stop = (iter == it.nstep) + showinfo = (iter == 1) || (iter % check_interval == 0) || stop + time1 = time() + if showinfo + Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) + @info @sprintf( + "NTU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, it.state.t, Δλ, time1 - time0 + ) + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + end + info0, time0 = info, time() + end + return +end -To enable convergence check (for imaginary time evolution of InfinitePEPS only), -provide the Hamiltonian `H`, CTMRG environment `env`, CTMRG algorithm `ctm_alg` -and setting `tol > 0`. +""" +$(SIGNATURES) -`check_interval` sets the number of iterations between energy checks -(for ground state search) and outputs of information. +Perform time evolution until convergence of `TimeEvolver` iterator `it`. +For `NeighbourUpdate`, convergence is determined by the change of energy +between two checks being smaller than `tol`. """ function MPSKit.time_evolve( - it::TimeEvolver{<:NeighbourUpdate}, + it::TimeEvolver{<:NeighbourUpdate, G, S}, H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm; tol::Float64 = 1.0e-7, check_interval::Int = 10 - ) + ) where {G, S <: NTUState{<:InfinitePEPS}} @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" time_start = time0 = time() psi0 = copy(it.state.psi) - @assert (psi0 isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + @assert it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." # initial energy env, = leading_boundary(env, psi0, ctm_alg) - energy = expectation_value(psi0, H, env) / prod(size(psi0)) + energy = real(expectation_value(psi0, H, env)) / prod(size(psi0)) @info @sprintf("NTU iter 0: E = %.4e", energy) info0 = (; energy, env) # start evolving @@ -201,7 +223,7 @@ function MPSKit.time_evolve( end env, = leading_boundary(env, psi, ctm_alg) # measure energy - energy = expectation_value(psi, H, env) / prod(size(psi)) + energy = real(expectation_value(psi, H, env)) / prod(size(psi)) ΔE = energy - energy0 info = @insert info.energy = energy info = @insert info.env = env @@ -239,57 +261,3 @@ function MPSKit.time_evolve( end return end - -function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 50) - time_start = time0 = time() - @info "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" - info0 = nothing - for (psi, info) in it - iter = it.state.iter - stop = (iter == it.nstep) - showinfo = (iter == 1) || (iter % check_interval == 0) || stop - time1 = time() - if showinfo - Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) - @info @sprintf( - "NTU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", - it.state.iter, it.state.t, Δλ, time1 - time0 - ) - end - if stop - time_end = time() - @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) - return psi, info - end - info0, time0 = info, time() - end - return -end - -""" - time_evolve( - psi0::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, - dt::Number, nstep::Int, alg::NeighbourUpdate; - t0::Number = 0.0, symmetrize_gates::Bool = false, - check_interval::Int = 10 - ) -> (psi, info) - -Perform time evolution on the initial state `psi0` with Hamiltonian `H`, -using `NeighbourUpdate` algorithm `alg`, time step `dt` for -`nstep` number of steps. - -- Set `symmetrize_gates = true` for second-order Trotter decomposition. -- Use `t0` to specify the initial time of `psi0`. -- `check_interval` sets the interval to output information (and check convergence). - Output during the evolution can be turned off by setting `check_interval <= 0`. -- `info` is a NamedTuple containing information of the evolution, - including the time `info.t` evolved since `psi0`. -""" -function MPSKit.time_evolve( - psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, - alg::NeighbourUpdate; symmetrize_gates::Bool = false, - t0::Number = 0.0, check_interval::Int = 10 - ) - it = TimeEvolver(psi0, H, dt, nstep, alg; t0, symmetrize_gates) - return time_evolve(it; check_interval) -end From c14fa161ee54c4c84f82aafda343075b85268c71 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 17 Apr 2026 20:09:33 +0800 Subject: [PATCH 12/27] Merge #360 --- src/algorithms/time_evolution/apply_gate.jl | 5 -- src/algorithms/time_evolution/apply_mpo.jl | 10 +-- src/algorithms/time_evolution/simpleupdate.jl | 79 +++++++++++++------ .../time_evolution/simpleupdate3site.jl | 9 ++- test/bondenv/benv_ctm.jl | 2 +- test/ctmrg/suweight.jl | 7 +- test/runtests.jl | 3 + test/timeevol/cluster_projectors.jl | 3 +- test/timeevol/fixedspacetruncation.jl | 32 ++++++++ 9 files changed, 109 insertions(+), 41 deletions(-) create mode 100644 test/timeevol/fixedspacetruncation.jl diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl index 363b23f8d..139b68196 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -134,11 +134,6 @@ function _apply_gate( else @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4] end - trunc = if trunc isa FixedSpaceTruncation - need_flip ? truncspace(flip(V)) : truncspace(V) - else - trunc - end a, s, b, ϵ = svd_trunc!(a2b2; trunc, alg = LAPACK_QRIteration()) a, b = absorb_s(a, s, b) if need_flip diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl index b7ae6f42f..60bfc6db7 100644 --- a/src/algorithms/time_evolution/apply_mpo.jl +++ b/src/algorithms/time_evolution/apply_mpo.jl @@ -229,14 +229,8 @@ function _get_allprojs( N = length(Ms) Rs, Ls = _get_allRLs(Ms) @assert length(truncs) == N - 1 - projs_errs = map(1:(N - 1)) do i - trunc = if isa(truncs[i], FixedSpaceTruncation) - tspace = space(Ms[i + 1], 1) - isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace) - else - truncs[i] - end - return _proj_from_RL(Rs[i], Ls[i]; trunc) + projs_errs = map(zip(Rs, Ls, truncs)) do (R, L, trunc) + return _proj_from_RL(R, L; trunc) end Pas = map(Base.Fix2(getindex, 1), projs_errs) wts = map(Base.Fix2(getindex, 2), projs_errs) diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index a28f86d8b..3bcfdc256 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -89,19 +89,23 @@ function _su_iter!( sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) Nr, Nc = size(state) - truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) - @assert length(sites) == 2 && length(truncs) == 1 + trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc))) + @assert length(sites) == 2 Ms, open_vaxs, = _get_cluster(state, sites, env; permute = false) normalize!.(Ms, Inf) # rotate bond, rev = _nn_bondrev(sites..., (Nr, Nc)) A, B = _bond_rotation.(Ms, bond[1], rev; inv = false) + if trunc isa FixedSpaceTruncation + V = west_virtualspace(B) + trunc = truncspace(isdual(V) ? flip(V) : V) + end # apply gate ϵ, s = 0.0, nothing gate_axs = alg.purified ? (1:1) : (1:2) for gate_ax in gate_axs X, a, b, Y = _qr_bond(A, B; gate_ax, positive = true) - a, s, b, ϵ′ = _apply_gate(a, b, gate, truncs[1]) + a, s, b, ϵ′ = _apply_gate(a, b, gate, trunc) ϵ = max(ϵ, ϵ′) A, B = _qr_bond_undo(X, a, b, Y) end @@ -202,10 +206,8 @@ function MPSKit.timestep( end """ - time_evolve( - it::TimeEvolver{<:SimpleUpdate}; - tol::Float64 = 0.0, check_interval::Int = 500 - ) -> (psi, env, info) + time_evolve(it; check_interval = 500) -> (psi, env, info) + time_evolve(it, H; tol = 1.0e-8, check_interval = 500) -> (psi, env, info) Perform time evolution to the end of `TimeEvolver` iterator `it`, or until convergence of `SUWeight` set by a positive `tol`. @@ -215,15 +217,41 @@ or until convergence of `SUWeight` set by a positive `tol`. - `check_interval` sets the number of iterations between outputs of information. """ function MPSKit.time_evolve( - it::TimeEvolver{<:SimpleUpdate}; - tol::Float64 = 0.0, check_interval::Int = 500 + it::TimeEvolver{<:SimpleUpdate}; check_interval::Int = 500 ) time_start = time() - check_convergence = (tol > 0) @info "--- Time evolution (simple update), dt = $(it.dt) ---" - if check_convergence - @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + env0, time0 = it.state.env, time() + for (psi, env, info) in it + iter = it.state.iter + diff = compare_weights(env0, env) + stop = (iter == it.nstep) + showinfo = (check_interval > 0) && + ((iter % check_interval == 0) || (iter == 1) || stop) + time1 = time() + if showinfo + @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))" + @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0) + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, env, info + else + env0 = env + end + time0 = time() end + return +end + +function MPSKit.time_evolve( + it::TimeEvolver{<:SimpleUpdate, G, S}, H::LocalOperator; + tol::Float64 = 1.0e-8, check_interval::Int = 500 + ) where {G, S <: SUState{<:InfinitePEPS}} + time_start = time() + @info "--- Time evolution (simple update), dt = $(it.dt) ---" + @assert it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." env0, time0 = it.state.env, time() for (psi, env, info) in it iter = it.state.iter @@ -233,16 +261,20 @@ function MPSKit.time_evolve( ((iter % check_interval == 0) || (iter == 1) || stop) time1 = time() if showinfo + # TODO: convert to BPEnv instead + ctmenv = CTMRGEnv(env) + energy = real(expectation_value(psi, H, ctmenv)) / prod(size(psi)) @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))" - @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0) + @info @sprintf( + "SU iter %-7d: E ≈ %.5f, |Δλ| = %.3e. Time = %.3f s/it", + iter, energy, diff, time1 - time0 + ) end - if check_convergence - if (iter == it.nstep) && (diff >= tol) - @warn "SU: bond weights have not converged." - end - if diff < tol - @info "SU: bond weights have converged." - end + if (iter == it.nstep) && (diff >= tol) + @warn "SU: bond weights have not converged." + end + if diff < tol + @info "SU: bond weights have converged." end if stop time_end = time() @@ -269,7 +301,6 @@ algorithm `alg`, time step `dt` for `nstep` number of steps. - Set `symmetrize_gates = true` for second-order Trotter decomposition. - Set `tol > 0` to enable convergence check (for imaginary time evolution of iPEPS only). - For other usages it should not be changed. - Use `t0` to specify the initial time of the evolution. - `check_interval` sets the interval to output information. Output during the evolution can be turned off by setting `check_interval <= 0`. - `info` is a NamedTuple containing information of the evolution, @@ -281,5 +312,9 @@ function MPSKit.time_evolve( tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500 ) it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0, symmetrize_gates) - return time_evolve(it; tol, check_interval) + return if tol == 0 + time_evolve(it; check_interval) + else + time_evolve(it, H; tol, check_interval) + end end diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index 4b091c830..3b16820f8 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -167,13 +167,20 @@ function _su_iter!( sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) where {T <: AbstractTensorMap} Nr, Nc = size(state) - truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) Ms, open_vaxs, invperms = _get_cluster(state, sites, env) flips = [isdual(space(M, 1)) for M in Ms[2:end]] Vphys = [codomain(M, 2) for M in Ms] normalize!.(Ms, Inf) # flip virtual arrows in `Ms` to ← _flip_virtuals!(Ms, flips) + truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) + truncs = map(enumerate(truncs)) do (i, trunc) + return if trunc isa FixedSpaceTruncation + truncspace(space(Ms[i + 1], 1)) + else + trunc + end + end # apply gate MPOs and truncate gate_axs = alg.purified ? (1:1) : (1:2) wts, ϵs = nothing, nothing diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index 7676c5139..95ece1658 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -19,7 +19,7 @@ function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0) wts = SUWeight(peps) alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4)) evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts) - peps, = time_evolve(evolver; tol = 1.0e-8, check_interval = 2000) + peps, = time_evolve(evolver, H; tol = 1.0e-8, check_interval = 2000) normalize!.(peps.A, Inf) return peps end diff --git a/test/ctmrg/suweight.jl b/test/ctmrg/suweight.jl index 3e90e28cc..12ff64e09 100644 --- a/test/ctmrg/suweight.jl +++ b/test/ctmrg/suweight.jl @@ -2,7 +2,7 @@ using Test using Random using TensorKit using PEPSKit -using PEPSKit: str, twistdual, _prev, _next +using PEPSKit: str, twistdual, _prev, _next, unitcell Vps = Dict( Z2Irrep => Vect[Z2Irrep](0 => 1, 1 => 2), @@ -43,12 +43,13 @@ end normalize!.(wts.data, Inf) end env = CTMRGEnv(wts) - for (r, c) in Tuple.(CartesianIndices(peps.A)) + for idx in CartesianIndices(unitcell(peps)) + r, c = Tuple(idx) ρ1 = su_rdm_1x1(r, c, peps, wts) if init == :trivial @test ρ1 ≈ su_rdm_1x1(r, c, peps, nothing) end - ρ2 = reduced_densitymatrix(((r, c),), peps, env) + ρ2 = reduced_densitymatrix([idx], peps, env) @test ρ1 ≈ ρ2 end end diff --git a/test/runtests.jl b/test/runtests.jl index e6627e258..bdc2f7801 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -110,6 +110,9 @@ end @time @safetestset "Time evolution with site-dependent truncation" begin include("timeevol/sitedep_truncation.jl") end + @time @safetestset "Time evolution with FixedSpaceTruncation" begin + include("timeevol/fixedspacetruncation.jl") + end @time @safetestset "Transverse field Ising model at finite temperature" begin include("timeevol/tf_ising_finiteT.jl") end diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl index 53fad0a68..1276ad2f8 100644 --- a/test/timeevol/cluster_projectors.jl +++ b/test/timeevol/cluster_projectors.jl @@ -35,7 +35,8 @@ Vspaces = [ flips = [isdual(space(M, 1)) for M in Ms1[2:end]] # no truncation Ms2 = _flip_virtuals!(deepcopy(Ms1), flips) - wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1)) + truncs = [truncrank(dim(space(M, 1))) for M in Ms2[2:end]] + wts2, ϵs, = _cluster_truncate!(Ms2, truncs) @test all((ϵ == 0) for ϵ in ϵs) normalize!.(Ms2, Inf) @test fidelity_cluster(Ms1, Ms2) ≈ 1.0 diff --git a/test/timeevol/fixedspacetruncation.jl b/test/timeevol/fixedspacetruncation.jl new file mode 100644 index 000000000..c77703fd7 --- /dev/null +++ b/test/timeevol/fixedspacetruncation.jl @@ -0,0 +1,32 @@ +using Test +using Random +using TensorKit +using PEPSKit + +elt = Float64 +ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(2, 2); J1 = 1.0, J2 = 0.5, sublattice = false) +Vphy = physicalspace(ham) +Vvir = U1Space(0 => 1, -1 / 2 => 1, 1 / 2 => 1) +Vns = [ + U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 1, -1 => 2); + U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 1, -1 => 2) +] +Ves1 = [ + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 2, -1 / 2 => 1, 3 / 2 => 1); + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1) +] +Ves2 = fill(U1Space(0 => 1, 1 => 1, -1 => 2), (2, 2)) +Venv = U1Space(0 => 2, 1 => 1, -1 => 1) +states = [ + InfinitePEPS(randn, elt, Vphy, Vns, Ves1), + InfinitePEPO(randn, elt, Vphy, Vns, Ves2), +] + +@testset "Simple update on $(typeof(state0).name.wrapper)" for state0 in states + alg = SimpleUpdate(; trunc = FixedSpaceTruncation()) + wts0 = SUWeight(state0) + state, wts, = time_evolve(state0, ham, 0.1, 1, alg, wts0) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) + end +end From 3a915c41f76770979e814d877c68cdca7e039f77 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 17 Apr 2026 20:21:26 +0800 Subject: [PATCH 13/27] FixedSpaceTruncation for NTU --- src/algorithms/time_evolution/ntupdate.jl | 24 +++++++----- .../time_evolution/ntupdate3site.jl | 38 +++++++++++++------ src/algorithms/time_evolution/simpleupdate.jl | 7 ++++ test/timeevol/fixedspacetruncation.jl | 10 +++++ 4 files changed, 57 insertions(+), 22 deletions(-) diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index cda69b436..ace6efee3 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -155,10 +155,21 @@ function MPSKit.timestep( end """ -$(SIGNATURES) + time_evolve( + it::TimeEvolver{<:NeighbourUpdate}, + [H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm]; + tol::Float64 = 1.0e-7, check_interval::Int = 10 + ) -> (psi, info) + +Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`, +or until convergence of energy set by a positive `tol`. + +To enable convergence check (for imaginary time evolution of InfinitePEPS only), +provide the Hamiltonian `H`, CTMRG environment `env`, CTMRG algorithm `ctm_alg` +and setting `tol > 0`. -Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`. -`check_interval` sets the number of iterations between outputs of information. +`check_interval` sets the number of iterations between energy checks +(for ground state search) and outputs of information. """ function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval::Int = 50) time_start = time0 = time() @@ -186,13 +197,6 @@ function MPSKit.time_evolve(it::TimeEvolver{<:NeighbourUpdate}; check_interval:: return end -""" -$(SIGNATURES) - -Perform time evolution until convergence of `TimeEvolver` iterator `it`. -For `NeighbourUpdate`, convergence is determined by the change of energy -between two checks being smaller than `tol`. -""" function MPSKit.time_evolve( it::TimeEvolver{<:NeighbourUpdate, G, S}, H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm; diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl index 00055f54d..92cc98081 100644 --- a/src/algorithms/time_evolution/ntupdate3site.jl +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -6,27 +6,37 @@ function _ntu_iter( sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate ) where {T <: AbstractTensorMap} Nr, Nc = size(state) + state, wts = copy(state), deepcopy(wts) - # apply gate MPO without truncation Ms, _, invperms = _get_cluster(state, sites) flips = [isdual(space(M, 1)) for M in Ms[2:end]] _flip_virtuals!(Ms, flips) # flip virtual arrows in `Ms` to ← + truncs = _get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc)) + truncs = map(enumerate(truncs)) do (i, trunc) + return if trunc isa FixedSpaceTruncation + truncspace(space(Ms[i + 1], 1)) + else + trunc + end + end + + # apply gate MPO without truncation _apply_gatempo!(Ms, gate) _flip_virtuals!(Ms, flips) # restore virtual arrows in `Ms` - state2, wts2 = deepcopy(state), deepcopy(wts) for (M, s, invperm) in zip(Ms, sites, invperms) s′ = CartesianIndex(mod1(s[1], Nr), mod1(s[2], Nc)) - state2[s′] = permute(M, invperm) + state[s′] = permute(M, invperm) end # truncate each bond sequentially along the path info = (; fid = 1.0) - for bondsites in zip(sites, Iterators.drop(sites, 1)) - state2, wts2, info′ = _bond_truncate(state2, wts2, bondsites, alg) + for (bondsites, trunc) in zip(zip(sites, Iterators.drop(sites, 1)), truncs) + alg′ = (@set alg.opt_alg.trunc = trunc) + state, wts, info′ = _bond_truncate(state, wts, bondsites, alg′) # record the worst fidelity (info′.fid < info.fid) && (info = info′) end - return state2, wts2, info + return state, wts, info end """ @@ -45,11 +55,7 @@ function _bond_truncate( wts2 = _bond_rotation(wts, bond[1], rev; inv = false) # rotated bond tensors - siteA = if bond[1] == 1 - rev ? siterot180(site1, ucell) : site1 - else - rev ? siterotl90(site1, ucell) : siterotr90(site1, ucell) - end + siteA = _bond_rotation(site1, bond[1], rev, ucell) row, col = mod1.(Tuple(siteA), size(state2)[1:2]) cp1 = _next(col, size(state2, 2)) A, B = state2[row, col], state2[row, cp1] @@ -68,14 +74,22 @@ function _bond_truncate( end # (optional) apply the NN gate + opt_alg = alg.opt_alg if !(gate === nothing) + trunc = if alg.opt_alg.trunc isa FixedSpaceTruncation + V = space(b, 1) + truncspace(isdual(V) ? flip(V) : V) + else + alg.opt_alg.trunc + end + @reset opt_alg.trunc = trunc a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) else a = permute(a, ((1, 2), (3,))) b = permute(b, ((1,), (2, 3))) end - a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) + a, s, b, info = bond_truncate(a, b, benv, opt_alg) A, B = _qr_bond_undo(X, a, b, Y) normalize!(A, Inf) normalize!(B, Inf) diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 3bcfdc256..bac4a3fc6 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -79,6 +79,13 @@ function _bond_rotation(x, bonddir::Int, rev::Bool; inv::Bool = false) error("`bonddir` must be 1 (for x-bonds) or 2 (for y-bonds).") end end +function _bond_rotation(x::CartesianIndex{2}, bonddir::Int, rev::Bool, unitcell::NTuple{2, Int}) + return if bonddir == 1 + rev ? siterot180(x, unitcell) : x + else + rev ? siterotl90(x, unitcell) : siterotr90(x, unitcell) + end +end """ Simple update optimized for nearest neighbor gates diff --git a/test/timeevol/fixedspacetruncation.jl b/test/timeevol/fixedspacetruncation.jl index c77703fd7..1777c0436 100644 --- a/test/timeevol/fixedspacetruncation.jl +++ b/test/timeevol/fixedspacetruncation.jl @@ -30,3 +30,13 @@ states = [ @test space(t) == space(t0) end end + +@testset "Neighborhood tensor update on $(typeof(state0).name.wrapper)" for state0 in states + opt_alg = ALSTruncation(; trunc = FixedSpaceTruncation()) + alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv()) + evolver = TimeEvolver(state0, ham, 0.1, 1, alg) + state, = time_evolve(evolver) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) + end +end From 6bfb34f4cf64934e52e49240c5e1674cb9f594e4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 17 Apr 2026 20:21:49 +0800 Subject: [PATCH 14/27] Fix time_evolve interface in tests --- test/timeevol/j1j2_finiteT.jl | 3 ++- test/timeevol/timestep.jl | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl index 1217333ec..1e3845f39 100644 --- a/test/timeevol/j1j2_finiteT.jl +++ b/test/timeevol/j1j2_finiteT.jl @@ -47,7 +47,8 @@ end pepo = deepcopy(pepo0) for (β, bme) in zip(βs, bm) t0 = β - βs[1] - pepo, info = time_evolve(pepo, ham, dt, nstep, alg; t0, check_interval) + evolver = TimeEvolver(pepo, ham, dt, nstep, alg; t0) + pepo, info = time_evolve(evolver; check_interval) # measure energy env = converge_env(InfinitePEPS(pepo), 16) energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl index b9c2ea9fb..ce86a4569 100644 --- a/test/timeevol/timestep.jl +++ b/test/timeevol/timestep.jl @@ -43,7 +43,8 @@ end ψ1, info1 = timestep(evolver, ψ1) end # time_evolve - ψ2, info2 = time_evolve(ψ0, H, dt, nstep, alg) + evolver = TimeEvolver(ψ0, H, dt, nstep, alg) + ψ2, info2 = time_evolve(evolver) # for-loop syntax ## manually reset internal state of evolver evolver.state = PEPSKit.NTUState(0, 0.0, ψ0) From 26344d209b1ca30435510327d752dcb387326f9a Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 18 Apr 2026 08:08:05 +0800 Subject: [PATCH 15/27] Fix realtime Ising test --- test/timeevol/tf_ising_realtime.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/timeevol/tf_ising_realtime.jl b/test/timeevol/tf_ising_realtime.jl index 977d5f283..b34318dcf 100644 --- a/test/timeevol/tf_ising_realtime.jl +++ b/test/timeevol/tf_ising_realtime.jl @@ -53,7 +53,7 @@ ntu_alg = NeighbourUpdate(; ) # do one step of NTU to match benchmark data -peps0, = time_evolve(peps0, ham, 0.01, 6, ntu_alg) +peps0, = time_evolve(TimeEvolver(peps0, ham, 0.01, 6, ntu_alg)) @info "Space of `peps0[1, 1]` = $(space(peps0[1, 1]))." env0 = CTMRGEnv(ones, ComplexF64, peps0, ℂ^1) env0, = leading_boundary(env0, peps0, ctm_alg) From 8ed7d13a9a0651a989a65bca11010984dc48370b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Apr 2026 15:44:07 +0800 Subject: [PATCH 16/27] Fix FixedSpaceTruncation for simple update --- src/PEPSKit.jl | 2 +- src/algorithms/bp/beliefpropagation.jl | 2 +- src/algorithms/bp/gaugefix.jl | 12 +- src/algorithms/time_evolution/apply_gate.jl | 6 +- src/algorithms/time_evolution/gaugefix_su.jl | 10 +- src/algorithms/time_evolution/simpleupdate.jl | 28 +++-- .../time_evolution/simpleupdate3site.jl | 22 +--- src/algorithms/time_evolution/time_evolve.jl | 36 +++--- src/algorithms/time_evolution/trotter_gate.jl | 2 - .../truncation/truncationschemes.jl | 56 +++++----- src/environments/bp_environments.jl | 9 ++ src/environments/suweight.jl | 9 ++ src/operators/infinitepepo.jl | 8 ++ test/bp/gaugefix.jl | 34 +++++- test/runtests.jl | 2 +- test/timeevol/sitedep_truncation.jl | 103 +++++++----------- 16 files changed, 183 insertions(+), 158 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4700303a5..44d3d42d8 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -106,9 +106,9 @@ include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_truncation.jl") -include("algorithms/time_evolution/trotter_gate.jl") include("algorithms/time_evolution/apply_gate.jl") include("algorithms/time_evolution/apply_mpo.jl") +include("algorithms/time_evolution/trotter_gate.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl index 0246d7e20..0aa900912 100644 --- a/src/algorithms/bp/beliefpropagation.jl +++ b/src/algorithms/bp/beliefpropagation.jl @@ -59,7 +59,7 @@ function leading_boundary(env₀::BPEnv, network::InfiniteSquareNetwork, alg::Be end function leading_boundary(env₀::BPEnv, state::InfiniteState, alg::BeliefPropagation) if alg.bipartite - @assert _state_bipartite_check(state) + _is_bipartite(state) || error("Input state is not bipartite") end return leading_boundary(env₀, InfiniteSquareNetwork(state), alg) end diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl index 561dd2183..5c98dc265 100644 --- a/src/algorithms/bp/gaugefix.jl +++ b/src/algorithms/bp/gaugefix.jl @@ -7,16 +7,6 @@ Algorithm for gauging PEPS with belief propagation fixed point messages. # TODO: add options end -function _bpenv_bipartite_check(env::BPEnv) - for (r, c) in Iterators.product(1:2, 1:2) - r′, c′ = _next(r, 2), _next(c, 2) - if !all(env[:, r, c] .== env[:, r′, c′]) - return false - end - end - return true -end - """ gauge_fix(psi::Union{InfinitePEPS, InfinitePEPO}, alg::BPGauge, env::BPEnv) @@ -25,7 +15,7 @@ an [`InfinitePEPO`](@ref) interpreted as purified state with two physical legs) using fixed point environment `env` of belief propagation. """ function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv) - bipartite = _state_bipartite_check(psi) && _bpenv_bipartite_check(env) + bipartite = _is_bipartite(psi) && _is_bipartite(env) psi′ = copy(psi) XXinv = map(eachcoordinate(psi, 1:2)) do I _, X, Xinv = _bp_gauge_fix!(CartesianIndex(I), psi′, env) diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl index 139b68196..5ffb98fee 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -1,3 +1,5 @@ +const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2} + """ Apply 1-site `gate` on the PEPS or PEPO tensor `a`. """ @@ -125,8 +127,8 @@ Apply 2-site `gate` on the reduced matrices `a`, `b` """ function _apply_gate( a::AbstractTensorMap, b::AbstractTensorMap, - gate::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy - ) where {T <: Number, S <: ElementarySpace} + gate::NNGate, trunc::TruncationStrategy + ) V = space(b, 1) need_flip = isdual(V) if isdual(space(a, 2)) diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl index 324a8d604..13a2782d6 100644 --- a/src/algorithms/time_evolution/gaugefix_su.jl +++ b/src/algorithms/time_evolution/gaugefix_su.jl @@ -37,22 +37,24 @@ end Fix the gauge of `psi` using trivial simple update. """ function gauge_fix(psi::InfiniteState, alg::SUGauge) + time0 = time() gates = _trivial_gates(scalartype(psi), physicalspace(psi)) - su_alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite = _state_bipartite_check(psi)) + trunc = _get_fixedspacetrunc(psi) + su_alg = SimpleUpdate(; trunc, bipartite = _is_bipartite(psi)) wts0 = SUWeight(psi) # use default constructor to avoid calculation of exp(-H * 0) evolver = TimeEvolver(su_alg, 0.0, alg.maxiter, gates, SUState(0, 0.0, psi, wts0)) for (i, (psi′, wts, info)) in enumerate(evolver) ϵ = compare_weights(wts, wts0) if i >= alg.miniter && ϵ < alg.tol - @info "Trivial SU conv $i: |Δλ| = $ϵ." + @info "Trivial SU conv $i: |Δλ| = $ϵ, time = $(time() - time0) s" return psi′, wts, ϵ end if i == alg.maxiter - @warn "Trivial SU cancel $i: |Δλ| = $ϵ." + @warn "Trivial SU cancel $i: |Δλ| = $ϵ, time = $(time() - time0) s" return psi′, wts, ϵ end - wts0 = deepcopy(wts) + wts0 = wts end return end diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index bac4a3fc6..ecebda1c5 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -7,19 +7,19 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO. $(TYPEDFIELDS) """ -@kwdef struct SimpleUpdate <: TimeEvolution +@kwdef struct SimpleUpdate{T <: TruncationStrategy} <: TimeEvolution "Truncation strategy for bonds updated by Trotter gates" - trunc::TruncationStrategy + trunc::T "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" imaginary_time::Bool = true "When true, force decomposition of nearest neighbor gates to MPOs." force_mpo::Bool = false "When true, assume bipartite unit cell structure" bipartite::Bool = false - "(Only applicable to InfinitePEPO) + "(Only applicable to InfinitePEPO) When true, the PEPO is regarded as a purified PEPS, and updated as `|ρ(t + dt)⟩ = exp(-H dt/2) |ρ(t)⟩`. - When false, the PEPO is updated as + When false, the PEPO is updated as `ρ(t + dt) = exp(-H dt/2) ρ(t) exp(-H dt/2)`." purified::Bool = true end @@ -62,6 +62,12 @@ function TimeEvolver( # create Trotter gates gate = trotterize(H, dt′; symmetrize_gates, force_mpo = alg.force_mpo) state = SUState(0, t0, psi0, env0) + # convert FixedSpaceTruncation to site-dependent `truncspace`s + if alg.trunc isa FixedSpaceTruncation + trunc = _get_fixedspacetrunc(psi0) + @reset alg.trunc = trunc + end + # TODO: bipartite check for alg.trunc after equality is defined for all kinds of truncation strategies # TODO: check gates for bipartite case return TimeEvolver(alg, dt, nstep, gate, state) end @@ -96,8 +102,8 @@ function _su_iter!( sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) Nr, Nc = size(state) - trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc))) @assert length(sites) == 2 + trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc))) Ms, open_vaxs, = _get_cluster(state, sites, env; permute = false) normalize!.(Ms, Inf) # rotate @@ -159,14 +165,14 @@ function su_iter( (!alg.bipartite) && continue if d == 1 rp1, cp1 = _next(r, Nr), _next(c, Nc) - state2[rp1, cp1] = deepcopy(state2[r, c]) - state2[rp1, c] = deepcopy(state2[r, cp1]) - env2[1, rp1, cp1] = deepcopy(env2[1, r, c]) + state2[rp1, cp1] = copy(state2[r, c]) + state2[rp1, c] = copy(state2[r, cp1]) + env2[1, rp1, cp1] = copy(env2[1, r, c]) else rm1, cm1 = _prev(r, Nr), _prev(c, Nc) - state2[rm1, cm1] = deepcopy(state2[r, c]) - state2[r, cm1] = deepcopy(state2[rm1, c]) - env2[2, rm1, cm1] = deepcopy(env2[2, r, c]) + state2[rm1, cm1] = copy(state2[r, c]) + state2[r, cm1] = copy(state2[rm1, c]) + env2[2, rm1, cm1] = copy(env2[2, r, c]) end else # N-site MPO gate (N ≥ 2) diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index 3b16820f8..df0d767ab 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -224,24 +224,14 @@ updated by the Trotter evolution MPO. """ function _get_cluster_trunc( trunc::TruncationStrategy, sites::Vector{CartesianIndex{2}}, - (Nrow, Ncol)::NTuple{2, Int} + unitcell::NTuple{2, Int} ) return map(zip(sites, Iterators.drop(sites, 1))) do (site1, site2) - diff = site2 - site1 - if diff == CartesianIndex(0, 1) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) - return truncation_strategy(trunc, 1, r, c) - elseif diff == CartesianIndex(0, -1) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) - return truncation_strategy(trunc, 1, r, c) - elseif diff == CartesianIndex(1, 0) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) - return truncation_strategy(trunc, 2, r, c) - elseif diff == CartesianIndex(-1, 0) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) - return truncation_strategy(trunc, 2, r, c) - else - error("The path `sites` contains a long-range bond.") + (d, r, c), rev = _nn_bondrev(site1, site2, unitcell) + t = truncation_strategy(trunc, d, r, c) + if rev && isa(t, TruncationSpace) + t = truncspace(flip(t.space)') end + return t end end diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl index ee5f236c1..dd51c2f54 100644 --- a/src/algorithms/time_evolution/time_evolve.jl +++ b/src/algorithms/time_evolution/time_evolve.jl @@ -30,22 +30,6 @@ end Base.iterate(it::TimeEvolver) = iterate(it, it.state) -function _state_bipartite_check(psi::InfiniteState) - if isa(psi, InfinitePEPO) - @assert size(psi, 3) == 1 "Input InfinitePEPO is expected to have only one layer." - end - if !(size(psi, 1) == size(psi, 2) == 2) - return false - end - for (r, c) in Iterators.product(1:2, 1:2) - r′, c′ = _next(r, 2), _next(c, 2) - if psi[r, c] != psi[r′, c′] - return false - end - end - return true -end - """ Process the Trotter time step `dt` according to the intended usage. """ @@ -79,7 +63,7 @@ function _timeevol_sanity_check( @assert ψ₀ isa InfinitePEPO "alg.purified = false is only applicable to PEPO." end if hasfield(typeof(alg), :bipartite) && alg.bipartite - @assert _state_bipartite_check(ψ₀) "Input state is not bipartite with 2 x 2 unit cell." + @assert _is_bipartite(ψ₀) "Input state is not bipartite with 2 x 2 unit cell." end return nothing end @@ -94,3 +78,21 @@ function MPSKit.infinite_temperature_density_matrix(H::LocalOperator) end return InfinitePEPO(cat(A; dims = 3)) end + +""" +Get the `SiteDependentTruncation` used by time evolution +that preserves virtual spaces of `state`. +""" +function _get_fixedspacetrunc(state::InfiniteState) + if state isa InfinitePEPO + size(state, 3) != 1 && error("Input InfinitePEPO is expect to have only one layer.") + end + Nr, Nc = size(state) + return SiteDependentTruncation( + map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c) + V = domain(state[r, c], (d == 1) ? EAST : NORTH) + isdual(V) && (V = flip(V)) + return truncspace(V) + end + ) +end diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl index 24248ad9c..67ae916d3 100644 --- a/src/algorithms/time_evolution/trotter_gate.jl +++ b/src/algorithms/time_evolution/trotter_gate.jl @@ -1,5 +1,3 @@ -const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2} - """ Convert an N-site gate (N ≥ 2) to MPO by SVD, in which the axes are ordered as diff --git a/src/algorithms/truncation/truncationschemes.jl b/src/algorithms/truncation/truncationschemes.jl index 71b5e9dbe..228f9065d 100644 --- a/src/algorithms/truncation/truncationschemes.jl +++ b/src/algorithms/truncation/truncationschemes.jl @@ -1,16 +1,41 @@ """ $(TYPEDEF) -CTMRG specific truncation strategy for `svd_trunc` which keeps the bond space on which the SVD -is performed fixed. Since different environment directions and unit cell entries might -have different spaces, this truncation style is different from `TruncationSpace`. +SVD truncation strategy which preserves the `CTMRGEnv` environment virtual spaces, +or `InfinitePEPS`, `InfinitePEPO` virtual spaces. """ struct FixedSpaceTruncation <: TruncationStrategy end +""" +$(TYPEDEF) + +SVD truncation strategy specified for each nearest neighbor bond +in `InfinitePEPS`, `InfinitePEPO`: +- `trunc[1, r, c]` applies to the x-bond between `[r, c]` and `[r, c+1]`. + If it is a `TruncationSpace`, the space refers to the east domain of + `[r, c]` or its `flip`, whichever is non-dual. +- `trunc[2, r, c]` applies to the y-bond between `[r, c]` and `[r-1, c]`. + If it is a `TruncationSpace`, the space refers to the north domain of + `[r, c]` or its `flip`, whichever is non-dual. +""" struct SiteDependentTruncation{T <: TruncationStrategy} <: TruncationStrategy truncs::Array{T, 3} + + function SiteDependentTruncation(truncs::Array{T, 3}) where {T} + # TODO: generalize it to CTMRGEnv + size(truncs, 1) != 2 && throw( + DimensionMismatch( + "The first dimension of `truncs` must have a size of 2. Got $(size(truncs, 1))." + ) + ) + return new{T}(truncs) + end end +Base.getindex(trunc::SiteDependentTruncation, args...) = Base.getindex(trunc.truncs, args...) + +# TODO: _is_bipartite(trunc::SiteDependentTruncation) + const TRUNCATION_STRATEGY_SYMBOLS = IdDict{Symbol, Type{<:TruncationStrategy}}( :notrunc => MatrixAlgebraKit.NoTruncation, :truncerror => MatrixAlgebraKit.TruncationByError, @@ -40,28 +65,5 @@ end function truncation_strategy( trunc::SiteDependentTruncation, direction::Int, row::Int, col::Int ) - return trunc.truncs[direction, row, col] -end - -# TODO: type piracy -Base.rotl90(trunc::TruncationStrategy) = trunc - -function Base.rotl90(trunc::SiteDependentTruncation) - directions, rows, cols = size(trunc.truncs) - truncs_rotated = similar(trunc.truncs, directions, cols, rows) - - if directions == 2 - truncs_rotated[NORTH, :, :] = circshift( - rotl90(trunc.truncs[EAST, :, :]), (0, -1) - ) - truncs_rotated[EAST, :, :] = rotl90(trunc.truncs[NORTH, :, :]) - elseif directions == 4 - for dir in 1:4 - dir′ = _prev(dir, 4) - truncs_rotated[dir′, :, :] = rotl90(trunc.truncs[dir, :, :]) - end - else - throw(ArgumentError("Unsupported number of directions for rotl90: $directions")) - end - return SiteDependentTruncation(truncs_rotated) + return trunc[direction, row, col] end diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl index 765d31c0c..ea7ae9301 100644 --- a/src/environments/bp_environments.jl +++ b/src/environments/bp_environments.jl @@ -158,6 +158,15 @@ function eachcoordinate(x::BPEnv, dirs) return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) end +## Bipartite check +function _is_bipartite(env::BPEnv) + (size(env, 2) == size(env, 3) == 2) || (return false) + for (d, c) in Iterators.product(axes(env, 1), axes(env, 3)) + (env[d, 1, c] == env[d, 2, _next(c, 2)]) || (return false) + end + return true +end + # conversion to CTMRGEnv """ CTMRGEnv(bp_env::BPEnv) diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 2966d4155..cd24cb111 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -138,6 +138,15 @@ TensorKit.spacetype(::Type{T}) where {E, T <: SUWeight{E}} = spacetype(E) TensorKit.sectortype(w::SUWeight) = sectortype(typeof(w)) TensorKit.sectortype(::Type{<:SUWeight{T}}) where {T} = sectortype(spacetype(T)) +## Bipartite check +function _is_bipartite(wts::SUWeight) + (size(wts, 2) == size(wts, 3) == 2) || (return false) + for (d, c) in Iterators.product(1:2, 1:2) + (wts[d, 1, c] == wts[d, 2, _next(c, 2)]) || (return false) + end + return true +end + ## (Approximate) equality function Base.:(==)(wts1::SUWeight, wts2::SUWeight) return wts1.data == wts2.data diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 172c56b05..0eb44444a 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -202,6 +202,14 @@ function InfinitePEPS(ρ::InfinitePEPO) ) end +## Bipartite check for PEPS/PEPO +function _is_bipartite(psi::InfiniteState) + (size(psi, 1) == size(psi, 2) == 2) || (return false) + for (c, h) in Iterators.product(1:2, 1:size(psi, 3)) + (psi[1, c, h] == psi[2, _next(c, 2), h]) || (return false) + end + return true +end ## Vector interface diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl index 941c8cbea..350a4f205 100644 --- a/test/bp/gaugefix.jl +++ b/test/bp/gaugefix.jl @@ -3,10 +3,13 @@ using Random using TensorKit using PEPSKit using PEPSKit: compare_weights, random_dual!, twistdual +using PEPSKit: _next, _is_bipartite -@testset "Compare BP and SU ($S, posdef msgs = $h)" for (S, h) in - Iterators.product([U1Irrep, FermionParity], [true, false]) - unitcell = (2, 3) +@testset "BP vs SU ($S, bipartite = $(bipartite), posdef msgs = $h)" for + (S, bipartite, h) in Iterators.product( + [U1Irrep, FermionParity], [true, false], [true, false] + ) + unitcell = bipartite ? (2, 2) : (2, 3) elt = ComplexF64 maxiter, tol = 100, 1.0e-9 Random.seed!(52840679) @@ -32,25 +35,48 @@ using PEPSKit: compare_weights, random_dual!, twistdual end end Nspaces, Espaces = random_dual!(Nspaces), random_dual!(Espaces) + if bipartite + for c in 1:2 + cp1 = _next(c, 2) + Pspaces[2, c] = Pspaces[1, cp1] + Nspaces[2, c] = Nspaces[1, cp1] + Espaces[2, c] = Espaces[1, cp1] + end + end peps0 = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces) + if bipartite + for c in 1:2 + peps0[2, c] = copy(peps0[1, _next(c, 2)]) + end + end # start by gauging with SU peps1, wts1 = gauge_fix(peps0, SUGauge(; maxiter, tol)) for (a0, a1) in zip(peps0.A, peps1.A) @test space(a0) == space(a1) end + if bipartite + @test _is_bipartite(peps1) + @test _is_bipartite(wts1) + end normalize!.(wts1.data) # find BP fixed point and SUWeight - bp_alg = BeliefPropagation(; maxiter, tol, project_hermitian = h) + bp_alg = BeliefPropagation(; maxiter, tol, bipartite, project_hermitian = h) env = BPEnv(randn, elt, peps1; posdef = h) env, err = leading_boundary(env, peps1, bp_alg) + if bipartite + @test _is_bipartite(env) + end wts2 = SUWeight(env) normalize!.(wts2.data) @test compare_weights(wts1, wts2) < 1.0e-9 bpg_alg = BPGauge() peps2, XXinv = @constinferred gauge_fix(peps1, bpg_alg, env) + if bipartite + @test _is_bipartite(peps2) + end for (a1, a2) in zip(peps1.A, peps2.A) @test space(a1) == space(a2) end diff --git a/test/runtests.jl b/test/runtests.jl index bdc2f7801..a50cdfe9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -107,7 +107,7 @@ end @time @safetestset "Cluster truncation with projectors" begin include("timeevol/cluster_projectors.jl") end - @time @safetestset "Time evolution with site-dependent truncation" begin + @time @safetestset "Fixed-space and site-dependent truncation" begin include("timeevol/sitedep_truncation.jl") end @time @safetestset "Time evolution with FixedSpaceTruncation" begin diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl index a36f3de75..6919e71b1 100644 --- a/test/timeevol/sitedep_truncation.jl +++ b/test/timeevol/sitedep_truncation.jl @@ -1,70 +1,51 @@ using Test -using LinearAlgebra using Random using TensorKit using PEPSKit -using PEPSKit: NORTH, EAST, _next +using PEPSKit: _is_bipartite -function get_bonddims(peps::InfinitePEPS) - xdims = collect(dim(domain(t, EAST)) for t in peps.A) - ydims = collect(dim(domain(t, NORTH)) for t in peps.A) - return stack([xdims, ydims]; dims = 1) -end - -function get_bonddims(wts::SUWeight) - xdims = collect(dim(space(wt, 1)) for wt in wts[1, :, :]) - ydims = collect(dim(space(wt, 1)) for wt in wts[2, :, :]) - return stack([xdims, ydims]; dims = 1) -end +elt = Float64 +Nr, Nc = 2, 2 +Vps = fill(U1Space(1 / 2 => 1, -1 / 2 => 1), (Nr, Nc)) +Vns = [ + U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 2, -1 => 1)'; + U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 2, -1 => 1) +] +Ves1 = [ + U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2); + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)' +] +Ves2 = [ + U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2); + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(0 => 1, 1 => 2, -1 => 1)' +] +Venv = U1Space(0 => 2, 1 => 1, -1 => 1) +states = ( + InfinitePEPS(randn, elt, Vps, Vns, Ves1), + InfinitePEPO(randn, elt, Vps, Vns, Ves2), +) -@testset "Simple update: bipartite 2-site" begin - Nr, Nc = 2, 2 - ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0)) - Random.seed!(100) - peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc)) - # make state bipartite - for r in 1:2 - peps0.A[_next(r, 2), 2] = copy(peps0.A[r, 1]) +@testset "Simple update on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for + (state0, bipartite) in Iterators.product(states, (true, false)) + J2 = 0.5 + if bipartite + state0[2, 1] = copy(state0[1, 2]) + state0[2, 2] = copy(state0[1, 1]) + J2 = 0.0 end - env0 = SUWeight(peps0) - normalize!.(peps0.A, Inf) - # set trunc to be compatible with bipartite structure - bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1) - trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims)) - alg = SimpleUpdate(; trunc, bipartite = true) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims - # check bipartite structure is preserved - for col in 1:2 - cp1 = PEPSKit._next(col, 2) - @test ( - peps.A[1, col] == peps.A[2, cp1] && - env[1, 1, col] == env[1, 2, cp1] && - env[2, 1, col] == env[2, 2, cp1] - ) + ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(Nr, Nc); J1 = 1.0, J2, sublattice = false) + # converted internally to SiteDependentTruncation + alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite) + wts0 = SUWeight(state0) + state, wts, = time_evolve(state0, ham, 0.1, 1, alg, wts0) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) + end + for (wt, wt0) in zip(wts.data, wts0.data) + @test space(wt) == space(wt0) + end + if bipartite + @test _is_bipartite(state) + @test _is_bipartite(wts) end -end - -@testset "Simple update: generic 2-site and 3-site" begin - Nr, Nc = 3, 4 - Random.seed!(100) - peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc)) - normalize!.(peps0.A, Inf) - env0 = SUWeight(peps0) - # Site dependent truncation - bonddims = rand(2:8, 2, Nr, Nc) - @show bonddims - trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims)) - alg = SimpleUpdate(; trunc) - # 2-site SU - ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.0, sublattice = false) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims - # 3-site SU - ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.2, sublattice = false) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims end From 3498fea2d4501ef21a3bf539ed7ece176e2892e3 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Apr 2026 15:44:49 +0800 Subject: [PATCH 17/27] Improve efficiency of bond truncation --- .../contractions/bondenv/als_solve.jl | 230 ++++++++++-------- src/algorithms/truncation/bond_truncation.jl | 75 +++--- .../truncation/fullenv_truncation.jl | 15 +- test/bondenv/bond_truncate.jl | 47 ++-- 4 files changed, 211 insertions(+), 156 deletions(-) diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index af083f618..745478073 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -2,139 +2,151 @@ In the following, the names `Ra`, `Sa` etc comes from the fast full update article Physical Review B 92, 035142 (2015) =# - """ -$(SIGNATURES) +Contract the virtual legs between +``` + -- DX --a-- D --b-- DY -- + ↓ ↓ + da db +``` +""" +function _combine_ket(a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2}) where {T, S} + return @tensor ket[DX DY; da db] := a[DX da; D] * b[D; db DY] +end +function _combine_ket(a::MPSTensor, b::MPSTensor) + return @tensor ket[DX DY; da db] := a[DX da; D] * b[D db; DY] +end -Construct the tensor -``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 Db0 - b -- DY0 -┘ - | | ↓ - |benv| db - | | ↓ - ┌---| |- DX1 Db1 - b† - DY1 -┐ - | └----┘ | - └-----------------------------------┘ -``` -""" -function _tensor_Ra(benv::BondEnv, b::MPSTensor) - return @autoopt @tensor Ra[DX1 Db1; DX0 Db0] := ( - benv[DX1 DY1; DX0 DY0] * b[Db0 db; DY0] * conj(b[Db1 db; DY1]) - ) +function _combine_ket_for_svd(a::MPSTensor, b::MPSTensor) + return @tensor ket[DX da; db DY] := a[DX da; D] * b[D db; DY] end """ -$(SIGNATURES) - -Construct the tensor +Construct the norm with bra bond tensors removed ``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 -- (a2 b2) -- DY0 --┘ - | | ↓ ↓ - |benv| da db - | | ↓ - ┌---| |- DX1 Db1 -- b† - DY1 --┐ - | └----┘ | - └-----------------------------------┘ + ┌benv-------┐ + ├---a---b---┤ + | ↓ ↓ | + ├-- --┤ + └-----------┘ ``` """ -function _tensor_Sa( - benv::BondEnv, b::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sa[DX1 da; Db1] := ( - benv[DX1 DY1; DX0 DY0] * conj(b[Db1 db; DY1]) * a2b2[DX0 DY0; da db] - ) +function _benv_ket(benv::BondEnv, ket::AbstractTensorMap{T, S, 2, 2}) where {T, S} + return benv * twistdual(ket, 1:2) end """ -$(SIGNATURES) + _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, i::Int) -Construct the tensor -``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 - a -- Da0 DY0 -┘ - | | ↓ - |benv| da - | | ↓ - ┌---| |- DX1 - a† - Da1 DY1 -┐ - | └----┘ | - └-----------------------------------┘ -``` -""" -function _tensor_Rb(benv::BondEnv, a::MPSTensor) - return @autoopt @tensor Rb[Da1 DY1; Da0 DY0] := ( - benv[DX1 DY1; DX0 DY0] * a[DX0 da; Da0] * conj(a[DX1 da; Da1]) - ) +Construct the bond environment around the `i`th bond tensor +in two-site ALS optimization. +``` + i = 1 i = 2 + ┌benv-------┐ ┌benv-------┐ + ├-- --b---┤ ├---a-- --┤ + | ↓ | | ↓ | + ├-- --b̄---┤ ├---ā-- --┤ + └-----------┘ └-----------┘ +``` +""" +_als_tensor_R(benv, xs, i::Int) = _als_tensor_R(benv, xs, Val(i)) +function _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, ::Val{1}) + return @tensor Ra[DX1 D1; DX0 D0] := + benv[DX1 DY1; DX0 DY0] * xs[2][D0 db; DY0] * conj(xs[2][D1 db; DY1]) +end +function _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, ::Val{2}) + return @tensor Rb[D1 DY1; D0 DY0] := + benv[DX1 DY1; DX0 DY0] * xs[1][DX0 da; D0] * conj(xs[1][DX1 da; D1]) end """ -$(SIGNATURES) - -Construct the tensor +Calculate the 2-site norm ``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 -- (a2 b2) -- DY0 --┘ - | | ↓ ↓ - |benv| da db - | | ↓ - ┌---| |- DX1 -- a† - Da1 DY1 --┐ - | └----┘ | - └-----------------------------------┘ + ┌benv-------┐ + ├---a---b---┤ + | ↓ ↓ | + ├---ā---b̄---┤ + └-----------┘ ``` +using pre-calcuated partial contraction results. """ -function _tensor_Sb( - benv::BondEnv, a::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sb[Da1 db; DY1] := ( - benv[DX1 DY1; DX0 DY0] * conj(a[DX1 da; Da1]) * a2b2[DX0 DY0; da db] - ) +function _als_norm( + ket::AbstractTensorMap{T, S, 2, 2}, benv_ket::AbstractTensorMap{T, S, 2, 2} + ) where {T, S} + return @tensor benv_ket[DX1 DY1; da db] * conj(ket[DX1 DY1; da db]) +end +function _als_norm(a::MPSTensor, Ra::BondEnv) + return @tensor Ra[DX1 D1; DX0 D0] * a[DX0 da; D0] * conj(a[DX1 da; D1]) end """ -$(SIGNATURES) + _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, i::Int + ) where {T <: Number, S <: ElementarySpace} -Calculate the inner product +Construct the overlap but with one of the bra bond tensor removed. ``` - ┌--------------------------------┐ - | ┌----┐ | - └---| |- DX0 - (a2 b2) - DY0 -┘ - | | ↓ ↓ - |benv| da db - | | ↓ ↓ - ┌---| |- DX1 - (a1 b1)†- DY1 -┐ - | └----┘ | - └--------------------------------┘ + i = 1 i = 2 + ┌benv-------┐ ┌benv-------┐ + ├---a₂==b₂--┤ ├---a₂==b₂--┤ + | ↓ ↓ | | ↓ ↓ | + ├-- --b̄---┤ ├---ā-- --┤ + └-----------┘ └-----------┘ ``` +The ket part is provided by the partial contraction `benv_ket`. """ -function inner_prod( - benv::BondEnv, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} +_als_tensor_S(benv_ket, xs, i::Int) = _als_tensor_S(benv_ket, xs, Val(i)) +function _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, ::Val{1} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * - conj(a1b1[DX1 DY1; da db]) * a2b2[DX0 DY0; da db] + return @tensor Sa[DX1 da; D1] := + benv_ket[DX1 DY1; da db] * conj(xs[2][D1 db; DY1]) +end +function _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, ::Val{2} + ) where {T <: Number, S <: ElementarySpace} + return @tensor contractcheck = true Sb[D1 db; DY1] := + benv_ket[DX1 DY1; da db] * conj(xs[1][DX1 da; D1]) end """ -$(SIGNATURES) +Calculate the inner product (overlap) +``` + ┌benv-------┐ + ├---a₂--b₂--┤ + | ↓ ↓ | + ├---ā---b̄---┤ + └-----------┘ +``` +using pre-calculated partial contraction results. +""" +function _als_overlap(a::MPSTensor, Sa::MPSTensor) + # applies to b, Sb as well + # @tensor Sb[D1 db; DY1] * conj(b[D1 db; DY1]) + return @tensor Sa[DX1 da; D1] * conj(a[DX1 da; D1]) +end -Contract the axis between `a` and `b` tensors +""" +Calculate the 2-site ALS inner product ⟨a₁,b₁|a₂,b₂⟩ ``` - -- DX - a - D - b - DY -- - ↓ ↓ - da db + ┌benv-------┐ + ├---a₂--b₂--┤ + | ↓ ↓ | + ├---ā₁--b̄₁--┤ + └-----------┘ ``` +where `|bra⟩ = |a₁,b₁⟩` and `|ket⟩ = |a₂,b₂⟩`, +with virtual leg between a, b contracted. """ -function _combine_ab( - a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2} +function inner_prod( + benv::BondEnv, bra::AbstractTensorMap{T, S, 2, 2}, + ket::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @tensor ab[DX DY; da db] := a[DX da; D] * b[D; db DY] -end -function _combine_ab(a::MPSTensor, b::MPSTensor) - return @tensor ab[DX DY; da db] := a[DX da; D] * b[D db; DY] + return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * + conj(bra[DX1 DY1; da db]) * ket[DX0 DY0; da db] end """ @@ -161,10 +173,30 @@ function cost_function_als(benv, ψ1, ψ2) return cost, fid end +# applies to Rb, Sb, b as well +# b22 is the pre-calculated untruncated norm +function cost_function_als(Ra::BondEnv, Sa::MPSTensor, a::MPSTensor, b22::Real) + b11 = real(_als_norm(a, Ra)) + b12 = _als_overlap(a, Sa) + cost = b11 + b22 - 2 * real(b12) + fid = abs2(b12) / abs(b11 * b22) + return cost, fid +end + """ $(SIGNATURES) Solve the equations `Rx x = Sx` with initial guess `x0`. + +In ALS over `a`, `b`, if we fix `b`, the cost function can +be expressed in the `Ra`, `Sa` tensors as +``` + f(a†,a) = a† Ra a - a† Sa - Sa† a + const +``` +Therefore `f` is minimized when +``` + ∂f/∂ā = Ra a - Sa = 0 +``` """ function _solve_als( Rx::AbstractTensorMap{T, S, N, N}, diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index 57fd68cc7..40d083e5e 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -18,8 +18,8 @@ The truncation algorithm can be constructed from the following keyword arguments * `tol::Float64=1e-9` : ALS converges when the relative change in bond SVD spectrum between two iterations is smaller than `tol`. * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. """ -@kwdef struct ALSTruncation - trunc::TruncationStrategy +@kwdef struct ALSTruncation{T <: TruncationStrategy} + trunc::T maxiter::Int = 50 tol::Float64 = 1.0e-9 check_interval::Int = 0 @@ -34,6 +34,20 @@ function _als_message( ) * @sprintf(" cost = %.3e, Δcost/cost0 = %.3e, |Δs| = %.4e.", cost, Δcost, Δs) end +""" +Initialize truncated bond tensors for 2-site ALS +""" +function _als_init_truncate( + ket2::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy + ) where {T, S} + a, s0, b = svd_trunc!(permute(ket2, ((1, 3), (4, 2)); copy = true); trunc) + a, b = absorb_s(a, s0, b) + # put b in MPS axis order + b = permute(b, ((1, 2), (3,))) + xs = [a, b] + return xs, s0 +end + """ bond_truncate(a::AbstractTensorMap{T,S,2,1}, b::AbstractTensorMap{T,S,1,2}, benv::BondEnv{T,S}, alg) -> U, S, V, info @@ -69,40 +83,39 @@ function bond_truncate( need_flip = isdual(space(b, 1)) time00 = time() verbose = (alg.check_interval > 0) - a2b2 = _combine_ab(a, b) - # initialize truncated a, b - perm_ab = ((1, 3), (4, 2)) - a, s0, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) - a, b = absorb_s(a, s0, b) - # put b in MPS axis order - b = permute(b, ((1, 2), (3,))) - ab = _combine_ab(a, b) + + # untruncated things + ket2 = _combine_ket(a, b) + benv_ket2 = _benv_ket(benv, ket2) + b22 = _als_norm(ket2, benv_ket2) + + # initialize truncated bond tensors and bond weight + xs, s0 = _als_init_truncate(ket2, alg.trunc) + + # initialize ALS cache + Rs = [_als_tensor_R(benv, xs, i) for i in 1:2] + Ss = [_als_tensor_S(benv_ket2, xs, i) for i in 1:2] + # cost function will be normalized by initial value - cost00, fid = cost_function_als(benv, ab, a2b2) + cost00, fid = cost_function_als(Rs[1], Ss[1], xs[1], b22) cost0, fid0, Δcost, Δfid, Δs = cost00, fid, NaN, NaN, NaN verbose && @info "ALS init" * _als_message(0, cost0, fid, Δcost, Δfid, Δs, 0.0) + for iter in 1:(alg.maxiter) time0 = time() - #= - Fixing `b`, the cost function can be expressed in the R, S tensors as - ``` - f(a†,a) = a† Ra a - a† Sa - Sa† a + const - ``` - `f` is minimized when - ∂f/∂ā = Ra a - Sa = 0 - =# - Ra = _tensor_Ra(benv, b) - Sa = _tensor_Sa(benv, b, a2b2) - a, info_a = _solve_als(Ra, Sa, a) - # Fixing `a`, solve for `b` from `Rb b = Sb` - Rb = _tensor_Rb(benv, a) - Sb = _tensor_Sb(benv, a, a2b2) - b, info_b = _solve_als(Rb, Sb, b) - @debug "Bond truncation info" info_a info_b - ab = _combine_ab(a, b) - cost, fid = cost_function_als(benv, ab, a2b2) + for (i, (Rx, Sx, x)) in enumerate(zip(Rs, Ss, xs)) + # TODO: option to use pinv + xs[i], info_x = _solve_als(Rx, Sx, x) + @debug "Bond truncation info $(i):" info_x + # update R, S for the next site + i_next = _next(i, 2) + Rs[i_next] = _als_tensor_R(benv, xs, i_next) + Ss[i_next] = _als_tensor_S(benv_ket2, xs, i_next) + end + # cost function and local fidelity + cost, fid = cost_function_als(Rs[1], Ss[1], xs[1], b22) # TODO: replace with truncated svdvals (without calculating u, vh) - _, s, _ = svd_trunc!(permute(ab, perm_ab); trunc = alg.trunc) + _, s, _ = svd_trunc!(_combine_ket_for_svd(xs...); trunc = alg.trunc) # fidelity, cost and normalized bond-s change s_nrm = norm(s0, Inf) Δs = _singular_value_distance(s, s0) / s_nrm @@ -129,7 +142,7 @@ function bond_truncate( end converge && break end - a, s, b = svd_trunc!(permute(_combine_ab(a, b), perm_ab); trunc = alg.trunc) + a, s, b = svd_trunc!(_combine_ket_for_svd(xs...); trunc = alg.trunc) a, b = absorb_s(a, s, b) if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl index fd9c685bd..4ff6bac2b 100644 --- a/src/algorithms/truncation/fullenv_truncation.jl +++ b/src/algorithms/truncation/fullenv_truncation.jl @@ -23,8 +23,8 @@ The truncation algorithm can be constructed from the following keyword arguments * [Glen Evenbly, Phys. Rev. B 98, 085155 (2018)](@cite evenbly_gauge_2018). """ -@kwdef struct FullEnvTruncation - trunc::TruncationStrategy +@kwdef struct FullEnvTruncation{T <: TruncationStrategy} + trunc::T maxiter::Int = 50 tol::Float64 = 1.0e-9 trunc_init::Bool = true @@ -75,13 +75,13 @@ function _fet_message( end """ - fullenv_truncate(benv::BondEnv{T,S}, b0::AbstractTensorMap{T,S,1,1}, alg::FullEnvTruncation) -> U, S, V, info + fullenv_truncate(b0, benv::BondEnv, alg::FullEnvTruncation) -> U, S, V, info Perform full environment truncation algorithm from [Phys. Rev. B 98, 085155 (2018)](@cite evenbly_gauge_2018) on `benv`. -Given a fixed state `|b0⟩` with bond matrix `b0` -and the corresponding positive-definite bond environment `benv`, +Given a fixed state `|b0⟩` with bond matrix `b0` and the +corresponding positive-definite bond environment `benv`, find the state `|b⟩` with truncated bond matrix `b = u s v†` that maximizes the fidelity (not normalized by `⟨b0|b0⟩`) ``` @@ -215,11 +215,12 @@ function fullenv_truncate( b1 = similar(b0) s0 = deepcopy(s) Δfid, Δs, fid, fid0 = NaN, NaN, 0.0, 0.0 + @tensor benv_b0[-1 -2] := benv[-1 -2; 3 4] * b0[3; 4] for iter in 1:(alg.maxiter) time0 = time() # update `← r - = ← s ← v† -` @tensor r[-1 -2] := s[-1; 1] * vh[1; -2] - @tensor p[-1 -2] := conj(u[1; -1]) * benv[1 -2; 3 4] * b0[3; 4] + @tensor p[-1 -2] := conj(u[1; -1]) * benv_b0[1 -2] @tensor B[-1 -2; -3 -4] := conj(u[1; -1]) * benv[1 -2; 3 -4] * u[3; -3] _linearmap_twist!(p) _linearmap_twist!(B) @@ -228,7 +229,7 @@ function fullenv_truncate( u, s, vh = svd_trunc(b1; trunc = alg.trunc) # update `- l ← = - u ← s ←` @tensor l[-1 -2] := u[-1; 1] * s[1; -2] - @tensor p[-1 -2] := conj(vh[-2; 2]) * benv[-1 2; 3 4] * b0[3; 4] + @tensor p[-1 -2] := conj(vh[-2; 2]) * benv_b0[-1 2] @tensor B[-1 -2; -3 -4] := conj(vh[-2; 2]) * benv[-1 2; -3 4] * vh[-4; 4] _linearmap_twist!(p) _linearmap_twist!(B) diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 25f56c200..684cd33c9 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -5,40 +5,49 @@ using TensorKit using PEPSKit using LinearAlgebra using KrylovKit -using PEPSKit: cost_function_als +using PEPSKit: bond_truncate, cost_function_als +using PEPSKit: _combine_ket, _combine_ket_for_svd Random.seed!(0) maxiter = 600 -check_interval = 20 -trunc = truncerror(; atol = 1.0e-10) & truncrank(8) -Vext = Vect[FermionParity](0 => 100, 1 => 100) -Vint = Vect[FermionParity](0 => 6, 1 => 6) -Vphy = Vect[FermionParity](0 => 1, 1 => 2) -perm_ab = ((1, 3), (4, 2)) -for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') - Vbond = Vbondl ⊗ Vbondr +check_interval = 30 +elt = Float64 +# simulating the situation of applying a 2-site gate +# to a bond with virtual dimension D, physical dimension d. +d, D = 2, 4 +trunc = truncerror(; atol = 1.0e-10) & truncrank(D) +Vphy = Vect[FermionParity](0 => div(d, 2), 1 => div(d, 2)) +Vqro = Vect[FermionParity](0 => div(d * D, 2), 1 => div(d * D, 2)) +# virtual dimension of gate MPO is d^2 +Vint = Vect[FermionParity](0 => div(d^2 * D, 2), 1 => div(d^2 * D, 2)) +for Vl in (Vqro, Vqro'), Vr in (Vqro, Vqro') # random positive-definite environment - Z = randn(Float64, Vext ← Vbond) + Vbond = Vl ⊗ Vr + Dext = dim(Vbond) + Vext = Vect[FermionParity](0 => div(Dext, 2) + 1, 1 => div(Dext, 2) + 1) + Z = randn(elt, Vext ← Vbond) + normalize!(Z, Inf) benv = Z' * Z - normalize!(benv, Inf) - # untruncated bond tensor - a2b2 = randn(Float64, Vbondl ⊗ Vbondr ← Vphy' ⊗ Vphy') - a2, s, b2 = svd_compact(permute(a2b2, perm_ab)) - a2, b2 = PEPSKit.absorb_s(a2, s, b2) + @info "Dimension of benv = $(Dext)" + # untruncated bond tensors + a2 = randn(elt, Vl ⊗ Vphy ← Vint) + b2 = randn(elt, Vint ← Vphy' ⊗ Vr') # bond tensor (truncated SVD initialization) - a0, s, b0 = svd_trunc(permute(a2b2, perm_ab); trunc = trunc) + a2b2 = _combine_ket(a2, b2) + a0, s, b0 = svd_trunc(permute(a2b2, ((1, 3), (4, 2))); trunc = trunc) a0, b0 = PEPSKit.absorb_s(a0, s, b0) - fid0 = cost_function_als(benv, PEPSKit._combine_ab(a0, b0), a2b2)[2] + fid0 = cost_function_als(benv, _combine_ket(a0, b0), a2b2)[2] @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() + # FET is slower when d is large for (label, alg) in ( ("ALS", ALSTruncation(; trunc, maxiter, check_interval)), ("FET", FullEnvTruncation(; trunc, maxiter, check_interval, trunc_init = false)), ) - a1, ss[label], b1, info = PEPSKit.bond_truncate(a2, b2, benv, alg) + a1, ss[label], b1, info = bond_truncate(a2, b2, benv, alg) @info "$label improved fidelity = $(info.fid)." # display(ss[label]) - @test info.fid ≈ cost_function_als(benv, PEPSKit._combine_ab(a1, b1), a2b2)[2] + @test info.fid ≈ cost_function_als(benv, _combine_ket(a1, b1), a2b2)[2] @test info.fid > fid0 end @test isapprox(ss["ALS"], ss["FET"], atol = 1.0e-3) From ea0a968ba9237a60981b7fb8fe6a4fd83fcb707e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Apr 2026 17:48:58 +0800 Subject: [PATCH 18/27] Make j1j2 finiteT test work for all allowed symmetries --- test/timeevol/j1j2_finiteT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl index 1e3845f39..b834b3522 100644 --- a/test/timeevol/j1j2_finiteT.jl +++ b/test/timeevol/j1j2_finiteT.jl @@ -10,7 +10,7 @@ const bm = [-0.08624893, -0.15688984, -0.21300888] function converge_env(state, χ::Int) trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12) - env0 = CTMRGEnv(ones, Float64, state, Vect[SU2Irrep](0 => 1)) + env0 = CTMRGEnv(ones, Float64, state, oneunit(spacetype(state))) env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10) return env end From 3cc776929dd3af694d1b00634e8fdd01e5b29488 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Apr 2026 21:21:25 +0800 Subject: [PATCH 19/27] Reorganize code to get bond tensors --- src/PEPSKit.jl | 1 + src/algorithms/time_evolution/apply_gate.jl | 93 +-------------------- src/algorithms/truncation/bond_tensor.jl | 90 ++++++++++++++++++++ 3 files changed, 92 insertions(+), 92 deletions(-) create mode 100644 src/algorithms/truncation/bond_tensor.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 44d3d42d8..cf2e89377 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -104,6 +104,7 @@ include("algorithms/ctmrg/c4v.jl") include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") +include("algorithms/truncation/bond_tensor.jl") include("algorithms/truncation/bond_truncation.jl") include("algorithms/time_evolution/apply_gate.jl") diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl index 5ffb98fee..7af55701a 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -23,98 +23,7 @@ end """ $(SIGNATURES) -Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. -When `A`, `B` are PEPSTensors, -``` - 2 1 1 - | | | - 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 - | ↘ | ↘ ↘ | - 4 1 3 2 2 3 -``` -When `A`, `B` are PEPOTensors, -- If `gate_ax = 1` -``` - 2 3 1 2 1 2 - ↘ | ↘ | ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ ↘ | - 5 1 4 2 2 4 -``` -- If `gate_ax = 2` -``` - 2 3 2 2 2 2 - ↘ | | ↘ ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ | ↘ - 5 1 4 1 4 1 -``` -""" -function _qr_bond(A::PT, B::PT; gate_ax::Int = 1, kwargs...) where {PT <: Union{PEPSTensor, PEPOTensor}} - @assert 1 <= gate_ax <= numout(A) - permA, permB, permX, permY = if A isa PEPSTensor - ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) - else - if gate_ax == 1 - ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - else - ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - end - end - X, a = left_orth!(permute(A, permA; copy = true); kwargs...) - Y, b = left_orth!(permute(B, permB; copy = true); kwargs...) - X, Y = permute(X, permX), permute(Y, permY) - b = permute(b, ((3, 2), (1,))) - return X, a, b, Y -end - -""" -$(SIGNATURES) - -Reconstruct the tensors connected by a bond from their `_qr_bond` results. -For PEPSTensors, -``` - -2 -2 - | | - -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 - | ↘ ↘ | - -4 -1 -1 -4 -``` -For PEPOTensors -``` - -2 -3 -2 -3 - ↘ | ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ ↘ | - -5 -1 -1 -5 - - -3 -2 -2 -3 - | ↘ ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ | ↘ - -5 -1 -5 -1 -``` -""" -function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) - @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] - @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] - return A, B -end -function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) - if !isdual(space(a, 2)) - @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] - else - @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] - end - return A, B -end - -""" -$(SIGNATURES) - -Apply 2-site `gate` on the reduced matrices `a`, `b` +Apply 2-site `gate` on the reduced bond tensors `a`, `b` ``` -1← a --- 3 --- b ← -4 -2 -3 ↓ ↓ ↓ ↓ diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl new file mode 100644 index 000000000..89bfc86ee --- /dev/null +++ b/src/algorithms/truncation/bond_tensor.jl @@ -0,0 +1,90 @@ +""" +$(SIGNATURES) + +Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. +When `A`, `B` are PEPSTensors, +``` + 2 1 1 + | | | + 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 + | ↘ | ↘ ↘ | + 4 1 3 2 2 3 +``` +When `A`, `B` are PEPOTensors, +- If `gate_ax = 1` +``` + 2 3 1 2 1 2 + ↘ | ↘ | ↘ | + 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 + | ↘ | ↘ ↘ | + 5 1 4 2 2 4 +``` +- If `gate_ax = 2` +``` + 2 3 2 2 2 2 + ↘ | | ↘ ↘ | + 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 + | ↘ | ↘ | ↘ + 5 1 4 1 4 1 +``` +""" +function _qr_bond(A::PT, B::PT; gate_ax::Int = 1, kwargs...) where {PT <: Union{PEPSTensor, PEPOTensor}} + @assert 1 <= gate_ax <= numout(A) + permA, permB, permX, permY = if A isa PEPSTensor + ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) + else + if gate_ax == 1 + ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) + else + ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) + end + end + X, a = left_orth!(permute(A, permA; copy = true); kwargs...) + Y, b = left_orth!(permute(B, permB; copy = true); kwargs...) + X, Y = permute(X, permX), permute(Y, permY) + b = permute(b, ((3, 2), (1,))) + return X, a, b, Y +end + +""" +$(SIGNATURES) + +Reconstruct the tensors connected by a bond from their `_qr_bond` results. +For PEPSTensors, +``` + -2 -2 + | | + -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 + | ↘ ↘ | + -4 -1 -1 -4 +``` +For PEPOTensors +``` + -2 -3 -2 -3 + ↘ | ↘ | + -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 + | ↘ ↘ | + -5 -1 -1 -5 + + -3 -2 -2 -3 + | ↘ ↘ | + -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 + | ↘ | ↘ + -5 -1 -5 -1 +``` +""" +function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) + @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] + @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] + return A, B +end +function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) + if !isdual(space(a, 2)) + @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] + @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] + else + @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] + @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] + end + return A, B +end From 5a901940cbc8cbcbb6b61654472a8b422f5a2f23 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Apr 2026 21:45:47 +0800 Subject: [PATCH 20/27] Improve NTU struct type stability --- src/algorithms/time_evolution/ntupdate.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index ace6efee3..b64e0ca98 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -11,14 +11,16 @@ Reference: - Physical Review B 104, 094411 (2021) - Physical Review B 106, 195105 (2022) """ -@kwdef struct NeighbourUpdate <: TimeEvolution +@kwdef struct NeighbourUpdate{ + TR <: Union{ALSTruncation, FullEnvTruncation}, + BE <: NeighbourEnv, + } <: TimeEvolution "Bond truncation algorithm after applying time evolution gate" - opt_alg::Union{ALSTruncation, FullEnvTruncation} = - ALSTruncation(; trunc = truncerror(; atol = 1.0e-10)) + opt_alg::TR = ALSTruncation(; trunc = truncerror(; atol = 1.0e-10)) "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" imaginary_time::Bool = true "Algorithm to construct NTU bond environment." - bondenv_alg::NeighbourEnv = NNEnv() + bondenv_alg::BE = NNEnv() "When true, fix gauge of bond environment" fixgauge::Bool = true "When true, assume bipartite unit cell structure" From 3d1c3f95bef3033c5de96a5e18b42f99e50a2f4c Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 25 Apr 2026 11:27:19 +0800 Subject: [PATCH 21/27] Fix FixedSpaceTruncation again --- src/algorithms/time_evolution/ntupdate.jl | 10 ++++- .../time_evolution/ntupdate3site.jl | 21 ++-------- src/algorithms/time_evolution/simpleupdate.jl | 4 -- .../time_evolution/simpleupdate3site.jl | 9 +--- test/runtests.jl | 3 -- test/timeevol/fixedspacetruncation.jl | 42 ------------------- test/timeevol/sitedep_truncation.jl | 22 ++++++++++ 7 files changed, 35 insertions(+), 76 deletions(-) delete mode 100644 test/timeevol/fixedspacetruncation.jl diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index b64e0ca98..3ac1507b5 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -58,6 +58,11 @@ function TimeEvolver( _timeevol_sanity_check(psi0, physicalspace(H), alg) dt′ = _get_dt(psi0, dt, alg.imaginary_time) gate = trotterize(H, dt′; symmetrize_gates) + # convert FixedSpaceTruncation to site-dependent `truncspace`s + if alg.opt_alg.trunc isa FixedSpaceTruncation + trunc = _get_fixedspacetrunc(psi0) + @reset alg.opt_alg.trunc = trunc + end state = NTUState(0, t0, psi0) return TimeEvolver(alg, dt, nstep, gate, state) end @@ -71,7 +76,10 @@ function _ntu_iter( sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate ) @assert length(sites) == 2 - return _bond_truncate(state, wts, Tuple(sites), alg; gate) + Nr, Nc = size(state) + trunc = only(_get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc))) + alg′ = (@set alg.opt_alg.trunc = trunc) + return _bond_truncate(state, wts, Tuple(sites), alg′; gate) end """ diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl index 92cc98081..ded2447f0 100644 --- a/src/algorithms/time_evolution/ntupdate3site.jl +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -6,19 +6,12 @@ function _ntu_iter( sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate ) where {T <: AbstractTensorMap} Nr, Nc = size(state) + truncs = _get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc)) state, wts = copy(state), deepcopy(wts) Ms, _, invperms = _get_cluster(state, sites) flips = [isdual(space(M, 1)) for M in Ms[2:end]] _flip_virtuals!(Ms, flips) # flip virtual arrows in `Ms` to ← - truncs = _get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc)) - truncs = map(enumerate(truncs)) do (i, trunc) - return if trunc isa FixedSpaceTruncation - truncspace(space(Ms[i + 1], 1)) - else - trunc - end - end # apply gate MPO without truncation _apply_gatempo!(Ms, gate) @@ -73,23 +66,15 @@ function _bond_truncate( @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))" end - # (optional) apply the NN gate - opt_alg = alg.opt_alg + # (optional) apply the NN gate without truncation if !(gate === nothing) - trunc = if alg.opt_alg.trunc isa FixedSpaceTruncation - V = space(b, 1) - truncspace(isdual(V) ? flip(V) : V) - else - alg.opt_alg.trunc - end - @reset opt_alg.trunc = trunc a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) else a = permute(a, ((1, 2), (3,))) b = permute(b, ((1,), (2, 3))) end - a, s, b, info = bond_truncate(a, b, benv, opt_alg) + a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) A, B = _qr_bond_undo(X, a, b, Y) normalize!(A, Inf) normalize!(B, Inf) diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index ecebda1c5..2a79ebfec 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -109,10 +109,6 @@ function _su_iter!( # rotate bond, rev = _nn_bondrev(sites..., (Nr, Nc)) A, B = _bond_rotation.(Ms, bond[1], rev; inv = false) - if trunc isa FixedSpaceTruncation - V = west_virtualspace(B) - trunc = truncspace(isdual(V) ? flip(V) : V) - end # apply gate ϵ, s = 0.0, nothing gate_axs = alg.purified ? (1:1) : (1:2) diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index df0d767ab..2fc9d096f 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -167,20 +167,13 @@ function _su_iter!( sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) where {T <: AbstractTensorMap} Nr, Nc = size(state) + truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) Ms, open_vaxs, invperms = _get_cluster(state, sites, env) flips = [isdual(space(M, 1)) for M in Ms[2:end]] Vphys = [codomain(M, 2) for M in Ms] normalize!.(Ms, Inf) # flip virtual arrows in `Ms` to ← _flip_virtuals!(Ms, flips) - truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) - truncs = map(enumerate(truncs)) do (i, trunc) - return if trunc isa FixedSpaceTruncation - truncspace(space(Ms[i + 1], 1)) - else - trunc - end - end # apply gate MPOs and truncate gate_axs = alg.purified ? (1:1) : (1:2) wts, ϵs = nothing, nothing diff --git a/test/runtests.jl b/test/runtests.jl index a50cdfe9c..fb3848f85 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -110,9 +110,6 @@ end @time @safetestset "Fixed-space and site-dependent truncation" begin include("timeevol/sitedep_truncation.jl") end - @time @safetestset "Time evolution with FixedSpaceTruncation" begin - include("timeevol/fixedspacetruncation.jl") - end @time @safetestset "Transverse field Ising model at finite temperature" begin include("timeevol/tf_ising_finiteT.jl") end diff --git a/test/timeevol/fixedspacetruncation.jl b/test/timeevol/fixedspacetruncation.jl deleted file mode 100644 index 1777c0436..000000000 --- a/test/timeevol/fixedspacetruncation.jl +++ /dev/null @@ -1,42 +0,0 @@ -using Test -using Random -using TensorKit -using PEPSKit - -elt = Float64 -ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(2, 2); J1 = 1.0, J2 = 0.5, sublattice = false) -Vphy = physicalspace(ham) -Vvir = U1Space(0 => 1, -1 / 2 => 1, 1 / 2 => 1) -Vns = [ - U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 1, -1 => 2); - U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 1, -1 => 2) -] -Ves1 = [ - U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 2, -1 / 2 => 1, 3 / 2 => 1); - U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1) -] -Ves2 = fill(U1Space(0 => 1, 1 => 1, -1 => 2), (2, 2)) -Venv = U1Space(0 => 2, 1 => 1, -1 => 1) -states = [ - InfinitePEPS(randn, elt, Vphy, Vns, Ves1), - InfinitePEPO(randn, elt, Vphy, Vns, Ves2), -] - -@testset "Simple update on $(typeof(state0).name.wrapper)" for state0 in states - alg = SimpleUpdate(; trunc = FixedSpaceTruncation()) - wts0 = SUWeight(state0) - state, wts, = time_evolve(state0, ham, 0.1, 1, alg, wts0) - for (t, t0) in zip(state.A, state0.A) - @test space(t) == space(t0) - end -end - -@testset "Neighborhood tensor update on $(typeof(state0).name.wrapper)" for state0 in states - opt_alg = ALSTruncation(; trunc = FixedSpaceTruncation()) - alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv()) - evolver = TimeEvolver(state0, ham, 0.1, 1, alg) - state, = time_evolve(evolver) - for (t, t0) in zip(state.A, state0.A) - @test space(t) == space(t0) - end -end diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl index 6919e71b1..0efde3c5c 100644 --- a/test/timeevol/sitedep_truncation.jl +++ b/test/timeevol/sitedep_truncation.jl @@ -49,3 +49,25 @@ states = ( @test _is_bipartite(wts) end end + +@testset "NTU on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for + (state0, bipartite) in Iterators.product(states, (true, false)) + J2 = 0.5 + if bipartite + state0[2, 1] = copy(state0[1, 2]) + state0[2, 2] = copy(state0[1, 1]) + J2 = 0.0 + end + ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(Nr, Nc); J1 = 1.0, J2, sublattice = false) + # converted internally to SiteDependentTruncation + opt_alg = ALSTruncation(; trunc = FixedSpaceTruncation()) + alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv(), bipartite) + state, info = time_evolve(TimeEvolver(state0, ham, 0.1, 1, alg)) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) + end + if bipartite + @test _is_bipartite(state) + @test _is_bipartite(info.wts) + end +end From 99e3e4865eb9646ea448478b9e18b918db3469f4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 25 Apr 2026 16:43:41 +0800 Subject: [PATCH 22/27] Require bond tensor to be an MPSTensor --- .../contractions/bondenv/als_solve.jl | 3 -- src/algorithms/truncation/bond_truncation.jl | 40 +++++++------------ test/bondenv/bond_truncate.jl | 3 +- 3 files changed, 17 insertions(+), 29 deletions(-) diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index 745478073..5d049f84c 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -10,9 +10,6 @@ Contract the virtual legs between da db ``` """ -function _combine_ket(a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2}) where {T, S} - return @tensor ket[DX DY; da db] := a[DX da; D] * b[D; db DY] -end function _combine_ket(a::MPSTensor, b::MPSTensor) return @tensor ket[DX DY; da db] := a[DX da; D] * b[D db; DY] end diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index 40d083e5e..1a95b8d34 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -49,7 +49,7 @@ function _als_init_truncate( end """ - bond_truncate(a::AbstractTensorMap{T,S,2,1}, b::AbstractTensorMap{T,S,1,2}, benv::BondEnv{T,S}, alg) -> U, S, V, info + bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg) -> U, S, V, info After time-evolving the reduced tensors `a` and `b` connected by a bond, truncate the bond dimension using the bond environment tensor `benv`. @@ -63,19 +63,14 @@ truncate the bond dimension using the bond environment tensor `benv`. └-----------------------┘ ``` The truncation algorithm `alg` can be either `FullEnvTruncation` or `ALSTruncation`. -The index order of `a` or `b` is +The index order of `a` or `b` follows `MPSTensor` convention. ``` 1 -a/b- 3 - ↓ a[1 2; 3] - 2 b[1; 2 3] + ↓ [1 2; 3] + 2 ``` """ -function bond_truncate( - a::AbstractTensorMap{T, S, 2, 1}, - b::AbstractTensorMap{T, S, 1, 2}, - benv::BondEnv{T, S}, - alg::ALSTruncation, - ) where {T <: Number, S <: ElementarySpace} +function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::ALSTruncation) # dual check of physical index @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @@ -144,18 +139,14 @@ function bond_truncate( end a, s, b = svd_trunc!(_combine_ket_for_svd(xs...); trunc = alg.trunc) a, b = absorb_s(a, s, b) + b = permute(b, ((1, 2), (3,))) if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end return a, s, b, (; fid, Δfid, Δs) end -function bond_truncate( - a::AbstractTensorMap{T, S, 2, 1}, - b::AbstractTensorMap{T, S, 1, 2}, - benv::BondEnv{T, S}, - alg::FullEnvTruncation, - ) where {T <: Number, S <: ElementarySpace} +function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::FullEnvTruncation) # dual check of physical index @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @@ -163,14 +154,14 @@ function bond_truncate( need_flip = isdual(space(b, 1)) #= initialize bond matrix using QR as `Ra Lb` - --- a == b --- ==> - Qa ← Ra == Rb ← Qb - + --- a == b --- ==> - Qa ← Ra == Rb → Qb - ↓ ↓ ↓ ↓ =# Qa, Ra = left_orth(a; positive = true) - Rb, Qb = right_orth(b; positive = true) - @assert !isdual(space(Ra, 1)) && !isdual(space(Qb, 1)) - @tensor b0[-1; -2] := Ra[-1 1] * Rb[1 -2] - #= initialize bond environment around `Ra Lb` + b = permute(b, ((3, 2), (1,)); copy = true) + Qb, Rb = left_orth!(b; positive = true) + @tensor b0[-1; -2] := Ra[-1 1] * Rb[-2 1] + #= initialize bond environment around `Ra Rb` ┌--------------------------------------┐ | ┌----┐ | @@ -182,15 +173,14 @@ function bond_truncate( | └----┘ | └--------------------------------------┘ =# - @tensor benv2[-1 -2; -3 -4] := ( - benv[1 2; 3 4] * conj(Qa[1 5 -1]) * conj(Qb[-2 6 2]) * Qa[3 5 -3] * Qb[-4 6 4] - ) + @tensor benv2[-1 -2; -3 -4] := benv[1 2; 3 4] * + conj(Qa[1 5 -1]) * conj(Qb[2 6 -2]) * Qa[3 5 -3] * Qb[4 6 -4] # optimize bond matrix u, s, vh, info = fullenv_truncate(b0, benv2, alg) u, vh = absorb_s(u, s, vh) # truncate a, b tensors with u, s, vh @tensor a[-1 -2; -3] := Qa[-1 -2 3] * u[3 -3] - @tensor b[-1; -2 -3] := vh[-1 1] * Qb[1 -2 -3] + @tensor b[-1 -2; -3] := vh[-1 1] * Qb[-3 -2 1] if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 684cd33c9..415d7c0d4 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -31,11 +31,12 @@ for Vl in (Vqro, Vqro'), Vr in (Vqro, Vqro') @info "Dimension of benv = $(Dext)" # untruncated bond tensors a2 = randn(elt, Vl ⊗ Vphy ← Vint) - b2 = randn(elt, Vint ← Vphy' ⊗ Vr') + b2 = randn(elt, Vint ⊗ Vphy ← Vr') # bond tensor (truncated SVD initialization) a2b2 = _combine_ket(a2, b2) a0, s, b0 = svd_trunc(permute(a2b2, ((1, 3), (4, 2))); trunc = trunc) a0, b0 = PEPSKit.absorb_s(a0, s, b0) + b0 = permute(b0, ((1, 2), (3,))) fid0 = cost_function_als(benv, _combine_ket(a0, b0), a2b2)[2] @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() From b6aeff3724cc14c0fd4dd664e592b721ed3a1302 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 25 Apr 2026 17:05:15 +0800 Subject: [PATCH 23/27] Replace `_qr_bond` with `bond_tensor` functions --- .../contractions/bondenv/benv_ctm.jl | 4 +- .../contractions/bondenv/benv_tools.jl | 2 +- .../contractions/bondenv/gaugefix.jl | 6 +- src/algorithms/time_evolution/apply_gate.jl | 6 +- .../time_evolution/ntupdate3site.jl | 6 +- src/algorithms/time_evolution/simpleupdate.jl | 6 +- src/algorithms/truncation/bond_tensor.jl | 158 ++++++++++-------- test/bondenv/benv_ctm.jl | 3 +- test/bondenv/benv_ntu.jl | 5 +- 9 files changed, 111 insertions(+), 85 deletions(-) diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 45262cc67..6910f30ea 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -9,8 +9,8 @@ Construct the environment (norm) tensor -1 0 1 2 ``` with `X` at position `[row, col]`. -`X, Y` are unitary tensors produced when finding the reduced site tensors -with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together). +`X, Y` are unitary tensors produced when finding the reduced site tensors, +and `XX = X' X` and `YY = Y' Y` (stacked together). Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index c3ae13ded..8b2ff362e 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -9,7 +9,7 @@ const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: Ele const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace} const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} # Orthogonal tensors obtained PEPSTensor/PEPOTensor -# with one physical leg factored out by `_qr_bond` +# with one physical leg factored out by `bond_tensor_...` const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace} const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace} diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 9df4c1599..a8edafc0a 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -44,9 +44,7 @@ Reference: - Physical Review B 92, 035142 (2015) """ function fixgauge_benv( - Z::AbstractTensorMap{T, S, 1, 2}, - a::AbstractTensorMap{T, S, 1, 2}, - b::AbstractTensorMap{T, S, 2, 1}, + Z::AbstractTensorMap{T, S, 1, 2}, a::MPSTensor, b::MPSTensor ) where {T <: Number, S <: ElementarySpace} @assert !isdual(space(Z, 1)) @assert !isdual(space(a, 2)) @@ -83,7 +81,7 @@ function fixgauge_benv( ↓ -1 =# - @plansor a[-1; -2 -3] := R[-1; 1] * a[1; -2 -3] + @plansor a[-1 -2; -3] := R[-1; 1] * a[1 -2; -3] @plansor b[-1 -2; -3] := b[-1 -2; 1] * L[-3; 1] @plansor Z[-1; -2 -3] := Z[-1; 1 2] * Rinv[1; -2] * Linv[2; -3] (isdual(space(R, 1)) == isdual(space(R, 2))) && twist!(a, 1) diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl index 7af55701a..67d4c63cf 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -34,10 +34,7 @@ Apply 2-site `gate` on the reduced bond tensors `a`, `b` -2 -3 -1← a --- 3 --- b ← -4 ``` """ -function _apply_gate( - a::AbstractTensorMap, b::AbstractTensorMap, - gate::NNGate, trunc::TruncationStrategy - ) +function _apply_gate(a::MPSTensor, b::MPSTensor, gate::NNGate, trunc::TruncationStrategy) V = space(b, 1) need_flip = isdual(V) if isdual(space(a, 2)) @@ -50,5 +47,6 @@ function _apply_gate( if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end + b = permute(b, ((1, 2), (3,))) return a, s, b, ϵ end diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl index ded2447f0..06d573774 100644 --- a/src/algorithms/time_evolution/ntupdate3site.jl +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -54,7 +54,8 @@ function _bond_truncate( A, B = state2[row, col], state2[row, cp1] # create bond environment - X, a, b, Y = _qr_bond(A, B; trunc = trunctol(; rtol = 1.0e-12)) + X, a = bond_tensor_first(A; trunc = trunctol(; rtol = 1.0e-12)) + Y, b = bond_tensor_last(B; trunc = trunctol(; rtol = 1.0e-12)) benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg) @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" if alg.fixgauge @@ -75,7 +76,8 @@ function _bond_truncate( end a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) - A, B = _qr_bond_undo(X, a, b, Y) + A = undo_bond_tensor_first(X, a) + B = undo_bond_tensor_last(Y, b) normalize!(A, Inf) normalize!(B, Inf) normalize!(s, Inf) diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 2a79ebfec..471bafe06 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -113,10 +113,12 @@ function _su_iter!( ϵ, s = 0.0, nothing gate_axs = alg.purified ? (1:1) : (1:2) for gate_ax in gate_axs - X, a, b, Y = _qr_bond(A, B; gate_ax, positive = true) + X, a = bond_tensor_first(A; gate_ax, positive = true) + Y, b = bond_tensor_last(B; gate_ax, positive = true) a, s, b, ϵ′ = _apply_gate(a, b, gate, trunc) ϵ = max(ϵ, ϵ′) - A, B = _qr_bond_undo(X, a, b, Y) + A = undo_bond_tensor_first(X, a; gate_ax) + B = undo_bond_tensor_last(Y, b; gate_ax) end # rotate back A = _bond_rotation(A, bond[1], rev; inv = true) diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl index 89bfc86ee..a1c02887b 100644 --- a/src/algorithms/truncation/bond_tensor.jl +++ b/src/algorithms/truncation/bond_tensor.jl @@ -1,90 +1,114 @@ """ -$(SIGNATURES) +Given the first tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its next bond. -Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. -When `A`, `B` are PEPSTensors, +For PEPSTensor, ``` - 2 1 1 - | | | - 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 - | ↘ | ↘ ↘ | - 4 1 3 2 2 3 + 1 + | + 4 - X ← 2 1 ← a - 3 + | ↘ + 3 2 ``` -When `A`, `B` are PEPOTensors, -- If `gate_ax = 1` +For PEPOTensor, ``` - 2 3 1 2 1 2 - ↘ | ↘ | ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ ↘ | - 5 1 4 2 2 4 -``` -- If `gate_ax = 2` -``` - 2 3 2 2 2 2 - ↘ | | ↘ ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ | ↘ - 5 1 4 1 4 1 + gate_ax = 1 gate_ax = 2 + + 1 2 2 2 + ↘ | | ↘ + 5 - X ← 3 1 ← a - 3 5 - X ← 3 1 ← a - 3 + | ↘ | ↘ + 4 2 4 1 ``` """ -function _qr_bond(A::PT, B::PT; gate_ax::Int = 1, kwargs...) where {PT <: Union{PEPSTensor, PEPOTensor}} - @assert 1 <= gate_ax <= numout(A) - permA, permB, permX, permY = if A isa PEPSTensor - ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) +function bond_tensor_first(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) + @assert gate_ax == 1 + X, a = left_orth!(permute(A, ((2, 4, 5), (1, 3)); copy = true); kwargs...) + X = permute(X, (1, 4, 2, 3)) + a = permute(a, ((1, 2), (3,))) + return X, a +end +function bond_tensor_first(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) + @assert 1 <= gate_ax <= 2 + X, a = if gate_ax == 1 + left_orth!(permute(A, ((2, 3, 5, 6), (1, 4)); copy = true); kwargs...) else - if gate_ax == 1 - ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - else - ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - end + left_orth!(permute(A, ((1, 3, 5, 6), (2, 4)); copy = true); kwargs...) end - X, a = left_orth!(permute(A, permA; copy = true); kwargs...) - Y, b = left_orth!(permute(B, permB; copy = true); kwargs...) - X, Y = permute(X, permX), permute(Y, permY) - b = permute(b, ((3, 2), (1,))) - return X, a, b, Y + X = permute(X, (1, 2, 5, 3, 4)) + a = permute(a, ((1, 2), (3,))) + return X, a end """ -$(SIGNATURES) +Undo the decomposition in `bond_tensor_first`. +""" +function undo_bond_tensor_first(X::PEPSOrth, a::MPSTensor; gate_ax::Integer = 1) + @assert gate_ax == 1 + return @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] +end +function undo_bond_tensor_first(X::PEPOOrth, a::MPSTensor; gate_ax::Integer = 1) + @assert 1 <= gate_ax <= 2 + if gate_ax == 1 + return @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] + else + return @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] + end +end -Reconstruct the tensors connected by a bond from their `_qr_bond` results. -For PEPSTensors, +""" +Given the last tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its previous bond. + +For PEPSTensor, ``` - -2 -2 - | | - -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 - | ↘ ↘ | - -4 -1 -1 -4 + 1 + | + 1 - b → 3 4 → Y - 2 + ↘ | + 2 3 ``` -For PEPOTensors +For PEPOTensor, ``` - -2 -3 -2 -3 - ↘ | ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ ↘ | - -5 -1 -1 -5 + gate_ax = 1 gate_ax = 2 - -3 -2 -2 -3 - | ↘ ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ | ↘ - -5 -1 -5 -1 + 1 2 2 2 + ↘ | ↘ | + 1 - b → 3 5 → Y - 3 1 - b → 3 5 → Y - 3 + ↘ | | ↘ + 2 4 4 1 ``` """ -function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) - @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] - @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] - return A, B +function bond_tensor_last(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) + @assert gate_ax == 1 + Y, b = left_orth!(permute(A, ((2, 3, 4), (1, 5)); copy = true); kwargs...) + Y = permute(Y, (1, 2, 3, 4)) + b = permute(b, ((3, 2), (1,))) + return Y, b +end +function bond_tensor_last(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) + @assert 1 <= gate_ax <= 2 + Y, b = if gate_ax == 1 + left_orth!(permute(A, ((2, 3, 4, 5), (1, 6)); copy = true); kwargs...) + else + left_orth!(permute(A, ((1, 3, 4, 5), (2, 6)); copy = true); kwargs...) + end + Y = permute(Y, (1, 2, 3, 4, 5)) + b = permute(b, ((3, 2), (1,))) + return Y, b +end + +""" +Undo the decomposition in `bond_tensor_last`. +""" +function undo_bond_tensor_last(Y::PEPSOrth, b::MPSTensor) + return @tensor A[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] end -function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) - if !isdual(space(a, 2)) - @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] +function undo_bond_tensor_last(Y::PEPOOrth, b::MPSTensor; gate_ax::Integer = 1) + @assert 1 <= gate_ax <= 2 + if gate_ax == 1 + return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] else - @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] + return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] end - return A, B end diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index 95ece1658..baf94e14a 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -42,7 +42,8 @@ function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) for row in 1:Nr, col in 1:Nc cp1 = PEPSKit._next(col, Nc) A, B = state.A[row, col], state.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) + X, a = PEPSKit.bond_tensor_first(A) + Y, b = PEPSKit.bond_tensor_last(B) benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) Z = PEPSKit.positive_approx(benv) # verify that gauge fixing can greatly reduce diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl index d3e2e4ffd..9b439fc85 100644 --- a/test/bondenv/benv_ntu.jl +++ b/test/bondenv/benv_ntu.jl @@ -15,10 +15,11 @@ function test_ntu_env( for row in 1:Nr, col in 1:Nc cp1 = PEPSKit._next(col, Nc) A, B = state.A[row, col], state.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) + X, a = PEPSKit.bond_tensor_first(A) + Y, b = PEPSKit.bond_tensor_last(B) @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY] benv = PEPSKit.bondenv_ntu(row, col, X, Y, state, env_alg) - # this is a result of `_qr_bond` + # this is a result of `bond_tensor_...` @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] # NTU bond environments are exact and should be positive definite @test benv' ≈ benv From bbc3bc4d65b6112197ef9c693ec1da5da34a02cb Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 25 Apr 2026 17:50:55 +0800 Subject: [PATCH 24/27] Fix `undo_bond_tensor_last` --- src/algorithms/truncation/bond_tensor.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl index a1c02887b..3ef7bdc0e 100644 --- a/src/algorithms/truncation/bond_tensor.jl +++ b/src/algorithms/truncation/bond_tensor.jl @@ -101,7 +101,8 @@ end """ Undo the decomposition in `bond_tensor_last`. """ -function undo_bond_tensor_last(Y::PEPSOrth, b::MPSTensor) +function undo_bond_tensor_last(Y::PEPSOrth, b::MPSTensor; gate_ax::Integer = 1) + @assert gate_ax == 1 return @tensor A[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] end function undo_bond_tensor_last(Y::PEPOOrth, b::MPSTensor; gate_ax::Integer = 1) From de5a5d5bd65c658195dc6db521dc47edbe8c35be Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 25 Apr 2026 20:59:44 +0800 Subject: [PATCH 25/27] Change bond_tensor return order --- .../time_evolution/ntupdate3site.jl | 8 +++---- src/algorithms/time_evolution/simpleupdate.jl | 8 +++---- src/algorithms/truncation/bond_tensor.jl | 24 +++++++++---------- test/bondenv/benv_ctm.jl | 4 ++-- test/bondenv/benv_ntu.jl | 4 ++-- 5 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl index 06d573774..713c789fc 100644 --- a/src/algorithms/time_evolution/ntupdate3site.jl +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -54,8 +54,8 @@ function _bond_truncate( A, B = state2[row, col], state2[row, cp1] # create bond environment - X, a = bond_tensor_first(A; trunc = trunctol(; rtol = 1.0e-12)) - Y, b = bond_tensor_last(B; trunc = trunctol(; rtol = 1.0e-12)) + a, X = bond_tensor_first(A; trunc = trunctol(; rtol = 1.0e-12)) + b, Y = bond_tensor_last(B; trunc = trunctol(; rtol = 1.0e-12)) benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg) @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" if alg.fixgauge @@ -76,8 +76,8 @@ function _bond_truncate( end a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) - A = undo_bond_tensor_first(X, a) - B = undo_bond_tensor_last(Y, b) + A = undo_bond_tensor_first(a, X) + B = undo_bond_tensor_last(b, Y) normalize!(A, Inf) normalize!(B, Inf) normalize!(s, Inf) diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 471bafe06..5c3d30a19 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -113,12 +113,12 @@ function _su_iter!( ϵ, s = 0.0, nothing gate_axs = alg.purified ? (1:1) : (1:2) for gate_ax in gate_axs - X, a = bond_tensor_first(A; gate_ax, positive = true) - Y, b = bond_tensor_last(B; gate_ax, positive = true) + a, X = bond_tensor_first(A; gate_ax, positive = true) + b, Y = bond_tensor_last(B; gate_ax, positive = true) a, s, b, ϵ′ = _apply_gate(a, b, gate, trunc) ϵ = max(ϵ, ϵ′) - A = undo_bond_tensor_first(X, a; gate_ax) - B = undo_bond_tensor_last(Y, b; gate_ax) + A = undo_bond_tensor_first(a, X; gate_ax) + B = undo_bond_tensor_last(b, Y; gate_ax) end # rotate back A = _bond_rotation(A, bond[1], rev; inv = true) diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl index 3ef7bdc0e..33df63d79 100644 --- a/src/algorithms/truncation/bond_tensor.jl +++ b/src/algorithms/truncation/bond_tensor.jl @@ -14,11 +14,11 @@ For PEPOTensor, ``` gate_ax = 1 gate_ax = 2 - 1 2 2 2 - ↘ | | ↘ + 1 2 2 2 + ↘ | | ↘ 5 - X ← 3 1 ← a - 3 5 - X ← 3 1 ← a - 3 - | ↘ | ↘ - 4 2 4 1 + | ↘ | ↘ + 4 2 4 1 ``` """ function bond_tensor_first(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) @@ -26,7 +26,7 @@ function bond_tensor_first(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) X, a = left_orth!(permute(A, ((2, 4, 5), (1, 3)); copy = true); kwargs...) X = permute(X, (1, 4, 2, 3)) a = permute(a, ((1, 2), (3,))) - return X, a + return a, X end function bond_tensor_first(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) @assert 1 <= gate_ax <= 2 @@ -37,17 +37,17 @@ function bond_tensor_first(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) end X = permute(X, (1, 2, 5, 3, 4)) a = permute(a, ((1, 2), (3,))) - return X, a + return a, X end """ Undo the decomposition in `bond_tensor_first`. """ -function undo_bond_tensor_first(X::PEPSOrth, a::MPSTensor; gate_ax::Integer = 1) +function undo_bond_tensor_first(a::MPSTensor, X::PEPSOrth; gate_ax::Integer = 1) @assert gate_ax == 1 return @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] end -function undo_bond_tensor_first(X::PEPOOrth, a::MPSTensor; gate_ax::Integer = 1) +function undo_bond_tensor_first(a::MPSTensor, X::PEPOOrth; gate_ax::Integer = 1) @assert 1 <= gate_ax <= 2 if gate_ax == 1 return @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] @@ -84,7 +84,7 @@ function bond_tensor_last(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) Y, b = left_orth!(permute(A, ((2, 3, 4), (1, 5)); copy = true); kwargs...) Y = permute(Y, (1, 2, 3, 4)) b = permute(b, ((3, 2), (1,))) - return Y, b + return b, Y end function bond_tensor_last(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) @assert 1 <= gate_ax <= 2 @@ -95,17 +95,17 @@ function bond_tensor_last(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) end Y = permute(Y, (1, 2, 3, 4, 5)) b = permute(b, ((3, 2), (1,))) - return Y, b + return b, Y end """ Undo the decomposition in `bond_tensor_last`. """ -function undo_bond_tensor_last(Y::PEPSOrth, b::MPSTensor; gate_ax::Integer = 1) +function undo_bond_tensor_last(b::MPSTensor, Y::PEPSOrth; gate_ax::Integer = 1) @assert gate_ax == 1 return @tensor A[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] end -function undo_bond_tensor_last(Y::PEPOOrth, b::MPSTensor; gate_ax::Integer = 1) +function undo_bond_tensor_last(b::MPSTensor, Y::PEPOOrth; gate_ax::Integer = 1) @assert 1 <= gate_ax <= 2 if gate_ax == 1 return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index baf94e14a..ee081520f 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -42,8 +42,8 @@ function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) for row in 1:Nr, col in 1:Nc cp1 = PEPSKit._next(col, Nc) A, B = state.A[row, col], state.A[row, cp1] - X, a = PEPSKit.bond_tensor_first(A) - Y, b = PEPSKit.bond_tensor_last(B) + a, X = PEPSKit.bond_tensor_first(A) + b, Y = PEPSKit.bond_tensor_last(B) benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) Z = PEPSKit.positive_approx(benv) # verify that gauge fixing can greatly reduce diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl index 9b439fc85..15c1654c4 100644 --- a/test/bondenv/benv_ntu.jl +++ b/test/bondenv/benv_ntu.jl @@ -15,8 +15,8 @@ function test_ntu_env( for row in 1:Nr, col in 1:Nc cp1 = PEPSKit._next(col, Nc) A, B = state.A[row, col], state.A[row, cp1] - X, a = PEPSKit.bond_tensor_first(A) - Y, b = PEPSKit.bond_tensor_last(B) + a, X = PEPSKit.bond_tensor_first(A) + b, Y = PEPSKit.bond_tensor_last(B) @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY] benv = PEPSKit.bondenv_ntu(row, col, X, Y, state, env_alg) # this is a result of `bond_tensor_...` From 2d90e0497db675970f036dfbdf70af1a9c06468e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 26 Apr 2026 10:01:19 +0800 Subject: [PATCH 26/27] Do not move phys leg to bond tensor for middle cluster sites --- .../contractions/bondenv/benv_ctm.jl | 4 +- .../contractions/bondenv/benv_ntu.jl | 12 +-- .../contractions/bondenv/gaugefix.jl | 75 ++++++++------ src/algorithms/time_evolution/ntupdate.jl | 2 +- .../time_evolution/ntupdate3site.jl | 46 +++++++-- src/algorithms/truncation/bond_tensor.jl | 98 +++++++++++++++++++ test/bondenv/benv_gaugefix.jl | 3 +- test/bondenv/benv_ntu.jl | 3 +- 8 files changed, 194 insertions(+), 49 deletions(-) diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 6910f30ea..8c6b1a199 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -24,8 +24,8 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` """ function bondenv_ctm( - row::Int, col::Int, X::T, Y::T, env::CTMRGEnv - ) where {T} + row::Int, col::Int, X::TX, Y::TY, env::CTMRGEnv + ) where {TX, TY} Nr, Nc = size(env.corners)[[2, 3]] cm1 = _prev(col, Nc) cp1 = _next(col, Nc) diff --git a/src/algorithms/contractions/bondenv/benv_ntu.jl b/src/algorithms/contractions/bondenv/benv_ntu.jl index 2c88c1cfe..ad4327f7b 100644 --- a/src/algorithms/contractions/bondenv/benv_ntu.jl +++ b/src/algorithms/contractions/bondenv/benv_ntu.jl @@ -34,8 +34,8 @@ Calculate the bond environment within "NTU-NN" approximation. ``` """ function bondenv_ntu( - row::Int, col::Int, X::T, Y::T, state::S, alg::NNEnv - ) where {T, S <: InfiniteState} + row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNEnv + ) where {TX, TY, S <: InfiniteState} neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)] m = collect_neighbors(state, row, col, neighbors) X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) @@ -78,8 +78,8 @@ Calculate the bond environment within "NTU-NN+" approximation. Dotted lines and ○ are splitted using SVD with `truncrank(1)`. """ function bondenv_ntu( - row::Int, col::Int, X::T, Y::T, state::S, alg::NNpEnv - ) where {T, S <: InfiniteState} + row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNpEnv + ) where {TX, TY, S <: InfiniteState} neighbors = [ (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (1, 2), (0, 2), (-1, 2), (-1, 1), (-1, 0), (0, -2), (2, 0), (2, 1), (0, 3), (-2, 1), (-2, 0), @@ -191,8 +191,8 @@ Calculates the bond environment within "NTU-NNN" approximation. ``` """ function bondenv_ntu( - row::Int, col::Int, X::T, Y::T, state::S, alg::NNNEnv - ) where {T, S <: InfiniteState} + row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNNEnv + ) where {TX, TY, S <: InfiniteState} neighbors = [ (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (1, 2), (0, 2), diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index a8edafc0a..8d964d395 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -90,40 +90,59 @@ function fixgauge_benv( end """ -When the (half) bond environment `Z` consists of -two `PEPSOrth` or `PEPOOrth` tensors `X`, `Y` as +Apply the gauge transformation `Rinv` for `Z` ``` ┌-----------------------┐ - | | - └---Z---(X)-- --(Y)---┘ + └---Z--(X)--Rinv-- ---┘ ↓ ``` -apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: +to `X`. For example, when `X` is a `PEPSTensor`, ``` - -1 -1 - | | - -4 - X - 1 - Rinv - -2 -4 - Linv - 1 - Y - -2 - | | - -3 -3 - - -2 -2 - | | - -5 - X - 1 - Rinv - -3 -5 - Linv - 1 - Y - -3 - | ╲ | ╲ - -4 -1 -4 -1 + -2 + | + -5 - X - 1 - Rinv - -3 + | ╲ + -4 -1 ``` """ -function _fixgauge_benvXY( - X::PEPSOrth, Y::PEPSOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, - ) - @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] - @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] - return X, Y +function _fixgauge_benvX(X::PEPSOrth, Rinv::MPSBondTensor) + return @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] end -function _fixgauge_benvXY( - X::PEPOOrth, Y::PEPOOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, - ) - @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3] - @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5] - return X, Y +function _fixgauge_benvX(X::PEPSTensor, Rinv::MPSBondTensor) + return @plansor X[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * Rinv[1; -3] +end +function _fixgauge_benvX(X::PEPOOrth, Rinv::MPSBondTensor) + return @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3] +end +function _fixgauge_benvX(X::PEPOTensor, Rinv::MPSBondTensor) + return @plansor X[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * Rinv[1; -4] +end + +""" +Apply the gauge transformation `Linv` for `Z` +``` + ┌-----------------------┐ + └---Z--- ---Linv--(Y)--┘ + ↓ +``` +to `Y`. For example, when `Y` is a `PEPSTensor`, +``` + -2 + | + -5 - Linv - 1 - Y - -3 + | ╲ + -4 -1 +``` +""" +function _fixgauge_benvY(Y::PEPSOrth, Linv::MPSBondTensor) + return @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] +end +function _fixgauge_benvY(Y::PEPSTensor, Linv::MPSBondTensor) + return @plansor Y[-1; -2 -3 -4 -5] := Y[-1; -2 -3 -4 1] * Linv[1; -5] +end +function _fixgauge_benvY(Y::PEPOOrth, Linv::MPSBondTensor) + return @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5] +end +function _fixgauge_benvY(Y::PEPOTensor, Linv::MPSBondTensor) + return @plansor Y[-1 -2; -3 -4 -5 -6] := Y[-1 -2; -3 -4 -5 1] * Linv[1; -6] end diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl index 3ac1507b5..7c4281762 100644 --- a/src/algorithms/time_evolution/ntupdate.jl +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -79,7 +79,7 @@ function _ntu_iter( Nr, Nc = size(state) trunc = only(_get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc))) alg′ = (@set alg.opt_alg.trunc = trunc) - return _bond_truncate(state, wts, Tuple(sites), alg′; gate) + return _bond_truncate(state, wts, Tuple(sites), (:first, :last), alg′; gate) end """ diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl index 713c789fc..fa50cdb21 100644 --- a/src/algorithms/time_evolution/ntupdate3site.jl +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -23,9 +23,13 @@ function _ntu_iter( # truncate each bond sequentially along the path info = (; fid = 1.0) - for (bondsites, trunc) in zip(zip(sites, Iterators.drop(sites, 1)), truncs) + nbond = length(sites) - 1 + for (i, bondsites) in enumerate(zip(sites, Iterators.drop(sites, 1))) + trunc = truncs[i] alg′ = (@set alg.opt_alg.trunc = trunc) - state, wts, info′ = _bond_truncate(state, wts, bondsites, alg′) + stype1 = (i == 1) ? :first : :middle + stype2 = (i == nbond) ? :last : :middle + state, wts, info′ = _bond_truncate(state, wts, bondsites, (stype1, stype2), alg′) # record the worst fidelity (info′.fid < info.fid) && (info = info′) end @@ -35,10 +39,14 @@ end """ Truncate a nearest neighbor bond between `site1` and `site2` after rotating the bond to standard x direction `A ← B`. + +`bondtype` takes values in (1, 2, 3), meaning that the current bond is +(the first, a middle, the last) bond in the updated cluster. """ function _bond_truncate( state::InfiniteState, wts::SUWeight, (site1, site2)::NTuple{2, CartesianIndex{2}}, + (stype1, stype2)::NTuple{2, Symbol}, alg::NeighbourUpdate; gate::Union{NNGate, Nothing} = nothing ) # rotate bond to standard x direction `A ← B` @@ -54,14 +62,26 @@ function _bond_truncate( A, B = state2[row, col], state2[row, cp1] # create bond environment - a, X = bond_tensor_first(A; trunc = trunctol(; rtol = 1.0e-12)) - b, Y = bond_tensor_last(B; trunc = trunctol(; rtol = 1.0e-12)) + qrtrunc = trunctol(; rtol = 1.0e-12) + a, X = if stype1 == :first + bond_tensor_first(A; trunc = qrtrunc) + else + @assert stype1 == :middle + bond_tensor_midnext(A; trunc = qrtrunc) + end + b, Y = if stype2 == :last + bond_tensor_last(B; trunc = qrtrunc) + else + @assert stype2 == :middle + bond_tensor_midprev(B; trunc = qrtrunc) + end benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg) @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" if alg.fixgauge Z = positive_approx(benv) Z, a, b, (Linv, Rinv) = fixgauge_benv(Z, a, b) - X, Y = _fixgauge_benvXY(X, Y, Linv, Rinv) + X = _fixgauge_benvX(X, Rinv) + Y = _fixgauge_benvY(Y, Linv) benv = Z' * Z @debug "cond(L) = $(LinearAlgebra.cond(Linv)); cond(R): $(LinearAlgebra.cond(Rinv))" @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))" @@ -70,14 +90,20 @@ function _bond_truncate( # (optional) apply the NN gate without truncation if !(gate === nothing) a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) + end + a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) + + A = if stype1 == :first + undo_bond_tensor_first(a, X) else - a = permute(a, ((1, 2), (3,))) - b = permute(b, ((1,), (2, 3))) + undo_bond_tensor_midnext(a, X) + end + B = if stype2 == :last + undo_bond_tensor_last(b, Y) + else + undo_bond_tensor_midprev(b, Y) end - a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) - A = undo_bond_tensor_first(a, X) - B = undo_bond_tensor_last(b, Y) normalize!(A, Inf) normalize!(B, Inf) normalize!(s, Inf) diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl index 33df63d79..6b11fea2b 100644 --- a/src/algorithms/truncation/bond_tensor.jl +++ b/src/algorithms/truncation/bond_tensor.jl @@ -113,3 +113,101 @@ function undo_bond_tensor_last(b::MPSTensor, Y::PEPOOrth; gate_ax::Integer = 1) return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] end end + +""" +Given a middle tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its next bond. + +For PEPSTensor, +``` + 2 + | + 5 - X ← 3 1 ← a - 3 + | ↘ ↘ + 4 1 (2) +``` +For PEPOTensor, +``` + 2 3 + ↘ | + 6 - X ← 4 1 ← a - 3 + | ↘ ↘ + 5 1 (2) +``` + +Here the physical leg on `a` is an auxiliary trivial leg. +""" +function bond_tensor_midnext(A::PEPSTensor; kwargs...) + X, a = left_orth!(permute(A, ((1, 2, 4, 5), (3,)); copy = true); kwargs...) + X = permute(X, ((1,), (2, 5, 3, 4))) + a = insertrightunit(a, 1) + return a, X +end +function bond_tensor_midnext(A::PEPOTensor; kwargs...) + X, a = left_orth!(permute(A, ((1, 2, 3, 5, 6), (4,)); copy = true); kwargs...) + X = permute(X, ((1, 2), (3, 6, 4, 5))) + a = insertrightunit(a, 1) + return a, X +end + +""" +Undo the decomposition in `bond_tensor_midprev`. +""" +function undo_bond_tensor_midnext(a::MPSTensor, X::PEPSTensor) + a = removeunit(a, 2) + return @tensor A[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * a[1; -3] +end +function undo_bond_tensor_midnext(a::MPSTensor, X::PEPOTensor) + a = removeunit(a, 2) + return @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * a[1; -4] +end + +""" +Given a middle tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its previous bond. + +For PEPSTensor, +``` + 2 + | + 1 - b → 3 5 → Y - 3 + ↘ | ↘ + (2) 4 1 +``` +For PEPOTensor, +``` + 2 3 + ↘ | + 1 - b → 3 6 → Y - 4 + ↘ | ↘ + (2) 5 1 +``` + +Here the physical leg on `a` is an auxiliary trivial leg. +""" +function bond_tensor_midprev(A::PEPSTensor; kwargs...) + Y, b = left_orth!(permute(A, ((1, 2, 3, 4), (5,)); copy = true); kwargs...) + Y = permute(Y, ((1,), (2, 3, 4, 5))) + b = permute(b, ((2,), (1,))) + b = insertrightunit(b, 1) + return b, Y +end +function bond_tensor_midprev(A::PEPOTensor; kwargs...) + Y, b = left_orth!(permute(A, ((1, 2, 3, 4, 5), (6,)); copy = true); kwargs...) + Y = permute(Y, ((1, 2), (3, 4, 5, 6))) + b = permute(b, ((2,), (1,))) + b = insertrightunit(b, 1) + return b, Y +end + +""" +Undo the decomposition in `bond_tensor_midprev`. +""" +function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPSTensor) + b = removeunit(b, 2) + return @tensor A[-1; -2 -3 -4 -5] := b[-5; 1] * Y[-1; -2 -3 -4 1] +end +function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPOTensor) + b = removeunit(b, 2) + return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6; 1] * Y[-1 -2; -3 -4 -5 1] +end diff --git a/test/bondenv/benv_gaugefix.jl b/test/bondenv/benv_gaugefix.jl index bdd5cecf7..0c5e2491d 100644 --- a/test/bondenv/benv_gaugefix.jl +++ b/test/bondenv/benv_gaugefix.jl @@ -37,7 +37,8 @@ for V1 in Vs, V2 in Vs, V3 in Vs @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] @test half ≈ half2 # test gauge transformation of X, Y - X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv) + X2 = _fixgauge_benvX(X, Rinv) + Y2 = _fixgauge_benvY(Y, Linv) @tensor Z2_[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X2[Xn Xe Xs Xw] * Y2[Yn Ye Ys Yw] @test Z2 ≈ Z2_ end diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl index 15c1654c4..6b5e11c81 100644 --- a/test/bondenv/benv_ntu.jl +++ b/test/bondenv/benv_ntu.jl @@ -41,7 +41,8 @@ function test_ntu_env( # verify gauge transformation of X, Y @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY] nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2) - X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv) + X2 = _fixgauge_benvX(X, Rinv) + Y2 = _fixgauge_benvY(Y, Linv) benv3 = PEPSKit.bondenv_ntu(row, col, X2, Y2, state, env_alg) benv3 *= norm(benv2, Inf) nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2) From 931a575b5615ead021d5754832bcdaf03799f9bf Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 26 Apr 2026 10:53:54 +0800 Subject: [PATCH 27/27] Fix bondenv tests --- test/bondenv/benv_gaugefix.jl | 10 +++++----- test/bondenv/benv_ntu.jl | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/test/bondenv/benv_gaugefix.jl b/test/bondenv/benv_gaugefix.jl index 0c5e2491d..4be536b0d 100644 --- a/test/bondenv/benv_gaugefix.jl +++ b/test/bondenv/benv_gaugefix.jl @@ -30,15 +30,15 @@ for V1 in Vs, V2 in Vs, V3 in Vs ↓ ↓ ↓ -1 -2 -3 =# - a = randn(ComplexF64, V1 ← Vphy' ⊗ V2) + a = randn(ComplexF64, V1 ⊗ Vphy ← V2) b = randn(ComplexF64, V2 ⊗ Vphy ← V3) - @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + @tensor half[:] := Z[-1; 1 3] * a[1 -2; 2] * b[2 -3; 3] Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) - @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @tensor half2[:] := Z2[-1; 1 3] * a2[1 -2; 2] * b2[2 -3; 3] @test half ≈ half2 # test gauge transformation of X, Y - X2 = _fixgauge_benvX(X, Rinv) - Y2 = _fixgauge_benvY(Y, Linv) + X2 = PEPSKit._fixgauge_benvX(X, Rinv) + Y2 = PEPSKit._fixgauge_benvY(Y, Linv) @tensor Z2_[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X2[Xn Xe Xs Xw] * Y2[Yn Ye Ys Yw] @test Z2 ≈ Z2_ end diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl index 6b5e11c81..24cb5ab0c 100644 --- a/test/bondenv/benv_ntu.jl +++ b/test/bondenv/benv_ntu.jl @@ -41,8 +41,8 @@ function test_ntu_env( # verify gauge transformation of X, Y @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY] nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2) - X2 = _fixgauge_benvX(X, Rinv) - Y2 = _fixgauge_benvY(Y, Linv) + X2 = PEPSKit._fixgauge_benvX(X, Rinv) + Y2 = PEPSKit._fixgauge_benvY(Y, Linv) benv3 = PEPSKit.bondenv_ntu(row, col, X2, Y2, state, env_alg) benv3 *= norm(benv2, Inf) nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2)