diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index e239b0989..699dc5ddd 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -498,6 +498,22 @@ @misc{shen2025 keywords = {Condensed Matter - Strongly Correlated Electrons,High Energy Physics - Theory} } +@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 = {cond-mat.str-el} +} + % NOTE: lecture notes; MPSKit citation not independently confirmed during review. @misc{sinha2025, title = {Lectures on Quantum Field Theory on a Quantum Computer}, 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 diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 5b0b2ba32..ba4a3aeea 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -39,6 +39,7 @@ export TDVP, TDVP2, WI, WII, TaylorCluster export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S +export Zipup export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -192,6 +193,7 @@ include("algorithms/statmech/idmrg.jl") include("algorithms/fidelity_susceptibility.jl") include("algorithms/approximate/approximate.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 03b90bb16..1c699397f 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -4,10 +4,12 @@ approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ) approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) + 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 `ψ₀`. 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. @@ -35,12 +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 | `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 new file mode 100644 index 000000000..5ce3ee04d --- /dev/null +++ b/src/algorithms/approximate/zipup.jl @@ -0,0 +1,182 @@ +""" +$(TYPEDEF) + +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 `ϕ`. +Both return the truncation error `ϵ` alongside the approximated state. + +## Fields + +$(TYPEDFIELDS) + +## Constructors + + 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`. +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. + +## References + +- [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010) +- [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) +""" +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::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 + +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(), left_to_right::Bool = true) + if trunc isa TruncationStrategy + 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; left_to_right) + else + throw(ArgumentError("`trunc` should be a truncation strategy or a tuple of two truncation strategies")) + end +end + +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")) + end + + 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 +end + +function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) + O, ϕ = Oϕ + T = TensorOperations.promote_contract(scalartype(O), scalartype(ϕ)) + return approximate!(similar(ϕ, T), Oϕ, alg) +end + +@doc """ + zip_left_right!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ + zip_right_left!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ + +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 and both sweeps. +""" +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: + # 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(ψ))) + + # 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_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 + +function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing) + N = length(ψ) + + Aϕs = map(i -> i == N ? ϕ.AC[N] : ϕ.AL[i], 1:N) + foreach(f -> fill!(f, missing), (ψ.ALs, ψ.ARs, ψ.ACs, ψ.Cs)) + + 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(ψ))) + + # 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 diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index ddda45ff1..c5c033fe4 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -278,10 +278,50 @@ 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 +""" +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 +""" +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ᵣ) + # 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 function Base.:*(mpo::FiniteMPO{<:MPOTensor}, x::AbstractTensorMap) @assert length(mpo) > 1 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) 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