From 4c5fe137d564abe35a81f142db1f1d72c6b598f0 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 22 Jul 2026 17:03:17 +0800 Subject: [PATCH 01/21] Fix formatting --- test/algorithms/dynamical_dmrg.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl index 146c43af3..8e07757f4 100644 --- a/test/algorithms/dynamical_dmrg.jl +++ b/test/algorithms/dynamical_dmrg.jl @@ -60,4 +60,3 @@ end @test data ≈ predicted atol = 1.0e-8 end end - From e013be9b250104d9e2dc10b4c942cdaff73e616f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 22 Jul 2026 17:17:51 +0800 Subject: [PATCH 02/21] Add finite MPO-MPS zipper compression --- src/MPSKit.jl | 3 +- src/algorithms/changebonds/zipper.jl | 59 +++++++++++++++++++++++++++ src/operators/mpo.jl | 29 ++++++++++++++ src/states/abstractmps.jl | 19 +++++++-- test/algorithms/zipper.jl | 60 ++++++++++++++++++++++++++++ 5 files changed, 166 insertions(+), 4 deletions(-) create mode 100644 src/algorithms/changebonds/zipper.jl create mode 100644 test/algorithms/zipper.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 2ef8cf900..ee55dbfb8 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -36,7 +36,7 @@ export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster -export changebonds, changebonds! +export changebonds, changebonds!, zipper export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann @@ -157,6 +157,7 @@ include("algorithms/changebonds/changebonds.jl") include("algorithms/changebonds/optimalexpand.jl") include("algorithms/changebonds/vumpssvd.jl") include("algorithms/changebonds/svdcut.jl") +include("algorithms/changebonds/zipper.jl") include("algorithms/changebonds/randexpand.jl") include("algorithms/changebonds/sketchedexpand.jl") diff --git a/src/algorithms/changebonds/zipper.jl b/src/algorithms/changebonds/zipper.jl new file mode 100644 index 000000000..15d2c35ce --- /dev/null +++ b/src/algorithms/changebonds/zipper.jl @@ -0,0 +1,59 @@ +@doc """ + zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, trscheme; alg_svd=Defaults.alg_svd()) -> ψ′ + zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, alg::SvdCut) -> ψ′ + +Apply a finite open-boundary MPO `O` to a finite MPS `ψ` using a right-to-left +zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged +virtual bond is truncated immediately using `trscheme`. + +This is an unnormalized compression of `O * ψ`, comparable to +`changebonds(O * ψ, SvdCut(; trscheme); normalize=false)`, but without storing +the fully enlarged product on every site. +""" +function zipper end + +function _mpo_input_physicalspace(O) + return prod(x -> dual(space(O, x)), (numout(O) + 1):(numind(O) - 1)) +end + +function zipper( + O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, trscheme::TruncationStrategy; + alg_svd = Defaults.alg_svd() + ) + return zipper(O, ψ, SvdCut(; alg_svd, trscheme)) +end + +function zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, alg::SvdCut) + N = check_length(O, ψ) + if !isunitspace(left_virtualspace(O, 1)) || !isunitspace(right_virtualspace(O, N)) + throw(ArgumentError("zipper is only implemented for open-boundary MPOs")) + end + + ψ′ = copy(ψ) + T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) + A = TensorKit.similarstoragetype(eltype(ψ), T) + + Fᵣ = fuser(A, right_virtualspace(ψ′, N), right_virtualspace(O, N)) + local carry + alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme) + + As = map(N:-1:1) do i + Aψ = i == 1 ? ψ′.AC[1] : ψ′.AR[i] + physicalspace(Aψ) == _mpo_input_physicalspace(O[i]) || + throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) + Fₗ = fuser(A, left_virtualspace(ψ′, i), left_virtualspace(O, i)) + Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ) + i < N && (Aᶻ = Aᶻ * carry) + + if i == 1 + return Aᶻ + else + C, AR = right_gauge(Aᶻ, alg_gauge) + carry = C + Fᵣ = Fₗ + return AR + end + end + + return FiniteMPS(reverse(As); normalize = false, overwrite = true) +end diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 7bf0d16ca..3d385d5c9 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -271,6 +271,35 @@ function _fuse_mpo_mps(O::MPOTensor, A::MPSTensor, Fₗ, Fᵣ) @plansor A′[-1 -2; -3] := Fₗ[-1; 1 3] * A[1 2; 4] * O[3 -2; 2 5] * conj(Fᵣ[-3; 4 5]) return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end +@generated function _fuse_mpo_mps( + O::GenericMPOTensor{S, Nₒ}, A::GenericMPSTensor{S, Nₐ}, Fₗ, Fᵣ + ) where {S, Nₒ, Nₐ} + Nₒ == Nₐ || + return :(throw(ArgumentError("MPO input physical legs should match MPS physical legs"))) + + mps_left = 1 + mpo_left = 2 + mps_right = 3 + mpo_right = 4 + phys_in = 5:(Nₐ + 3) + + out_phys = -(2:Nₒ) + out_right = -(Nₒ + 1) + + t_out = tensorexpr(:A′, (-(1:Nₒ)...,), out_right) + t_left = tensorexpr(:Fₗ, -1, (mps_left, mpo_left)) + t_mps = tensorexpr(:A, (mps_left, phys_in...), mps_right) + t_mpo = tensorexpr(:O, (mpo_left, out_phys...), (phys_in..., mpo_right)) + t_right = tensorexpr(:Fᵣ, out_right, (mps_right, mpo_right)) + ex = macroexpand( + @__MODULE__, :(@plansor $t_out ≔ $t_left * $t_mps * $t_mpo * conj($t_right)) + ) + + return quote + $ex + return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ + end +end function Base.:*(mpo::FiniteMPO{<:MPOTensor}, x::AbstractTensorMap) @assert length(mpo) > 1 diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 40785968b..40a9a6643 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -2,13 +2,23 @@ Tensor types ===========================================================================================# +""" + GenericMPOTensor{S,N} + +Tensor type for representing local MPO tensors with matching numbers of codomain +and domain legs. The first codomain leg is the left virtual leg, the last domain +leg is the right virtual leg, the remaining codomain legs are output physical +legs, and the remaining domain legs are input physical legs. +""" +const GenericMPOTensor{S, N} = AbstractTensorMap{T, S, N, N} where {T} """ MPOTensor{S} -Tensor type for representing local MPO tensors, with the index convention `W ⊗ S ← N ⊗ E`, -where `N`, `E`, `S` and `W` denote the north, east, south and west virtual spaces respectively. +Tensor type for representing single-physical-leg local MPO tensors, with the +index convention `W ⊗ S ← N ⊗ E`, where `N`, `E`, `S` and `W` denote the north, +east, south and west virtual spaces respectively. """ -const MPOTensor{S} = AbstractTensorMap{T, S, 2, 2} where {T} +const MPOTensor{S} = GenericMPOTensor{S, 2} const MPSBondTensor{S} = AbstractTensorMap{T, S, 1, 1} where {T} const GenericMPSTensor{S, N} = AbstractTensorMap{T, S, N, 1} where {T} # some functions are also defined for "general mps tensors" (used in peps code) const MPSTensor{S} = GenericMPSTensor{S, 2} # the usual mps tensors on which we work @@ -214,6 +224,7 @@ Return the virtual space of the bond to the left of sites `pos`. function left_virtualspace end left_virtualspace(A::GenericMPSTensor) = space(A, 1) left_virtualspace(O::MPOTensor) = space(O, 1) +left_virtualspace(O::GenericMPOTensor) = space(O, 1) left_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(left_virtualspace, ψ), eachsite(ψ)) """ @@ -228,6 +239,7 @@ Return the virtual space of the bond to the right of site(s) `pos`. function right_virtualspace end right_virtualspace(A::GenericMPSTensor) = space(A, numind(A))' right_virtualspace(O::MPOTensor) = space(O, 4)' +right_virtualspace(O::GenericMPOTensor) = space(O, numind(O))' right_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(right_virtualspace, ψ), eachsite(ψ)) """ @@ -240,6 +252,7 @@ physicalspace(A::MPSTensor) = space(A, 2) physicalspace(A::GenericMPSTensor) = prod(x -> space(A, x), 2:(numind(A) - 1)) physicalspace(O::MPOTensor) = space(O, 2) physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = only(space(O, 2)) +physicalspace(O::GenericMPOTensor) = prod(x -> space(O, x), 2:numout(O)) physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ)) """ diff --git a/test/algorithms/zipper.jl b/test/algorithms/zipper.jl new file mode 100644 index 000000000..c64ecfbf5 --- /dev/null +++ b/test/algorithms/zipper.jl @@ -0,0 +1,60 @@ +println(" +----------------------------- +| Zipper tests | +----------------------------- +") + +using .TestSetup +using Test +using MPSKit +using TensorKit +using TensorKit: ℙ +using Random + +spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))] + +@testset "Finite MPO-MPS zipper $(spacetype(pspace))" for (pspace, Dspace) in spacelist + Random.seed!(1357) + L = 6 + Wspace = Dspace + Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)] + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) + O_copy = copy(O) + ψ_copy = copy(ψ) + + trscheme = trunctol(; atol = 1.0e-10) + ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) + got = zipper(O, ψ, trscheme) + + @test norm(ref - got) / norm(ref) < 1.0e-10 + @test norm(ψ - ψ_copy) < 1.0e-12 + @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) + + Dcut = 4 + got_tr = zipper(O, ψ, truncrank(Dcut)) + @test maximum(i -> dim(left_virtualspace(got_tr, i)), 2:length(got_tr)) <= Dcut +end + +@testset "Density-matrix FiniteMPS zipper $(spacetype(pcomp))" for (pcomp, Dspace) in [ + (ℙ^2 ⊗ (ℙ^2)', ℙ^6), + (Rep[SU₂](1 // 2 => 1) ⊗ Rep[SU₂](1 // 2 => 1)', Rep[SU₂](0 => 4, 1 => 3)), + ] + Random.seed!(9753) + L = 6 + Wspace = Dspace + Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)] + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pcomp ← pcomp ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, fill(pcomp, L), Dspace) + + trscheme = trunctol(; atol = 1.0e-10) + ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) + got = zipper(O, ψ, trscheme) + + @test numind(got.AC[L ÷ 2]) == 4 + @test norm(ref - got) / norm(ref) < 1.0e-10 +end From 85f6dbd83ceda54ac41d0776d9f86aa761decb59 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 23 Jul 2026 06:36:10 +0800 Subject: [PATCH 03/21] Change to `approximate` interface --- src/MPSKit.jl | 5 ++- src/algorithms/approximate/approximate.jl | 6 ++- .../{changebonds => approximate}/zipper.jl | 38 +++++++++---------- src/operators/mpo.jl | 29 -------------- src/states/abstractmps.jl | 19 ++-------- test/algorithms/zipper.jl | 25 +----------- 6 files changed, 30 insertions(+), 92 deletions(-) rename src/algorithms/{changebonds => approximate}/zipper.jl (57%) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index ee55dbfb8..ed369c5ea 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -36,8 +36,9 @@ export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster -export changebonds, changebonds!, zipper +export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand +export Zipper export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -157,7 +158,6 @@ include("algorithms/changebonds/changebonds.jl") include("algorithms/changebonds/optimalexpand.jl") include("algorithms/changebonds/vumpssvd.jl") include("algorithms/changebonds/svdcut.jl") -include("algorithms/changebonds/zipper.jl") include("algorithms/changebonds/randexpand.jl") include("algorithms/changebonds/sketchedexpand.jl") @@ -189,6 +189,7 @@ include("algorithms/statmech/idmrg.jl") include("algorithms/fidelity_susceptibility.jl") include("algorithms/approximate/approximate.jl") +include("algorithms/approximate/zipper.jl") include("algorithms/approximate/vomps.jl") include("algorithms/approximate/fvomps.jl") include("algorithms/approximate/idmrg.jl") diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 6967dcef0..b90f64770 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -4,10 +4,11 @@ approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ) approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) + approximate((O, ψ), algorithm) -> ψ′ Compute an approximation to the application of an operator `O` to the state `ψ` in the form -of an MPS `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, `ψ₀` is -approximated directly to `ψ` (i.e. `O` is taken to be the identity). +of an MPS, using initial guess `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, +`ψ₀` is approximated directly to `ψ` (i.e. `O` is taken to be the identity). **Not every algorithm supports every combination of arguments below** — see the per-algorithm notes at the end of this docstring before picking one. @@ -41,6 +42,7 @@ infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` ha |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:| | `DMRG` | single-site, fixes bond dim | `AbstractFiniteMPS` | ✅ | ✅ | | `DMRG2` | two-site, truncates via `trscheme` | `AbstractFiniteMPS` | ✅ | ✅ | +| `Zipper` | right-to-left streaming MPO-MPS compression | none, uses `(O, ψ)` directly | ❌ | ❌ | | `IDMRG` | single-site, thermodynamic limit | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `IDMRG2` | two-site, thermodynamic limit, needs unit cell ≥ 2 | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `VOMPS` | tangent-space truncation | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ❌ (out-of-place only) | diff --git a/src/algorithms/changebonds/zipper.jl b/src/algorithms/approximate/zipper.jl similarity index 57% rename from src/algorithms/changebonds/zipper.jl rename to src/algorithms/approximate/zipper.jl index 15d2c35ce..43167a471 100644 --- a/src/algorithms/changebonds/zipper.jl +++ b/src/algorithms/approximate/zipper.jl @@ -1,32 +1,30 @@ -@doc """ - zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, trscheme; alg_svd=Defaults.alg_svd()) -> ψ′ - zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, alg::SvdCut) -> ψ′ +""" + Zipper(; alg_svd=Defaults.alg_svd(), trscheme) + +Algorithm that approximates an open-boundary finite MPO-MPS product using a right-to-left +zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual +bond is truncated immediately using `trscheme`. + +Use as: -Apply a finite open-boundary MPO `O` to a finite MPS `ψ` using a right-to-left -zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged -virtual bond is truncated immediately using `trscheme`. + approximate((O, ψ), Zipper(; trscheme)) -This is an unnormalized compression of `O * ψ`, comparable to +This returns an unnormalized compression of `O * ψ`, comparable to `changebonds(O * ψ, SvdCut(; trscheme); normalize=false)`, but without storing the fully enlarged product on every site. """ -function zipper end - -function _mpo_input_physicalspace(O) - return prod(x -> dual(space(O, x)), (numout(O) + 1):(numind(O) - 1)) -end +@kwdef struct Zipper{S} <: Algorithm + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() -function zipper( - O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, trscheme::TruncationStrategy; - alg_svd = Defaults.alg_svd() - ) - return zipper(O, ψ, SvdCut(; alg_svd, trscheme)) + "algorithm used for truncation of the local gauge tensors" + trscheme::TruncationStrategy end -function zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, alg::SvdCut) +function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) N = check_length(O, ψ) if !isunitspace(left_virtualspace(O, 1)) || !isunitspace(right_virtualspace(O, N)) - throw(ArgumentError("zipper is only implemented for open-boundary MPOs")) + throw(ArgumentError("Zipper is only implemented for open-boundary MPOs")) end ψ′ = copy(ψ) @@ -39,7 +37,7 @@ function zipper(O::FiniteMPO{<:GenericMPOTensor}, ψ::FiniteMPS, alg::SvdCut) As = map(N:-1:1) do i Aψ = i == 1 ? ψ′.AC[1] : ψ′.AR[i] - physicalspace(Aψ) == _mpo_input_physicalspace(O[i]) || + physicalspace(Aψ) == physicalspace(O[i]) || throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) Fₗ = fuser(A, left_virtualspace(ψ′, i), left_virtualspace(O, i)) Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 3d385d5c9..7bf0d16ca 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -271,35 +271,6 @@ function _fuse_mpo_mps(O::MPOTensor, A::MPSTensor, Fₗ, Fᵣ) @plansor A′[-1 -2; -3] := Fₗ[-1; 1 3] * A[1 2; 4] * O[3 -2; 2 5] * conj(Fᵣ[-3; 4 5]) return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end -@generated function _fuse_mpo_mps( - O::GenericMPOTensor{S, Nₒ}, A::GenericMPSTensor{S, Nₐ}, Fₗ, Fᵣ - ) where {S, Nₒ, Nₐ} - Nₒ == Nₐ || - return :(throw(ArgumentError("MPO input physical legs should match MPS physical legs"))) - - mps_left = 1 - mpo_left = 2 - mps_right = 3 - mpo_right = 4 - phys_in = 5:(Nₐ + 3) - - out_phys = -(2:Nₒ) - out_right = -(Nₒ + 1) - - t_out = tensorexpr(:A′, (-(1:Nₒ)...,), out_right) - t_left = tensorexpr(:Fₗ, -1, (mps_left, mpo_left)) - t_mps = tensorexpr(:A, (mps_left, phys_in...), mps_right) - t_mpo = tensorexpr(:O, (mpo_left, out_phys...), (phys_in..., mpo_right)) - t_right = tensorexpr(:Fᵣ, out_right, (mps_right, mpo_right)) - ex = macroexpand( - @__MODULE__, :(@plansor $t_out ≔ $t_left * $t_mps * $t_mpo * conj($t_right)) - ) - - return quote - $ex - return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ - end -end function Base.:*(mpo::FiniteMPO{<:MPOTensor}, x::AbstractTensorMap) @assert length(mpo) > 1 diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 40a9a6643..40785968b 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -2,23 +2,13 @@ Tensor types ===========================================================================================# -""" - GenericMPOTensor{S,N} - -Tensor type for representing local MPO tensors with matching numbers of codomain -and domain legs. The first codomain leg is the left virtual leg, the last domain -leg is the right virtual leg, the remaining codomain legs are output physical -legs, and the remaining domain legs are input physical legs. -""" -const GenericMPOTensor{S, N} = AbstractTensorMap{T, S, N, N} where {T} """ MPOTensor{S} -Tensor type for representing single-physical-leg local MPO tensors, with the -index convention `W ⊗ S ← N ⊗ E`, where `N`, `E`, `S` and `W` denote the north, -east, south and west virtual spaces respectively. +Tensor type for representing local MPO tensors, with the index convention `W ⊗ S ← N ⊗ E`, +where `N`, `E`, `S` and `W` denote the north, east, south and west virtual spaces respectively. """ -const MPOTensor{S} = GenericMPOTensor{S, 2} +const MPOTensor{S} = AbstractTensorMap{T, S, 2, 2} where {T} const MPSBondTensor{S} = AbstractTensorMap{T, S, 1, 1} where {T} const GenericMPSTensor{S, N} = AbstractTensorMap{T, S, N, 1} where {T} # some functions are also defined for "general mps tensors" (used in peps code) const MPSTensor{S} = GenericMPSTensor{S, 2} # the usual mps tensors on which we work @@ -224,7 +214,6 @@ Return the virtual space of the bond to the left of sites `pos`. function left_virtualspace end left_virtualspace(A::GenericMPSTensor) = space(A, 1) left_virtualspace(O::MPOTensor) = space(O, 1) -left_virtualspace(O::GenericMPOTensor) = space(O, 1) left_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(left_virtualspace, ψ), eachsite(ψ)) """ @@ -239,7 +228,6 @@ Return the virtual space of the bond to the right of site(s) `pos`. function right_virtualspace end right_virtualspace(A::GenericMPSTensor) = space(A, numind(A))' right_virtualspace(O::MPOTensor) = space(O, 4)' -right_virtualspace(O::GenericMPOTensor) = space(O, numind(O))' right_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(right_virtualspace, ψ), eachsite(ψ)) """ @@ -252,7 +240,6 @@ physicalspace(A::MPSTensor) = space(A, 2) physicalspace(A::GenericMPSTensor) = prod(x -> space(A, x), 2:(numind(A) - 1)) physicalspace(O::MPOTensor) = space(O, 2) physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = only(space(O, 2)) -physicalspace(O::GenericMPOTensor) = prod(x -> space(O, x), 2:numout(O)) physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ)) """ diff --git a/test/algorithms/zipper.jl b/test/algorithms/zipper.jl index c64ecfbf5..84735af3f 100644 --- a/test/algorithms/zipper.jl +++ b/test/algorithms/zipper.jl @@ -27,34 +27,13 @@ spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 = trscheme = trunctol(; atol = 1.0e-10) ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) - got = zipper(O, ψ, trscheme) + got = approximate((O, ψ), Zipper(; trscheme)) @test norm(ref - got) / norm(ref) < 1.0e-10 @test norm(ψ - ψ_copy) < 1.0e-12 @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) Dcut = 4 - got_tr = zipper(O, ψ, truncrank(Dcut)) + got_tr = approximate((O, ψ), Zipper(; trscheme = truncrank(Dcut))) @test maximum(i -> dim(left_virtualspace(got_tr, i)), 2:length(got_tr)) <= Dcut end - -@testset "Density-matrix FiniteMPS zipper $(spacetype(pcomp))" for (pcomp, Dspace) in [ - (ℙ^2 ⊗ (ℙ^2)', ℙ^6), - (Rep[SU₂](1 // 2 => 1) ⊗ Rep[SU₂](1 // 2 => 1)', Rep[SU₂](0 => 4, 1 => 3)), - ] - Random.seed!(9753) - L = 6 - Wspace = Dspace - Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)] - O = FiniteMPO( - [rand(ComplexF64, Vspaces[i] ⊗ pcomp ← pcomp ⊗ Vspaces[i + 1]) for i in 1:L] - ) - ψ = FiniteMPS(rand, ComplexF64, fill(pcomp, L), Dspace) - - trscheme = trunctol(; atol = 1.0e-10) - ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) - got = zipper(O, ψ, trscheme) - - @test numind(got.AC[L ÷ 2]) == 4 - @test norm(ref - got) / norm(ref) < 1.0e-10 -end From 46ca0b743d05b5fbdb068dae7178d8169c7650bd Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Jul 2026 07:08:17 +0800 Subject: [PATCH 04/21] Remove unnecessary copy --- src/algorithms/approximate/zipper.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/algorithms/approximate/zipper.jl b/src/algorithms/approximate/zipper.jl index 43167a471..671b2cce0 100644 --- a/src/algorithms/approximate/zipper.jl +++ b/src/algorithms/approximate/zipper.jl @@ -27,19 +27,18 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) throw(ArgumentError("Zipper is only implemented for open-boundary MPOs")) end - ψ′ = copy(ψ) T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) A = TensorKit.similarstoragetype(eltype(ψ), T) - Fᵣ = fuser(A, right_virtualspace(ψ′, N), right_virtualspace(O, N)) + Fᵣ = fuser(A, right_virtualspace(ψ, N), right_virtualspace(O, N)) local carry alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme) As = map(N:-1:1) do i - Aψ = i == 1 ? ψ′.AC[1] : ψ′.AR[i] + Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] physicalspace(Aψ) == physicalspace(O[i]) || throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) - Fₗ = fuser(A, left_virtualspace(ψ′, i), left_virtualspace(O, i)) + Fₗ = fuser(A, left_virtualspace(ψ, i), left_virtualspace(O, i)) Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ) i < N && (Aᶻ = Aᶻ * carry) From 404fdbb9f3bdcc73a6f9f081770dfb50815a8848 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Jul 2026 07:22:06 +0800 Subject: [PATCH 05/21] Switch to left-to-right zipper --- src/algorithms/approximate/zipper.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/algorithms/approximate/zipper.jl b/src/algorithms/approximate/zipper.jl index 671b2cce0..1881cf053 100644 --- a/src/algorithms/approximate/zipper.jl +++ b/src/algorithms/approximate/zipper.jl @@ -1,7 +1,7 @@ """ Zipper(; alg_svd=Defaults.alg_svd(), trscheme) -Algorithm that approximates an open-boundary finite MPO-MPS product using a right-to-left +Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual bond is truncated immediately using `trscheme`. @@ -30,27 +30,27 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) A = TensorKit.similarstoragetype(eltype(ψ), T) - Fᵣ = fuser(A, right_virtualspace(ψ, N), right_virtualspace(O, N)) + Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) local carry alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme) - As = map(N:-1:1) do i + As = map(1:N) do i Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] physicalspace(Aψ) == physicalspace(O[i]) || throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) - Fₗ = fuser(A, left_virtualspace(ψ, i), left_virtualspace(O, i)) + Fᵣ = fuser(A, right_virtualspace(ψ, i), right_virtualspace(O, i)) Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ) - i < N && (Aᶻ = Aᶻ * carry) + i > 1 && (Aᶻ = _mul_front(carry, Aᶻ)) - if i == 1 + if i == N return Aᶻ else - C, AR = right_gauge(Aᶻ, alg_gauge) + AL, C, _ = left_gauge(Aᶻ, alg_gauge) carry = C - Fᵣ = Fₗ - return AR + Fₗ = Fᵣ + return AL end end - return FiniteMPS(reverse(As); normalize = false, overwrite = true) + return FiniteMPS(As; normalize = false, overwrite = true) end From 8757f954788f4143ab5278119d87b83122c97b71 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Jul 2026 07:24:03 +0800 Subject: [PATCH 06/21] Add reference --- docs/src/assets/mpskit.bib | 17 +++++++++++++++++ src/algorithms/approximate/zipper.jl | 12 ++++++------ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 296f4194b..cf9e3fd0b 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -385,6 +385,23 @@ @misc{shen2025 keywords = {Condensed Matter - Strongly Correlated Electrons,High Energy Physics - Theory} } +@article{sinha2024, + title = {Efficient Representation of Minimally Entangled Typical Thermal States in Two Dimensions via Projected Entangled Pair States}, + author = {Sinha, Aritra and Rams, Marek M. and Dziarmaga, Jacek}, + year = {2024}, + month = jan, + journal = {Physical Review B}, + volume = {109}, + number = {4}, + pages = {045136}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.109.045136}, + url = {https://link.aps.org/doi/10.1103/PhysRevB.109.045136}, + eprint = {2310.08533}, + archiveprefix = {arXiv}, + primaryclass = {quant-ph} +} + @article{tang2025, title = {Matrix Product State Fixed Points of Non-{{Hermitian}} Transfer Matrices}, author = {Tang, Wei and Verstraete, Frank and Haegeman, Jutho}, diff --git a/src/algorithms/approximate/zipper.jl b/src/algorithms/approximate/zipper.jl index 1881cf053..f8f12581a 100644 --- a/src/algorithms/approximate/zipper.jl +++ b/src/algorithms/approximate/zipper.jl @@ -1,17 +1,17 @@ """ - Zipper(; alg_svd=Defaults.alg_svd(), trscheme) +$(TYPEDEF) Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual bond is truncated immediately using `trscheme`. -Use as: +## Fields - approximate((O, ψ), Zipper(; trscheme)) +$(TYPEDFIELDS) -This returns an unnormalized compression of `O * ψ`, comparable to -`changebonds(O * ψ, SvdCut(; trscheme); normalize=false)`, but without storing -the fully enlarged product on every site. +## References + +* [Sinha et al. Phys. Rev. B 109 (2024)](@cite sinha2024) """ @kwdef struct Zipper{S} <: Algorithm "algorithm used for the singular value decomposition" From 2b8c493f9995cceab2cc2fb6eaf6bda82c0bfc44 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Jul 2026 07:33:56 +0800 Subject: [PATCH 07/21] Change constructors of `Zipper` --- src/algorithms/approximate/zipper.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/algorithms/approximate/zipper.jl b/src/algorithms/approximate/zipper.jl index f8f12581a..ca9d89fd9 100644 --- a/src/algorithms/approximate/zipper.jl +++ b/src/algorithms/approximate/zipper.jl @@ -9,16 +9,25 @@ bond is truncated immediately using `trscheme`. $(TYPEDFIELDS) +## Constructors + + Zipper(; trscheme, alg_svd=Defaults.alg_svd()) + Zipper(alg_gauge) + +Create a `Zipper` algorithm with the given truncated gauge algorithm, or by passing a +truncation scheme and singular value decomposition algorithm. + ## References -* [Sinha et al. Phys. Rev. B 109 (2024)](@cite sinha2024) +- [Sinha et al. Phys. Rev. B 109 (2024)](@cite sinha2024) """ -@kwdef struct Zipper{S} <: Algorithm - "algorithm used for the singular value decomposition" - alg_svd::S = Defaults.alg_svd() +struct Zipper{G} <: Algorithm + "algorithm used for gauging and truncating the local tensors" + alg_gauge::G +end - "algorithm used for truncation of the local gauge tensors" - trscheme::TruncationStrategy +function Zipper(; trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd()) + return Zipper(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) end function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) @@ -32,7 +41,6 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) local carry - alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme) As = map(1:N) do i Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] @@ -45,7 +53,7 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) if i == N return Aᶻ else - AL, C, _ = left_gauge(Aᶻ, alg_gauge) + AL, C, _ = left_gauge(Aᶻ, alg.alg_gauge) carry = C Fₗ = Fᵣ return AL From 9779e6289e6d341036fd92eedfa1a289122df8d1 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Jul 2026 08:52:18 +0800 Subject: [PATCH 08/21] Rename to `Zipup` and cite Stoudenmire instead --- docs/src/assets/mpskit.bib | 27 +++++++++---------- src/MPSKit.jl | 4 +-- src/algorithms/approximate/approximate.jl | 2 +- .../approximate/{zipper.jl => zipup.jl} | 20 +++++++------- test/algorithms/{zipper.jl => zipup.jl} | 8 +++--- 5 files changed, 30 insertions(+), 31 deletions(-) rename src/algorithms/approximate/{zipper.jl => zipup.jl} (69%) rename test/algorithms/{zipper.jl => zipup.jl} (81%) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index cf9e3fd0b..0bd289f7e 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -385,21 +385,20 @@ @misc{shen2025 keywords = {Condensed Matter - Strongly Correlated Electrons,High Energy Physics - Theory} } -@article{sinha2024, - title = {Efficient Representation of Minimally Entangled Typical Thermal States in Two Dimensions via Projected Entangled Pair States}, - author = {Sinha, Aritra and Rams, Marek M. and Dziarmaga, Jacek}, - year = {2024}, - month = jan, - journal = {Physical Review B}, - volume = {109}, - number = {4}, - pages = {045136}, - publisher = {American Physical Society}, - doi = {10.1103/PhysRevB.109.045136}, - url = {https://link.aps.org/doi/10.1103/PhysRevB.109.045136}, - eprint = {2310.08533}, +@article{stoudenmire2010, + title = {Minimally Entangled Typical Thermal State Algorithms}, + author = {Stoudenmire, E. M. and White, Steven R.}, + year = {2010}, + month = may, + journal = {New Journal of Physics}, + volume = {12}, + number = {5}, + pages = {055026}, + doi = {10.1088/1367-2630/12/5/055026}, + url = {https://doi.org/10.1088/1367-2630/12/5/055026}, + eprint = {1002.1305}, archiveprefix = {arXiv}, - primaryclass = {quant-ph} + primaryclass = {cond-mat.str-el} } @article{tang2025, diff --git a/src/MPSKit.jl b/src/MPSKit.jl index ed369c5ea..33c6b3c4f 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -38,7 +38,7 @@ export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand -export Zipper +export Zipup export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -189,7 +189,7 @@ include("algorithms/statmech/idmrg.jl") include("algorithms/fidelity_susceptibility.jl") include("algorithms/approximate/approximate.jl") -include("algorithms/approximate/zipper.jl") +include("algorithms/approximate/zipup.jl") include("algorithms/approximate/vomps.jl") include("algorithms/approximate/fvomps.jl") include("algorithms/approximate/idmrg.jl") diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index b90f64770..e7716208f 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -42,7 +42,7 @@ infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` ha |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:| | `DMRG` | single-site, fixes bond dim | `AbstractFiniteMPS` | ✅ | ✅ | | `DMRG2` | two-site, truncates via `trscheme` | `AbstractFiniteMPS` | ✅ | ✅ | -| `Zipper` | right-to-left streaming MPO-MPS compression | none, uses `(O, ψ)` directly | ❌ | ❌ | +| `Zipup` | left-to-right streaming MPO-MPS compression | none, uses `(O, ψ)` directly | ❌ | ❌ | | `IDMRG` | single-site, thermodynamic limit | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `IDMRG2` | two-site, thermodynamic limit, needs unit cell ≥ 2 | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `VOMPS` | tangent-space truncation | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ❌ (out-of-place only) | diff --git a/src/algorithms/approximate/zipper.jl b/src/algorithms/approximate/zipup.jl similarity index 69% rename from src/algorithms/approximate/zipper.jl rename to src/algorithms/approximate/zipup.jl index ca9d89fd9..f288a77e7 100644 --- a/src/algorithms/approximate/zipper.jl +++ b/src/algorithms/approximate/zipup.jl @@ -2,7 +2,7 @@ $(TYPEDEF) Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right -zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual +zip-up sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual bond is truncated immediately using `trscheme`. ## Fields @@ -11,29 +11,29 @@ $(TYPEDFIELDS) ## Constructors - Zipper(; trscheme, alg_svd=Defaults.alg_svd()) - Zipper(alg_gauge) + Zipup(; trscheme, alg_svd=Defaults.alg_svd()) + Zipup(alg_gauge) -Create a `Zipper` algorithm with the given truncated gauge algorithm, or by passing a +Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation scheme and singular value decomposition algorithm. ## References -- [Sinha et al. Phys. Rev. B 109 (2024)](@cite sinha2024) +- [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010) """ -struct Zipper{G} <: Algorithm +struct Zipup{G} <: Algorithm "algorithm used for gauging and truncating the local tensors" alg_gauge::G end -function Zipper(; trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd()) - return Zipper(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) +function Zipup(; trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd()) + return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) end -function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper) +function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) N = check_length(O, ψ) if !isunitspace(left_virtualspace(O, 1)) || !isunitspace(right_virtualspace(O, N)) - throw(ArgumentError("Zipper is only implemented for open-boundary MPOs")) + throw(ArgumentError("Zipup is only implemented for open-boundary MPOs")) end T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) diff --git a/test/algorithms/zipper.jl b/test/algorithms/zipup.jl similarity index 81% rename from test/algorithms/zipper.jl rename to test/algorithms/zipup.jl index 84735af3f..a3691b7a1 100644 --- a/test/algorithms/zipper.jl +++ b/test/algorithms/zipup.jl @@ -1,6 +1,6 @@ println(" ----------------------------- -| Zipper tests | +| Zipup tests | ----------------------------- ") @@ -13,7 +13,7 @@ using Random spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))] -@testset "Finite MPO-MPS zipper $(spacetype(pspace))" for (pspace, Dspace) in spacelist +@testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace) in spacelist Random.seed!(1357) L = 6 Wspace = Dspace @@ -27,13 +27,13 @@ spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 = trscheme = trunctol(; atol = 1.0e-10) ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) - got = approximate((O, ψ), Zipper(; trscheme)) + got = approximate((O, ψ), Zipup(; trscheme)) @test norm(ref - got) / norm(ref) < 1.0e-10 @test norm(ψ - ψ_copy) < 1.0e-12 @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) Dcut = 4 - got_tr = approximate((O, ψ), Zipper(; trscheme = truncrank(Dcut))) + got_tr = approximate((O, ψ), Zipup(; trscheme = truncrank(Dcut))) @test maximum(i -> dim(left_virtualspace(got_tr, i)), 2:length(got_tr)) <= Dcut end From 35e31f5a00e863642dd68ed7f712adaed590d548 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 28 Jul 2026 15:36:03 +0800 Subject: [PATCH 09/21] Two-sweep version of zipup --- src/algorithms/approximate/zipup.jl | 41 ++++++++++++++++++++--------- test/algorithms/zipup.jl | 33 +++++++++++++++++++---- 2 files changed, 57 insertions(+), 17 deletions(-) diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index f288a77e7..56b27e0f8 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -2,8 +2,8 @@ $(TYPEDEF) Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right -zip-up sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual -bond is truncated immediately using `trscheme`. +zip-up sweep, optionally followed by a second `changebonds` sweep. The MPO and MPS are +contracted one site at a time, and the enlarged virtual bond is truncated immediately. ## Fields @@ -12,22 +12,34 @@ $(TYPEDFIELDS) ## Constructors Zipup(; trscheme, alg_svd=Defaults.alg_svd()) - Zipup(alg_gauge) + Zipup(alg_zipup, [alg_finalize]) -Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a -truncation scheme and singular value decomposition algorithm. +Create a `Zipup` algorithm with the given [`SvdCut`](@ref), or by passing a truncation scheme +and singular value decomposition algorithm. If `alg_finalize` is provided, the state obtained +after the zip-up sweep is further compressed with `changebonds!`. + +Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more +permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_finalize` +to impose the final truncation. ## References - [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010) +- [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) """ -struct Zipup{G} <: Algorithm - "algorithm used for gauging and truncating the local tensors" - alg_gauge::G +struct Zipup{G <: SvdCut, F} <: Algorithm + "algorithm used for gauging and truncating the local tensors during the zip-up sweep" + alg_zipup::G + "algorithm used for the final locally gauged truncation pass; `nothing` skips this pass" + alg_finalize::F end -function Zipup(; trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd()) - return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) +Zipup(alg_zipup) = Zipup(alg_zipup, nothing) + +function Zipup(; + trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd() + ) + return Zipup(SvdCut(; alg_svd, trscheme)) end function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) @@ -41,6 +53,9 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) local carry + alg_zipup = MatrixAlgebraKit.TruncatedAlgorithm( + alg.alg_zipup.alg_svd, alg.alg_zipup.trscheme + ) As = map(1:N) do i Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] @@ -53,12 +68,14 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) if i == N return Aᶻ else - AL, C, _ = left_gauge(Aᶻ, alg.alg_gauge) + AL, C, _ = left_gauge(Aᶻ, alg_zipup) carry = C Fₗ = Fᵣ return AL end end - return FiniteMPS(As; normalize = false, overwrite = true) + ψ′ = FiniteMPS(As; normalize = false, overwrite = true) + return isnothing(alg.alg_finalize) ? + ψ′ : changebonds!(ψ′, alg.alg_finalize; normalize = false) end diff --git a/test/algorithms/zipup.jl b/test/algorithms/zipup.jl index a3691b7a1..4a8c1adac 100644 --- a/test/algorithms/zipup.jl +++ b/test/algorithms/zipup.jl @@ -11,9 +11,12 @@ using TensorKit using TensorKit: ℙ using Random -spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))] +spacelist = [ + (ℙ^4, ℙ^3, 4), + (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8), +] -@testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace) in spacelist +function _random_mpo_mps(pspace, Dspace) Random.seed!(1357) L = 6 Wspace = Dspace @@ -22,6 +25,11 @@ spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 = [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] ) ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) + return O, ψ +end + +@testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist + O, ψ = _random_mpo_mps(pspace, Dspace) O_copy = copy(O) ψ_copy = copy(ψ) @@ -32,8 +40,23 @@ spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 = @test norm(ref - got) / norm(ref) < 1.0e-10 @test norm(ψ - ψ_copy) < 1.0e-12 @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) +end + +@testset "Paeckel two-stage zip-up $(spacetype(pspace))" for (pspace, Dspace, Dcut) in spacelist + O, ψ = _random_mpo_mps(pspace, Dspace) + rtol = 1.0e-8 + final_trscheme = truncrank(Dcut) & truncerror(; rtol) + zipup_trscheme = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) + + ref_tr = changebonds(O * ψ, SvdCut(; trscheme = final_trscheme); normalize = false) + got_one_sweep = approximate((O, ψ), Zipup(; trscheme = final_trscheme)) + alg_zipup = SvdCut(; trscheme = zipup_trscheme) + alg_finalize = SvdCut(; trscheme = final_trscheme) + got_two_sweep = approximate((O, ψ), Zipup(alg_zipup, alg_finalize)) - Dcut = 4 - got_tr = approximate((O, ψ), Zipup(; trscheme = truncrank(Dcut))) - @test maximum(i -> dim(left_virtualspace(got_tr, i)), 2:length(got_tr)) <= Dcut + err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) + err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) + @test err_two_sweep < err_one_sweep / 2 + @test maximum(i -> dim(left_virtualspace(got_one_sweep, i)), 2:length(got_one_sweep)) <= Dcut + @test maximum(i -> dim(left_virtualspace(got_two_sweep, i)), 2:length(got_two_sweep)) <= Dcut end From 58bd0e55099c0925f00f84e2fdb1c18add06a2b5 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 29 Jul 2026 16:40:45 +0800 Subject: [PATCH 10/21] Manual zip-down instead of calling `changebonds` --- src/algorithms/approximate/zipup.jl | 56 ++++++++++++++++++++--------- test/algorithms/zipup.jl | 6 ++-- 2 files changed, 42 insertions(+), 20 deletions(-) diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index 56b27e0f8..f388184b9 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -2,7 +2,7 @@ $(TYPEDEF) Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right -zip-up sweep, optionally followed by a second `changebonds` sweep. The MPO and MPS are +zip-up sweep, optionally followed by a right-to-left zip-down sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual bond is truncated immediately. ## Fields @@ -12,14 +12,14 @@ $(TYPEDFIELDS) ## Constructors Zipup(; trscheme, alg_svd=Defaults.alg_svd()) - Zipup(alg_zipup, [alg_finalize]) + Zipup(alg_zipup, [alg_zipdown]) -Create a `Zipup` algorithm with the given [`SvdCut`](@ref), or by passing a truncation scheme -and singular value decomposition algorithm. If `alg_finalize` is provided, the state obtained -after the zip-up sweep is further compressed with `changebonds!`. +Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation +scheme and singular value decomposition algorithm. If `alg_zipdown` is provided, the state +obtained after the zip-up sweep is further compressed by sweeping back from right to left. Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more -permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_finalize` +permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` to impose the final truncation. ## References @@ -27,11 +27,14 @@ to impose the final truncation. - [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010) - [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) """ -struct Zipup{G <: SvdCut, F} <: Algorithm +struct Zipup{ + U <: MatrixAlgebraKit.TruncatedAlgorithm, + D <: Union{Nothing, MatrixAlgebraKit.TruncatedAlgorithm}, + } <: Algorithm "algorithm used for gauging and truncating the local tensors during the zip-up sweep" - alg_zipup::G + alg_zipup::U "algorithm used for the final locally gauged truncation pass; `nothing` skips this pass" - alg_finalize::F + alg_zipdown::D end Zipup(alg_zipup) = Zipup(alg_zipup, nothing) @@ -39,7 +42,7 @@ Zipup(alg_zipup) = Zipup(alg_zipup, nothing) function Zipup(; trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd() ) - return Zipup(SvdCut(; alg_svd, trscheme)) + return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) end function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) @@ -53,9 +56,6 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) local carry - alg_zipup = MatrixAlgebraKit.TruncatedAlgorithm( - alg.alg_zipup.alg_svd, alg.alg_zipup.trscheme - ) As = map(1:N) do i Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] @@ -68,14 +68,36 @@ function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) if i == N return Aᶻ else - AL, C, _ = left_gauge(Aᶻ, alg_zipup) + AL, C, _ = left_gauge(Aᶻ, alg.alg_zipup) carry = C Fₗ = Fᵣ return AL end end - ψ′ = FiniteMPS(As; normalize = false, overwrite = true) - return isnothing(alg.alg_finalize) ? - ψ′ : changebonds!(ψ′, alg.alg_finalize; normalize = false) + return isnothing(alg.alg_zipdown) ? + FiniteMPS(As; normalize = false, overwrite = true) : + _zipdown(As, alg.alg_zipdown) +end + +function _zipdown(As::Vector{A}, alg::MatrixAlgebraKit.TruncatedAlgorithm) where {A} + N = length(As) + N == 1 && return FiniteMPS(As; normalize = false, overwrite = true) + + ARs = Vector{Union{Missing, A}}(missing, N) + ALs = Vector{Union{Missing, A}}(missing, N) + ACs = Vector{Union{Missing, A}}(missing, N) + + local C + AC = As[N] + for i in N:-1:2 + C, AR, _ = right_gauge(AC, alg) + ARs[i] = AR + AC = _mul_tail(As[i - 1], C) + end + + B = typeof(C) + ACs[1] = AC + Cs = Vector{Union{Missing, B}}(missing, N + 1) + return FiniteMPS(ALs, ARs, ACs, Cs) end diff --git a/test/algorithms/zipup.jl b/test/algorithms/zipup.jl index 4a8c1adac..a4e515835 100644 --- a/test/algorithms/zipup.jl +++ b/test/algorithms/zipup.jl @@ -50,9 +50,9 @@ end ref_tr = changebonds(O * ψ, SvdCut(; trscheme = final_trscheme); normalize = false) got_one_sweep = approximate((O, ψ), Zipup(; trscheme = final_trscheme)) - alg_zipup = SvdCut(; trscheme = zipup_trscheme) - alg_finalize = SvdCut(; trscheme = final_trscheme) - got_two_sweep = approximate((O, ψ), Zipup(alg_zipup, alg_finalize)) + alg_zipup = Zipup(; trscheme = zipup_trscheme).alg_zipup + alg_zipdown = Zipup(; trscheme = final_trscheme).alg_zipup + got_two_sweep = approximate((O, ψ), Zipup(alg_zipup, alg_zipdown)) err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) From e84aa7e99750bf5b9e977609188ce6cbb99dd622 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 29 Jul 2026 16:47:59 +0800 Subject: [PATCH 11/21] Change trscheme to trunc --- src/algorithms/approximate/zipup.jl | 22 +++++++++++++++------- test/algorithms/zipup.jl | 20 +++++++++----------- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index f388184b9..b7a0f63d9 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -11,12 +11,14 @@ $(TYPEDFIELDS) ## Constructors - Zipup(; trscheme, alg_svd=Defaults.alg_svd()) + Zipup(; trunc, alg_svd=Defaults.alg_svd()) Zipup(alg_zipup, [alg_zipdown]) Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation -scheme and singular value decomposition algorithm. If `alg_zipdown` is provided, the state -obtained after the zip-up sweep is further compressed by sweeping back from right to left. +scheme and singular value decomposition algorithm. The keyword `trunc` can be either one +truncation strategy for a single zip-up sweep, or a tuple `(zipup_trunc, zipdown_trunc)` for a +zip-up sweep followed by a zip-down sweep. Equivalently, one can pass the corresponding truncated +gauge algorithms directly as `alg_zipup` and `alg_zipdown`. Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` @@ -39,10 +41,16 @@ end Zipup(alg_zipup) = Zipup(alg_zipup, nothing) -function Zipup(; - trscheme::TruncationStrategy, alg_svd = Defaults.alg_svd() - ) - return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)) +function Zipup(; trunc, alg_svd = Defaults.alg_svd()) + if trunc isa TruncationStrategy + return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc)) + elseif trunc isa Tuple{<:TruncationStrategy, <:TruncationStrategy} + alg_zipup = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[1]) + alg_zipdown = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[2]) + return Zipup(alg_zipup, alg_zipdown) + else + throw(ArgumentError("`trunc` should be a truncation strategy or a tuple of two truncation strategies")) + end end function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) diff --git a/test/algorithms/zipup.jl b/test/algorithms/zipup.jl index a4e515835..dbbdc8e3d 100644 --- a/test/algorithms/zipup.jl +++ b/test/algorithms/zipup.jl @@ -33,9 +33,9 @@ end O_copy = copy(O) ψ_copy = copy(ψ) - trscheme = trunctol(; atol = 1.0e-10) - ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false) - got = approximate((O, ψ), Zipup(; trscheme)) + trunc = trunctol(; atol = 1.0e-10) + ref = changebonds(O * ψ, SvdCut(; trscheme = trunc); normalize = false) + got = approximate((O, ψ), Zipup(; trunc)) @test norm(ref - got) / norm(ref) < 1.0e-10 @test norm(ψ - ψ_copy) < 1.0e-12 @@ -45,14 +45,12 @@ end @testset "Paeckel two-stage zip-up $(spacetype(pspace))" for (pspace, Dspace, Dcut) in spacelist O, ψ = _random_mpo_mps(pspace, Dspace) rtol = 1.0e-8 - final_trscheme = truncrank(Dcut) & truncerror(; rtol) - zipup_trscheme = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) - - ref_tr = changebonds(O * ψ, SvdCut(; trscheme = final_trscheme); normalize = false) - got_one_sweep = approximate((O, ψ), Zipup(; trscheme = final_trscheme)) - alg_zipup = Zipup(; trscheme = zipup_trscheme).alg_zipup - alg_zipdown = Zipup(; trscheme = final_trscheme).alg_zipup - got_two_sweep = approximate((O, ψ), Zipup(alg_zipup, alg_zipdown)) + final_trunc = truncrank(Dcut) & truncerror(; rtol) + zipup_trunc = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) + + ref_tr = changebonds(O * ψ, SvdCut(; trscheme = final_trunc); normalize = false) + got_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc)) + got_two_sweep = approximate((O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc))) err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) From 3529fc66180bd5e40d014704c08fa0a168c3a63f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 29 Jul 2026 17:19:26 +0800 Subject: [PATCH 12/21] Update `approximate` docstring --- src/algorithms/approximate/approximate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index e7716208f..502f62fe0 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -42,7 +42,7 @@ infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` ha |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:| | `DMRG` | single-site, fixes bond dim | `AbstractFiniteMPS` | ✅ | ✅ | | `DMRG2` | two-site, truncates via `trscheme` | `AbstractFiniteMPS` | ✅ | ✅ | -| `Zipup` | left-to-right streaming MPO-MPS compression | none, uses `(O, ψ)` directly | ❌ | ❌ | +| `Zipup` | streaming MPO-MPS compression | None | ❌ (tuple only) | ❌ (out-of-place only) | | `IDMRG` | single-site, thermodynamic limit | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `IDMRG2` | two-site, thermodynamic limit, needs unit cell ≥ 2 | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `VOMPS` | tangent-space truncation | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ❌ (out-of-place only) | From 8aa1fad02c975dc59cc040a5330d0476667ed8d0 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 30 Jul 2026 16:16:26 +0800 Subject: [PATCH 13/21] Stop requiring left/right-most virtual space to be trivial --- src/algorithms/approximate/zipup.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index b7a0f63d9..88a221282 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -55,13 +55,8 @@ end function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) N = check_length(O, ψ) - if !isunitspace(left_virtualspace(O, 1)) || !isunitspace(right_virtualspace(O, N)) - throw(ArgumentError("Zipup is only implemented for open-boundary MPOs")) - end - T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) A = TensorKit.similarstoragetype(eltype(ψ), T) - Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) local carry From 38beddf1a0222e8b2927abf671b0da4a3d818110 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:05:07 -0400 Subject: [PATCH 14/21] Add `_input_physicalspace` for MPO tensors `physicalspace` of an MPO tensor returns its output space, which is the wrong one to compare against a state that the MPO is applied to whenever the input and output physical spaces differ. Co-Authored-By: Claude Opus 5 (1M context) --- src/states/abstractmps.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 71e7fa81c..0b03e85b4 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -169,6 +169,11 @@ physicalspace(O::MPOTensor) = space(O, 2) physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = only(space(O, 2)) physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ)) +# the input physical space of an MPO tensor, i.e. the space of the state it can be applied to. +# This need not equal its output physical space `physicalspace`. +_input_physicalspace(O::MPOTensor) = dual(space(O, 3)) +_input_physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = dual(only(space(O, 3))) + """ eachsite(state::AbstractMPS) From 94a94c052b5bdd9947540d0bc4052c48873440d9 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:05:18 -0400 Subject: [PATCH 15/21] Add partially fused MPO-MPS contractions `_fuse_mpo_mps_left` and `_fuse_mpo_mps_right` leave the virtual legs on one side of the local MPO-MPS product unfused, such that a subsequent factorization across that side supplies the fuser of the next site. Co-Authored-By: Claude Opus 5 (1M context) --- src/operators/mpo.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index ddda45ff1..d4d2d83b7 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -283,6 +283,22 @@ function _fuse_mpo_mps(O::MPOTensor, A::MPSTensor, Fₗ, Fᵣ) return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end +# Variants that leave the virtual legs on one side unfused, for zip-up sweeps where the subsequent +# factorization supplies the fuser of the next site. The unfused legs are ordered as in +# `fuser(_, virtualspace(mps, i), virtualspace(mpo, i))`, such that the factor that spans them can be +# used as the fuser on the next site: `Fₗ` is oriented as `Vfused ← (Vmps ⊗ Vmpo)`, `Fᵣ` as +# `(Vmps ⊗ Vmpo) ← Vfused`. +function _fuse_mpo_mps_left(O::MPOTensor, A::MPSTensor, Fₗ) + @plansor A′[-1 -2; -3 -4] := Fₗ[-1; 1 3] * A[1 2; -3] * O[3 -2; 2 -4] + return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ +end +# note the index order of the physical leg `-3`: it is the last one, such that the resulting tensor is +# partitioned across the new bond and `_transpose_front` of the right factor is again an MPS tensor +function _fuse_mpo_mps_right(O::MPOTensor, A::MPSTensor, Fᵣ) + @plansor A′[-1 -2; -4 -3] := A[-1 1; 2] * O[-2 -3; 1 3] * Fᵣ[2 3; -4] + return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ +end + function Base.:*(mpo::FiniteMPO{<:MPOTensor}, x::AbstractTensorMap) @assert length(mpo) > 1 @assert numout(x) == length(mpo) From 44499990c532ad7fd594206fa08471770af9e13f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:07:19 -0400 Subject: [PATCH 16/21] Make `Zipup` in-place and drop the redundant right fuser MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `approximate!(ψ, (O, ϕ), alg::Zipup)` now writes the product into `ψ`, which is a destination rather than an initial guess and may alias `ϕ`; the out-of-place method allocates it as `similar(ϕ, T)` with the promoted scalar type. Both return `(ψ, ϵ)`, reporting the truncation error that was previously discarded. The zip-up sweep no longer fuses the right virtual bond of every site: after multiplying with the left fuser, the right factor of the truncated decomposition is simultaneously the truncation carry and the fuser of the next site, so the object with both bonds fused is never constructed. The sweep is split into in-place `zipup!` and `zipdown!` helpers, and all compatibility checks happen in the `approximate!` entry point. Co-Authored-By: Claude Opus 5 (1M context) --- src/algorithms/approximate/approximate.jl | 9 +- src/algorithms/approximate/zipup.jl | 125 +++++++++++++++------- 2 files changed, 90 insertions(+), 44 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 15f29e4e5..1c699397f 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -4,7 +4,8 @@ approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ) approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) - approximate((O, ψ), algorithm) -> ψ′ + approximate((O, ψ), algorithm) -> (ψ′, ϵ) + approximate!(ψ₀, (O, ψ), algorithm) -> (ψ, ϵ) Compute an approximation to the application of an operator `O` to the state `ψ` in the form of an MPS, using initial guess `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, @@ -36,13 +37,15 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). Each algorithm below only supports a subset of the general interface. Check this table before picking one — in particular, note that **only `DMRG`/`DMRG2` accept a bare state `ψ`**; the infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` has no in-place -`approximate!`** at all. +`approximate!`** at all. `Zipup` is a single sweep rather than an iterative optimization, so it uses +no environments and returns `(ψ, ϵ)`; its `ψ₀` is a write destination, not an initial guess, and it +may be omitted. | Algorithm | Scheme | State `ψ₀` | bare `ψ` allowed? | `approximate!` | |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:| | `DMRG` | single-site, fixes bond dim | `AbstractFiniteMPS` | ✅ | ✅ | | `DMRG2` | two-site, truncates via `trunc` | `AbstractFiniteMPS` | ✅ | ✅ | -| `Zipup` | streaming MPO-MPS compression | None | ❌ (tuple only) | ❌ (out-of-place only) | +| `Zipup` | streaming MPO-MPS compression | `FiniteMPS` destination, optional | ❌ (tuple only) | ✅ | | `IDMRG` | single-site, thermodynamic limit | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `IDMRG2` | two-site, thermodynamic limit, needs unit cell ≥ 2 | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ | | `VOMPS` | tangent-space truncation | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ❌ (out-of-place only) | diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index 88a221282..2a66c636c 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -5,6 +5,17 @@ Algorithm that approximates an open-boundary finite MPO-MPS product using a left zip-up sweep, optionally followed by a right-to-left zip-down sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual bond is truncated immediately. + approximate((O, ϕ), alg::Zipup) -> ψ, ϵ + approximate!(ψ, (O, ϕ), alg::Zipup) -> ψ, ϵ + +Contrary to the variational algorithms, this algorithm requires no initial guess: the in-place +version simply uses `ψ` as the destination of the sweep, overwriting its contents, and may alias +`ϕ`. The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`. +Both return the truncation error `ϵ` alongside the approximated state. + +Only the input physical spaces of `O` have to match those of `ϕ`: for an `O` whose output physical +spaces differ, the result is a state with the output physical spaces of `O`. + ## Fields $(TYPEDFIELDS) @@ -53,54 +64,86 @@ function Zipup(; trunc, alg_svd = Defaults.alg_svd()) end end -function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) - N = check_length(O, ψ) - T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ)) - A = TensorKit.similarstoragetype(eltype(ψ), T) - Fₗ = fuser(A, left_virtualspace(ψ, 1), left_virtualspace(O, 1)) - local carry - - As = map(1:N) do i - Aψ = i == 1 ? ψ.AC[1] : ψ.AR[i] - physicalspace(Aψ) == physicalspace(O[i]) || +function approximate!(ψ::FiniteMPS, (O, ϕ)::Tuple{Any, <:FiniteMPS}, alg::Zipup) + N = check_length(ψ, O, ϕ) + T = TensorOperations.promote_contract(scalartype(O), scalartype(ϕ)) + promote_type(T, scalartype(ψ)) === scalartype(ψ) || + throw(ArgumentError("destination state with scalartype $(scalartype(ψ)) cannot hold the result with scalartype $T")) + for i in 1:N + physicalspace(ϕ, i) == _input_physicalspace(O[i]) || throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) - Fᵣ = fuser(A, right_virtualspace(ψ, i), right_virtualspace(O, i)) - Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ) - i > 1 && (Aᶻ = _mul_front(carry, Aᶻ)) - - if i == N - return Aᶻ - else - AL, C, _ = left_gauge(Aᶻ, alg.alg_zipup) - carry = C - Fₗ = Fᵣ - return AL - end end - return isnothing(alg.alg_zipdown) ? - FiniteMPS(As; normalize = false, overwrite = true) : - _zipdown(As, alg.alg_zipdown) + ψ, ϵ = zipup!(ψ, O, ϕ, alg.alg_zipup) + if !isnothing(alg.alg_zipdown) + ψ, ϵ′ = zipdown!(ψ, alg.alg_zipdown) + ϵ = max(ϵ, ϵ′) + end + return ψ, ϵ end -function _zipdown(As::Vector{A}, alg::MatrixAlgebraKit.TruncatedAlgorithm) where {A} - N = length(As) - N == 1 && return FiniteMPS(As; normalize = false, overwrite = true) +function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) + O, ϕ = Oϕ + T = TensorOperations.promote_contract(scalartype(O), scalartype(ϕ)) + return approximate!(similar(ϕ, T), Oϕ, alg) +end + +""" + zipup!(ψ, O, ϕ, alg) -> ψ, ϵ + +Contract the MPO `O` with the MPS `ϕ` in a single left-to-right sweep, truncating the enlarged +virtual bond at every site with `alg`, and write the result into `ψ`. The destination is left with +its gauge center on the last site, and may alias `ϕ`. - ARs = Vector{Union{Missing, A}}(missing, N) - ALs = Vector{Union{Missing, A}}(missing, N) - ACs = Vector{Union{Missing, A}}(missing, N) +Instead of fusing both virtual bonds of every site, only the left bond is fused: the right factor of +the truncated decomposition is simultaneously the truncation carry and the fuser of the next site, so +the enlarged object is never constructed. - local C - AC = As[N] - for i in N:-1:2 - C, AR, _ = right_gauge(AC, alg) - ARs[i] = AR - AC = _mul_tail(As[i - 1], C) +Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all +bonds. + +Compatibility of the lengths, physical spaces and scalar types of `ψ`, `O` and `ϕ` is assumed, and +checked in [`approximate!`](@ref). +""" +function zipup!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg) + N = length(ψ) + + # obtain all input tensors before overwriting the destination, such that `ψ === ϕ` is allowed: + # from here on, the input is only queried through `Aϕs`, never through `ϕ` itself + Aϕs = map(i -> i == 1 ? ϕ.AC[1] : ϕ.AR[i], 1:N) + + # the sweep re-derives the entire state: discard all cached tensors, as their spaces are stale + # TODO: "reallocate" tensors?" + foreach(f -> fill!(f, missing), (ψ.ALs, ψ.ARs, ψ.ACs, ψ.Cs)) + + A = storagetype(eltype(ψ)) + Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1)) + ϵ = zero(real(scalartype(ψ))) + for i in 1:(N - 1) + Aᶻ = _fuse_mpo_mps_left(O[i], Aϕs[i], Fₗ) + AL, Fₗ, ϵᵢ = left_gauge(Aᶻ, alg) # right factor doubles as the next left fuser + ψ.ALs[i] = AL + ϵ = max(ϵ, ϵᵢ) end + Fᵣ = fuser(A, right_virtualspace(Aϕs[N]), right_virtualspace(O, N)) + ψ.ACs[N] = _fuse_mpo_mps(O[N], Aϕs[N], Fₗ, Fᵣ) - B = typeof(C) - ACs[1] = AC - Cs = Vector{Union{Missing, B}}(missing, N + 1) - return FiniteMPS(ALs, ARs, ACs, Cs) + return ψ, ϵ +end + +""" + zipdown!(ψ, alg) -> ψ, ϵ + +Sweep `ψ` from right to left, truncating every bond with `alg` in a locally gauged basis, and moving +the gauge center to the leftmost bond in the process. + +Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all bonds. +""" +function zipdown!(ψ::AbstractFiniteMPS, alg) + ϵ = zero(real(scalartype(ψ))) + for i in length(ψ):-1:2 + ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg) + ϵ = max(ϵ, ϵᵢ) + end + return ψ, ϵ end From 5dce6cc82d37cba36aad72ce973e5ebf0383c6ea Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:07:31 -0400 Subject: [PATCH 17/21] Merge the zip sweeps and make the direction selectable `zipup!` and `zipdown!` are merged into `zip_left_right!`, which zips up from left to right and then truncates from right to left, and its mirror image `zip_right_left!`. The new `left_to_right` field of `Zipup` selects between them. The right-to-left sweep needs the mirrored partially fused contraction, whose local tensors are already partitioned across the new bond, so the leg permutation that `right_gauge` applies to MPS tensors has to be skipped. Co-Authored-By: Claude Opus 5 (1M context) --- src/algorithms/approximate/zipup.jl | 137 +++++++++++++++++----------- 1 file changed, 85 insertions(+), 52 deletions(-) diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index 2a66c636c..5ce3ee04d 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -1,39 +1,34 @@ """ $(TYPEDEF) -Algorithm that approximates an open-boundary finite MPO-MPS product using a left-to-right -zip-up sweep, optionally followed by a right-to-left zip-down sweep. The MPO and MPS are -contracted one site at a time, and the enlarged virtual bond is truncated immediately. +Algorithm that approximates an open-boundary finite MPO-MPS product using a zip-up sweep, optionally +followed by a zip-down sweep in the opposite direction. The MPO and MPS are contracted one site at a +time, and the enlarged virtual bond is truncated immediately. The sweep direction is selected by +`left_to_right`. approximate((O, ϕ), alg::Zipup) -> ψ, ϵ approximate!(ψ, (O, ϕ), alg::Zipup) -> ψ, ϵ -Contrary to the variational algorithms, this algorithm requires no initial guess: the in-place -version simply uses `ψ` as the destination of the sweep, overwriting its contents, and may alias -`ϕ`. The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`. +Contrary to the variational algorithms, this algorithm requires no initial guess: +the in-place version simply uses `ψ` as the destination of the sweep, overwriting its contents, and may alias `ϕ`. +The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`. Both return the truncation error `ϵ` alongside the approximated state. -Only the input physical spaces of `O` have to match those of `ϕ`: for an `O` whose output physical -spaces differ, the result is a state with the output physical spaces of `O`. - ## Fields $(TYPEDFIELDS) ## Constructors - Zipup(; trunc, alg_svd=Defaults.alg_svd()) - Zipup(alg_zipup, [alg_zipdown]) + Zipup(; trunc, alg_svd=Defaults.alg_svd(), left_to_right=true) + Zipup(alg_zipup, [alg_zipdown]; left_to_right=true) -Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation -scheme and singular value decomposition algorithm. The keyword `trunc` can be either one -truncation strategy for a single zip-up sweep, or a tuple `(zipup_trunc, zipdown_trunc)` for a -zip-up sweep followed by a zip-down sweep. Equivalently, one can pass the corresponding truncated -gauge algorithms directly as `alg_zipup` and `alg_zipdown`. +Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation scheme and singular value decomposition algorithm. +The keyword `trunc` can be either one truncation strategy for a single zip-up sweep, or a tuple `(zipup_trunc, zipdown_trunc)` for a zip-up sweep followed by a zip-down sweep. +Equivalently, one can pass the corresponding truncated gauge algorithms directly as `alg_zipup` and `alg_zipdown`. +The keyword `left_to_right` selects the direction of the zip-up sweep, the zip-down sweep always running in the opposite direction. -Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more -permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` -to impose the final truncation. +Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` to impose the final truncation. ## References @@ -48,17 +43,21 @@ struct Zipup{ alg_zipup::U "algorithm used for the final locally gauged truncation pass; `nothing` skips this pass" alg_zipdown::D + "if `true`, zip up from left to right and truncate from right to left, and vice versa" + left_to_right::Bool end -Zipup(alg_zipup) = Zipup(alg_zipup, nothing) +function Zipup(alg_zipup, alg_zipdown = nothing; left_to_right::Bool = true) + return Zipup(alg_zipup, alg_zipdown, left_to_right) +end -function Zipup(; trunc, alg_svd = Defaults.alg_svd()) +function Zipup(; trunc, alg_svd = Defaults.alg_svd(), left_to_right::Bool = true) if trunc isa TruncationStrategy - return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc)) + return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc); left_to_right) elseif trunc isa Tuple{<:TruncationStrategy, <:TruncationStrategy} alg_zipup = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[1]) alg_zipdown = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[2]) - return Zipup(alg_zipup, alg_zipdown) + return Zipup(alg_zipup, alg_zipdown; left_to_right) else throw(ArgumentError("`trunc` should be a truncation strategy or a tuple of two truncation strategies")) end @@ -74,12 +73,11 @@ function approximate!(ψ::FiniteMPS, (O, ϕ)::Tuple{Any, <:FiniteMPS}, alg::Zipu throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) end - ψ, ϵ = zipup!(ψ, O, ϕ, alg.alg_zipup) - if !isnothing(alg.alg_zipdown) - ψ, ϵ′ = zipdown!(ψ, alg.alg_zipdown) - ϵ = max(ϵ, ϵ′) + return if alg.left_to_right + zip_left_right!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown) + else + zip_right_left!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown) end - return ψ, ϵ end function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) @@ -88,24 +86,23 @@ function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) return approximate!(similar(ϕ, T), Oϕ, alg) end -""" - zipup!(ψ, O, ϕ, alg) -> ψ, ϵ - -Contract the MPO `O` with the MPS `ϕ` in a single left-to-right sweep, truncating the enlarged -virtual bond at every site with `alg`, and write the result into `ψ`. The destination is left with -its gauge center on the last site, and may alias `ϕ`. +@doc """ + zip_left_right!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ + zip_right_left!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ -Instead of fusing both virtual bonds of every site, only the left bond is fused: the right factor of -the truncated decomposition is simultaneously the truncation carry and the fuser of the next site, so -the enlarged object is never constructed. +Contract the MPO `O` with the MPS `ϕ` in a single sweep, truncating the enlarged virtual bond at every +site with `alg_zipup`, and write the result into `ψ`. `zip_left_right!` zips up from left to right, +`zip_right_left!` from right to left. Unless `alg_zipdown` is `nothing`, a second sweep in the +opposite direction imposes a final truncation with `alg_zipdown` in a locally gauged basis, leaving +the gauge center of `ψ` at the far end. The destination may alias `ϕ`. Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all -bonds. - -Compatibility of the lengths, physical spaces and scalar types of `ψ`, `O` and `ϕ` is assumed, and -checked in [`approximate!`](@ref). +bonds and both sweeps. """ -function zipup!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg) +zip_left_right! +@doc (@doc zip_left_right!) zip_right_left! + +function zip_left_right!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing) N = length(ψ) # obtain all input tensors before overwriting the destination, such that `ψ === ϕ` is allowed: @@ -119,31 +116,67 @@ function zipup!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg) A = storagetype(eltype(ψ)) Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1)) ϵ = zero(real(scalartype(ψ))) + + # zip up from left to right, leaving the gauge center on the last site for i in 1:(N - 1) Aᶻ = _fuse_mpo_mps_left(O[i], Aϕs[i], Fₗ) - AL, Fₗ, ϵᵢ = left_gauge(Aᶻ, alg) # right factor doubles as the next left fuser + AL, Fₗ, ϵᵢ = left_gauge(Aᶻ, alg_zipup) # right factor doubles as the next left fuser ψ.ALs[i] = AL ϵ = max(ϵ, ϵᵢ) end Fᵣ = fuser(A, right_virtualspace(Aϕs[N]), right_virtualspace(O, N)) ψ.ACs[N] = _fuse_mpo_mps(O[N], Aϕs[N], Fₗ, Fᵣ) + # zip down from right to left, truncating in a locally gauged basis + if !isnothing(alg_zipdown) + for i in N:-1:2 + ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg_zipdown) + ϵ = max(ϵ, ϵᵢ) + end + end + return ψ, ϵ end -""" - zipdown!(ψ, alg) -> ψ, ϵ +function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing) + N = length(ψ) -Sweep `ψ` from right to left, truncating every bond with `alg` in a locally gauged basis, and moving -the gauge center to the leftmost bond in the process. + Aϕs = map(i -> i == N ? ϕ.AC[N] : ϕ.AL[i], 1:N) + foreach(f -> fill!(f, missing), (ψ.ALs, ψ.ARs, ψ.ACs, ψ.Cs)) -Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all bonds. -""" -function zipdown!(ψ::AbstractFiniteMPS, alg) + A = storagetype(eltype(ψ)) + # the right-hand fusers are oriented as `(Vmps ⊗ Vmpo) ← Vfused`, matching the factor that + # replaces them on the next site + Vᵣ = right_virtualspace(Aϕs[N]) ⊗ right_virtualspace(O, N) + Fᵣ = isomorphism(A, Vᵣ, fuse(Vᵣ)) ϵ = zero(real(scalartype(ψ))) - for i in length(ψ):-1:2 - ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg) + + # zip up from right to left, leaving the gauge center on the first site + for i in N:-1:2 + Aᶻ = _fuse_mpo_mps_right(O[i], Aϕs[i], Fᵣ) + Fᵣ, AR, ϵᵢ = _right_gauge_zip(Aᶻ, alg_zipup) # left factor doubles as the next right fuser + ψ.ARs[i] = AR ϵ = max(ϵ, ϵᵢ) end + # the carry is oriented such that it can simply be composed with the last local tensor + Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1)) + ψ.ACs[1] = _fuse_mpo_mps_left(O[1], Aϕs[1], Fₗ) * Fᵣ + + # zip down from left to right, truncating in a locally gauged basis + if !isnothing(alg_zipdown) + for i in 1:(N - 1) + ψ, ϵᵢ = left_gauge!(ψ, i, ψ.AC[i], alg_zipdown) + ϵ = max(ϵ, ϵᵢ) + end + end + return ψ, ϵ end + +# `right_gauge` for the tensors of a right-to-left zip-up sweep: these are already partitioned across +# the new bond, so the leg permutation that `right_gauge` applies to MPS tensors has to be skipped +function _right_gauge_zip(Aᶻ, alg::MatrixAlgebraKit.TruncatedAlgorithm) + U, S, Vᴴ, ϵ = svd_trunc(Aᶻ, alg) + C = LinearAlgebra.rmul!(U, S) # C = U * S, matching `RightOrthViaSVD` + return C, _transpose_front(Vᴴ), ϵ +end From 0249b1f9642e5e64cad059d135e68dd44230802d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:07:45 -0400 Subject: [PATCH 18/21] test: cover the in-place and directional `Zipup` interface Both sweep directions, with and without the zip-down pass, are checked against the dense MPO-MPS product, along with an empty destination, a destination with unrelated contents, a destination aliasing the input, scalar type promotion, non-trivial boundary virtual spaces, and an MPO whose input and output physical spaces differ. Co-Authored-By: Claude Opus 5 (1M context) --- test/algorithms/zipup.jl | 107 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 100 insertions(+), 7 deletions(-) diff --git a/test/algorithms/zipup.jl b/test/algorithms/zipup.jl index 309063abb..bd8b5479f 100644 --- a/test/algorithms/zipup.jl +++ b/test/algorithms/zipup.jl @@ -16,7 +16,7 @@ spacelist = [ (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8), ] -function _random_mpo_mps(pspace, Dspace) +function _random_mpo_mps(pspace, Dspace; elt = ComplexF64) Random.seed!(1357) L = 6 Wspace = Dspace @@ -24,7 +24,7 @@ function _random_mpo_mps(pspace, Dspace) O = FiniteMPO( [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] ) - ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) + ψ = FiniteMPS(rand, elt, L, pspace, Dspace) return O, ψ end @@ -35,22 +35,31 @@ end trunc = trunctol(; atol = 1.0e-10) ref = changebonds(O * ψ, SvdCut(; trunc); normalize = false) - got = approximate((O, ψ), Zipup(; trunc)) - @test norm(ref - got) / norm(ref) < 1.0e-10 + # both sweep directions, with and without the zip-down pass + for left_to_right in (true, false), trunc′ in (trunc, (notrunc(), trunc)) + got, ϵ = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) + @test norm(ref - got) / norm(ref) < 1.0e-10 + @test ϵ < 1.0e-10 + end + @test norm(ψ - ψ_copy) < 1.0e-12 @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) end -@testset "Paeckel two-stage zip-up $(spacetype(pspace))" for (pspace, Dspace, Dcut) in spacelist +@testset "Paeckel two-stage zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for + (pspace, Dspace, Dcut) in spacelist, left_to_right in (true, false) O, ψ = _random_mpo_mps(pspace, Dspace) rtol = 1.0e-8 final_trunc = truncrank(Dcut) & truncerror(; rtol) zipup_trunc = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) ref_tr = changebonds(O * ψ, SvdCut(; trunc = final_trunc); normalize = false) - got_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc)) - got_two_sweep = approximate((O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc))) + got_one_sweep, ϵ_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc, left_to_right)) + got_two_sweep, _ = approximate( + (O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc), left_to_right) + ) + @test ϵ_one_sweep > 0 err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) @@ -58,3 +67,87 @@ end @test maximum(i -> dim(left_virtualspace(got_one_sweep, i)), 2:length(got_one_sweep)) <= Dcut @test maximum(i -> dim(left_virtualspace(got_two_sweep, i)), 2:length(got_two_sweep)) <= Dcut end + +@testset "In-place zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for + (pspace, Dspace, Dcut) in spacelist, left_to_right in (true, false) + O, ψ = _random_mpo_mps(pspace, Dspace) + alg = Zipup(; trunc = (truncrank(2Dcut), truncrank(Dcut)), left_to_right) + ref, ϵ_ref = approximate((O, ψ), alg) + + # empty destination + dst = similar(ψ, ComplexF64) + got, ϵ = approximate!(dst, (O, ψ), alg) + @test got === dst + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref + + # a destination with unrelated contents is overwritten entirely + dst = FiniteMPS(rand, ComplexF64, length(ψ), pspace, oneunit(Dspace) ⊕ Dspace ⊕ Dspace) + got, ϵ = approximate!(dst, (O, ψ), alg) + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref + + # the input may serve as its own destination + got, ϵ = approximate!(ψ, (O, ψ), alg) + @test got === ψ + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref +end + +@testset "Zip-up with non-trivial boundary spaces $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist + Random.seed!(1357) + L = 4 + Vspaces = fill(Dspace, L + 1) + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace; left = Dspace, right = Dspace) + + trunc = trunctol(; atol = 1.0e-12) + got, _ = approximate((O, ψ), Zipup(; trunc)) + got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) + + # the boundary virtual spaces of the product are the fused ones, in both variants + for ψ′ in (got, got_two_sweep) + @test left_virtualspace(ψ′, 1) == fuse(left_virtualspace(ψ, 1) ⊗ left_virtualspace(O, 1)) + @test right_virtualspace(ψ′, L) == fuse(right_virtualspace(ψ, L) ⊗ right_virtualspace(O, L)) + end + @test norm(got - got_two_sweep) / norm(got) < 1.0e-10 + @test norm(got) ≈ norm(O * ψ) +end + +@testset "Zip-up with non-square MPO $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist + Random.seed!(1357) + L = 4 + pspace′ = pspace ⊕ oneunit(pspace) # output physical space, different from the input one + Vspaces = [oneunit(Dspace); fill(Dspace, L - 1); oneunit(Dspace)] + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace′ ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) + + trunc = trunctol(; atol = 1.0e-12) + ref = O * ψ + got, _ = approximate((O, ψ), Zipup(; trunc)) + got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) + + @test all(i -> physicalspace(got, i) == pspace′, 1:L) + @test norm(ref - got) / norm(ref) < 1.0e-10 + @test norm(ref - got_two_sweep) / norm(ref) < 1.0e-10 + + # the reverse product is not defined + @test_throws SpaceMismatch approximate((O, got), Zipup(; trunc)) +end + +@testset "Zip-up scalar type promotion $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist + O, ψ = _random_mpo_mps(pspace, Dspace; elt = Float64) + alg = Zipup(; trunc = trunctol(; atol = 1.0e-10)) + + got, _ = approximate((O, ψ), alg) + @test scalartype(got) === ComplexF64 + ref = changebonds(O * ψ, SvdCut(; trunc = trunctol(; atol = 1.0e-10)); normalize = false) + @test norm(ref - got) / norm(ref) < 1.0e-10 + + # a real destination cannot hold the complex result + @test_throws ArgumentError approximate!(ψ, (O, ψ), alg) +end From 0e81bd2f74b74f5b693cd764974dc0638f8fcacd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:07:46 -0400 Subject: [PATCH 19/21] Add changelog entry for `Zipup` Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/changelog.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 21b5caddc..a920c6437 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -23,6 +23,11 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Addition of `FiniteMPS`/`FiniteMPO` with different scalar types, through a new `Base.similar(ψ, ::Type{S})` for `S <: Number` on `FiniteMPS`. ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484)) +- `Zipup`, an algorithm for `approximate`/`approximate!` that compresses a finite MPO-MPS product in + a single sweep, optionally followed by a sweep in the opposite direction that imposes the final + truncation. The sweep direction is selected by the `left_to_right` keyword. Both + `approximate((O, ϕ), alg)` and `approximate!(ψ, (O, ϕ), alg)` are supported, where the destination + `ψ` is a write target rather than an initial guess and may alias `ϕ`; they return `(ψ, ϵ)`. ### Changed From 7e9ebb559b094cfa990a3a99674cebfed2cbf887 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 30 Jul 2026 11:14:28 -0400 Subject: [PATCH 20/21] test: fold the `Zipup` tests into the `approximate` testset The zip-up sweeps are `approximate` algorithms, so their tests belong with the other ones rather than in a file of their own. Co-Authored-By: Claude Opus 5 (1M context) --- test/algorithms/approximate.jl | 143 ++++++++++++++++++++++++++++++ test/algorithms/zipup.jl | 153 --------------------------------- 2 files changed, 143 insertions(+), 153 deletions(-) delete mode 100644 test/algorithms/zipup.jl diff --git a/test/algorithms/approximate.jl b/test/algorithms/approximate.jl index 6c3b046f3..7af88ee7c 100644 --- a/test/algorithms/approximate.jl +++ b/test/algorithms/approximate.jl @@ -9,9 +9,28 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ +using Random verbosity_conv = 1 +# fixtures for the `Zipup` testsets +zipup_spacelist = [ + (ℙ^4, ℙ^3, 4), + (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8), +] + +function _random_mpo_mps(pspace, Dspace; elt = ComplexF64) + Random.seed!(1357) + L = 6 + Wspace = Dspace + Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)] + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, elt, L, pspace, Dspace) + return O, ψ +end + @testset "approximate" verbose = true begin verbosity = verbosity_conv @testset "mpo * infinite ≈ infinite" begin @@ -84,4 +103,128 @@ verbosity_conv = 1 @test norm(O * ψ₁ - ψ₂) ≈ 0 atol = 0.001 end + + @testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist + O, ψ = _random_mpo_mps(pspace, Dspace) + O_copy = copy(O) + ψ_copy = copy(ψ) + + trunc = trunctol(; atol = 1.0e-10) + ref = changebonds(O * ψ, SvdCut(; trunc); normalize = false) + + # both sweep directions, with and without the zip-down pass + for left_to_right in (true, false), trunc′ in (trunc, (notrunc(), trunc)) + got, ϵ = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) + @test norm(ref - got) / norm(ref) < 1.0e-10 + @test ϵ < 1.0e-10 + end + + @test norm(ψ - ψ_copy) < 1.0e-12 + @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) + end + + @testset "Paeckel two-stage zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for + (pspace, Dspace, Dcut) in zipup_spacelist, left_to_right in (true, false) + O, ψ = _random_mpo_mps(pspace, Dspace) + rtol = 1.0e-8 + final_trunc = truncrank(Dcut) & truncerror(; rtol) + zipup_trunc = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) + + ref_tr = changebonds(O * ψ, SvdCut(; trunc = final_trunc); normalize = false) + got_one_sweep, ϵ_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc, left_to_right)) + got_two_sweep, _ = approximate( + (O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc), left_to_right) + ) + @test ϵ_one_sweep > 0 + + err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) + err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) + @test err_two_sweep < err_one_sweep / 2 + @test maximum(i -> dim(left_virtualspace(got_one_sweep, i)), 2:length(got_one_sweep)) <= Dcut + @test maximum(i -> dim(left_virtualspace(got_two_sweep, i)), 2:length(got_two_sweep)) <= Dcut + end + + @testset "In-place zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for + (pspace, Dspace, Dcut) in zipup_spacelist, left_to_right in (true, false) + O, ψ = _random_mpo_mps(pspace, Dspace) + alg = Zipup(; trunc = (truncrank(2Dcut), truncrank(Dcut)), left_to_right) + ref, ϵ_ref = approximate((O, ψ), alg) + + # empty destination + dst = similar(ψ, ComplexF64) + got, ϵ = approximate!(dst, (O, ψ), alg) + @test got === dst + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref + + # a destination with unrelated contents is overwritten entirely + dst = FiniteMPS(rand, ComplexF64, length(ψ), pspace, oneunit(Dspace) ⊕ Dspace ⊕ Dspace) + got, ϵ = approximate!(dst, (O, ψ), alg) + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref + + # the input may serve as its own destination + got, ϵ = approximate!(ψ, (O, ψ), alg) + @test got === ψ + @test norm(ref - got) / norm(ref) < 1.0e-12 + @test ϵ ≈ ϵ_ref + end + + @testset "Zip-up with non-trivial boundary spaces $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist + Random.seed!(1357) + L = 4 + Vspaces = fill(Dspace, L + 1) + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace; left = Dspace, right = Dspace) + + trunc = trunctol(; atol = 1.0e-12) + got, _ = approximate((O, ψ), Zipup(; trunc)) + got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) + + # the boundary virtual spaces of the product are the fused ones, in both variants + for ψ′ in (got, got_two_sweep) + @test left_virtualspace(ψ′, 1) == fuse(left_virtualspace(ψ, 1) ⊗ left_virtualspace(O, 1)) + @test right_virtualspace(ψ′, L) == fuse(right_virtualspace(ψ, L) ⊗ right_virtualspace(O, L)) + end + @test norm(got - got_two_sweep) / norm(got) < 1.0e-10 + @test norm(got) ≈ norm(O * ψ) + end + + @testset "Zip-up with non-square MPO $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist + Random.seed!(1357) + L = 4 + pspace′ = pspace ⊕ oneunit(pspace) # output physical space, different from the input one + Vspaces = [oneunit(Dspace); fill(Dspace, L - 1); oneunit(Dspace)] + O = FiniteMPO( + [rand(ComplexF64, Vspaces[i] ⊗ pspace′ ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] + ) + ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) + + trunc = trunctol(; atol = 1.0e-12) + ref = O * ψ + got, _ = approximate((O, ψ), Zipup(; trunc)) + got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) + + @test all(i -> physicalspace(got, i) == pspace′, 1:L) + @test norm(ref - got) / norm(ref) < 1.0e-10 + @test norm(ref - got_two_sweep) / norm(ref) < 1.0e-10 + + # the reverse product is not defined + @test_throws SpaceMismatch approximate((O, got), Zipup(; trunc)) + end + + @testset "Zip-up scalar type promotion $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist + O, ψ = _random_mpo_mps(pspace, Dspace; elt = Float64) + alg = Zipup(; trunc = trunctol(; atol = 1.0e-10)) + + got, _ = approximate((O, ψ), alg) + @test scalartype(got) === ComplexF64 + ref = changebonds(O * ψ, SvdCut(; trunc = trunctol(; atol = 1.0e-10)); normalize = false) + @test norm(ref - got) / norm(ref) < 1.0e-10 + + # a real destination cannot hold the complex result + @test_throws ArgumentError approximate!(ψ, (O, ψ), alg) + end end diff --git a/test/algorithms/zipup.jl b/test/algorithms/zipup.jl deleted file mode 100644 index bd8b5479f..000000000 --- a/test/algorithms/zipup.jl +++ /dev/null @@ -1,153 +0,0 @@ -println(" ------------------------------ -| Zipup tests | ------------------------------ -") - -using .TestSetup -using Test -using MPSKit -using TensorKit -using TensorKit: ℙ -using Random - -spacelist = [ - (ℙ^4, ℙ^3, 4), - (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8), -] - -function _random_mpo_mps(pspace, Dspace; elt = ComplexF64) - Random.seed!(1357) - L = 6 - Wspace = Dspace - Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)] - O = FiniteMPO( - [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] - ) - ψ = FiniteMPS(rand, elt, L, pspace, Dspace) - return O, ψ -end - -@testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist - O, ψ = _random_mpo_mps(pspace, Dspace) - O_copy = copy(O) - ψ_copy = copy(ψ) - - trunc = trunctol(; atol = 1.0e-10) - ref = changebonds(O * ψ, SvdCut(; trunc); normalize = false) - - # both sweep directions, with and without the zip-down pass - for left_to_right in (true, false), trunc′ in (trunc, (notrunc(), trunc)) - got, ϵ = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) - @test norm(ref - got) / norm(ref) < 1.0e-10 - @test ϵ < 1.0e-10 - end - - @test norm(ψ - ψ_copy) < 1.0e-12 - @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O)) -end - -@testset "Paeckel two-stage zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for - (pspace, Dspace, Dcut) in spacelist, left_to_right in (true, false) - O, ψ = _random_mpo_mps(pspace, Dspace) - rtol = 1.0e-8 - final_trunc = truncrank(Dcut) & truncerror(; rtol) - zipup_trunc = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) - - ref_tr = changebonds(O * ψ, SvdCut(; trunc = final_trunc); normalize = false) - got_one_sweep, ϵ_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc, left_to_right)) - got_two_sweep, _ = approximate( - (O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc), left_to_right) - ) - @test ϵ_one_sweep > 0 - - err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) - err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) - @test err_two_sweep < err_one_sweep / 2 - @test maximum(i -> dim(left_virtualspace(got_one_sweep, i)), 2:length(got_one_sweep)) <= Dcut - @test maximum(i -> dim(left_virtualspace(got_two_sweep, i)), 2:length(got_two_sweep)) <= Dcut -end - -@testset "In-place zip-up $(spacetype(pspace)), left_to_right = $left_to_right" for - (pspace, Dspace, Dcut) in spacelist, left_to_right in (true, false) - O, ψ = _random_mpo_mps(pspace, Dspace) - alg = Zipup(; trunc = (truncrank(2Dcut), truncrank(Dcut)), left_to_right) - ref, ϵ_ref = approximate((O, ψ), alg) - - # empty destination - dst = similar(ψ, ComplexF64) - got, ϵ = approximate!(dst, (O, ψ), alg) - @test got === dst - @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref - - # a destination with unrelated contents is overwritten entirely - dst = FiniteMPS(rand, ComplexF64, length(ψ), pspace, oneunit(Dspace) ⊕ Dspace ⊕ Dspace) - got, ϵ = approximate!(dst, (O, ψ), alg) - @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref - - # the input may serve as its own destination - got, ϵ = approximate!(ψ, (O, ψ), alg) - @test got === ψ - @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref -end - -@testset "Zip-up with non-trivial boundary spaces $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist - Random.seed!(1357) - L = 4 - Vspaces = fill(Dspace, L + 1) - O = FiniteMPO( - [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] - ) - ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace; left = Dspace, right = Dspace) - - trunc = trunctol(; atol = 1.0e-12) - got, _ = approximate((O, ψ), Zipup(; trunc)) - got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) - - # the boundary virtual spaces of the product are the fused ones, in both variants - for ψ′ in (got, got_two_sweep) - @test left_virtualspace(ψ′, 1) == fuse(left_virtualspace(ψ, 1) ⊗ left_virtualspace(O, 1)) - @test right_virtualspace(ψ′, L) == fuse(right_virtualspace(ψ, L) ⊗ right_virtualspace(O, L)) - end - @test norm(got - got_two_sweep) / norm(got) < 1.0e-10 - @test norm(got) ≈ norm(O * ψ) -end - -@testset "Zip-up with non-square MPO $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist - Random.seed!(1357) - L = 4 - pspace′ = pspace ⊕ oneunit(pspace) # output physical space, different from the input one - Vspaces = [oneunit(Dspace); fill(Dspace, L - 1); oneunit(Dspace)] - O = FiniteMPO( - [rand(ComplexF64, Vspaces[i] ⊗ pspace′ ← pspace ⊗ Vspaces[i + 1]) for i in 1:L] - ) - ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) - - trunc = trunctol(; atol = 1.0e-12) - ref = O * ψ - got, _ = approximate((O, ψ), Zipup(; trunc)) - got_two_sweep, _ = approximate((O, ψ), Zipup(; trunc = (notrunc(), trunc))) - - @test all(i -> physicalspace(got, i) == pspace′, 1:L) - @test norm(ref - got) / norm(ref) < 1.0e-10 - @test norm(ref - got_two_sweep) / norm(ref) < 1.0e-10 - - # the reverse product is not defined - @test_throws SpaceMismatch approximate((O, got), Zipup(; trunc)) -end - -@testset "Zip-up scalar type promotion $(spacetype(pspace))" for (pspace, Dspace, _) in spacelist - O, ψ = _random_mpo_mps(pspace, Dspace; elt = Float64) - alg = Zipup(; trunc = trunctol(; atol = 1.0e-10)) - - got, _ = approximate((O, ψ), alg) - @test scalartype(got) === ComplexF64 - ref = changebonds(O * ψ, SvdCut(; trunc = trunctol(; atol = 1.0e-10)); normalize = false) - @test norm(ref - got) / norm(ref) < 1.0e-10 - - # a real destination cannot hold the complex result - @test_throws ArgumentError approximate!(ψ, (O, ψ), alg) -end From 64abd627cb2fd332d2fe23413db613c055824f2b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 31 Jul 2026 12:17:39 +0800 Subject: [PATCH 21/21] Unicode drawing for `fuse_mpo_mps` --- src/operators/mpo.jl | 42 +++++++++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index d4d2d83b7..c5c033fe4 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -278,24 +278,48 @@ function Base.:*(mpo::InfiniteMPO, mps::InfiniteMPS) return changebonds(InfiniteMPS(As), SvdCut(; trunc = notrunc())) end +""" +Fuse the left and right virtual legs of the product of MPO-MPS tensors +``` + ┌---A---┐ + 1 --Fl | Fr-- 3 => A′[1 2; 3] + └---O---┘ + | + 2 +``` +""" function _fuse_mpo_mps(O::MPOTensor, A::MPSTensor, Fₗ, Fᵣ) @plansor A′[-1 -2; -3] := Fₗ[-1; 1 3] * A[1 2; 4] * O[3 -2; 2 5] * conj(Fᵣ[-3; 4 5]) return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end - -# Variants that leave the virtual legs on one side unfused, for zip-up sweeps where the subsequent -# factorization supplies the fuser of the next site. The unfused legs are ordered as in -# `fuser(_, virtualspace(mps, i), virtualspace(mpo, i))`, such that the factor that spans them can be -# used as the fuser on the next site: `Fₗ` is oriented as `Vfused ← (Vmps ⊗ Vmpo)`, `Fᵣ` as -# `(Vmps ⊗ Vmpo) ← Vfused`. +""" +Fuse the left virtual legs of the product of MPO-MPS tensors +``` + ┌---A--- 3 + 1 --Fl | => A′[1 2; 3 4] + └---O--- 4 + | + 2 +``` +""" function _fuse_mpo_mps_left(O::MPOTensor, A::MPSTensor, Fₗ) @plansor A′[-1 -2; -3 -4] := Fₗ[-1; 1 3] * A[1 2; -3] * O[3 -2; 2 -4] return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end -# note the index order of the physical leg `-3`: it is the last one, such that the resulting tensor is -# partitioned across the new bond and `_transpose_front` of the right factor is again an MPS tensor +""" +Fuse the right virtual legs of the product of MPO-MPS tensors +``` + 1 --A---┐ + | Fr-- 3 => A′[1 2; 3 4] + 2 --O---┘ + | + 4 +``` +""" function _fuse_mpo_mps_right(O::MPOTensor, A::MPSTensor, Fᵣ) - @plansor A′[-1 -2; -4 -3] := A[-1 1; 2] * O[-2 -3; 1 3] * Fᵣ[2 3; -4] + # the resulting tensor is partitioned across the new bond + # `_transpose_front` of the right factor is again an MPS tensor + @plansor A′[-1 -2; -3 -4] := A[-1 1; 2] * O[-2 -4; 1 3] * Fᵣ[2 3; -3] return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′ end