From 2fcd29a42b855ab125161ceedd1adba0b347db71 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 6 Oct 2026 15:45:59 -0400 Subject: [PATCH 1/8] Add nonperturbative cluster expansion for periodic Hamiltonians Construct time-evolution MPOs by matching exact finite-cluster exponentials, including all translations of mixed-space periodic unit cells. Use TensorKit primitives and generated planar contractions, with a single dispatch on the cluster size and inferred internal stages. Complete the SVD factors with complementary directions to preserve environment support for rectangular and rank-deficient residuals. Retain below-cutoff components and use untruncated LQ/QR bases to avoid redundant virtual channels. Document the API and tensor layouts. Validate exact clusters, symmetry sectors, dual and planar spaces, repeated cells, disconnected dimers, time-step edge cases, error scaling, and inference: 1,062 focused tests pass. --- docs/src/man/algorithms.md | 1 + src/MPSKit.jl | 3 +- src/algorithms/timestep/clusterexpansion.jl | 371 ++++++++++++++++++++ test/timeevolution/clusterexpansion.jl | 234 ++++++++++++ 4 files changed, 608 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/timestep/clusterexpansion.jl create mode 100644 test/timeevolution/clusterexpansion.jl diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 7b4a1e1fb..89c0d9ab2 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -123,6 +123,7 @@ Here, there are several algorithms available WI WII TaylorCluster +ClusterExpansion ``` See also [Time evolution accuracy](@ref). diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 4111185cc..ea209f449 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 AlgorithmInfo, convergence_measure -export TDVP, TDVP2, BUG, WI, WII, TaylorCluster +export TDVP, TDVP2, BUG, WI, WII, TaylorCluster, ClusterExpansion export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S @@ -174,6 +174,7 @@ include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/bug.jl") include("algorithms/timestep/taylorcluster.jl") include("algorithms/timestep/wii.jl") +include("algorithms/timestep/clusterexpansion.jl") include("algorithms/timestep/integrators.jl") include("algorithms/timestep/time_evolve.jl") diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl new file mode 100644 index 000000000..4aeaf1db2 --- /dev/null +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -0,0 +1,371 @@ +""" + ClusterExpansion(N; tol = 1.0e-12) + ClusterExpansion(; N = 2, tol = 1.0e-12) + +Construct an infinite time evolution MPO for a nearest-neighbor Hamiltonian by +matching exact exponentials on clusters of up to `N` sites. This is the nonperturbative construction of +[Vanhecke, Vanderstraeten and Verstraete (2021)](https://doi.org/10.1103/PhysRevA.103.L020402). +For nearest-neighbor Hamiltonians the error is `O(dt^N)`. + +The unit cell may contain different physical spaces of the same TensorKit space type. +Symmetry sectors are preserved throughout the construction. Every translated cluster +is matched, including those crossing the unit-cell boundary. + +For full-rank uniform physical spaces of dimension `d`, the virtual dimension grows +as `1 + d^2 + ⋯ + d^(2 floor(N/2))`. Mixed physical spaces and rank-deficient residuals +can require additional complementary virtual channels. Exact cluster exponentials +also grow exponentially with `N`. `tol` is the relative singular-value cutoff used +to identify shared SVD support and compute environment pseudoinverses; components +below the SVD cutoff are retained in the complementary channels. +""" +struct ClusterExpansion <: Algorithm + N::Int + tol::Float64 + function ClusterExpansion(N::Integer; tol::Real = 1.0e-12) + N ≥ 1 || throw(ArgumentError("cluster size must be positive")) + isfinite(tol) && 0 < tol < 1 || throw(ArgumentError("tol must be between zero and one")) + return new(N, tol) + end +end +ClusterExpansion(; N::Integer = 2, kwargs...) = ClusterExpansion(N; kwargs...) + +function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::ClusterExpansion) + # The only dispatch on a runtime cluster size. All subsequent stages and + # contraction shapes are determined by this Val parameter. + return make_cluster_mpo(H, -im * dt, Val(alg.N), alg.tol) +end + +function make_cluster_mpo(H::InfiniteMPOHamiltonian, τ::Number, n::Val, tol::Real) + T = promote_type(scalartype(H), typeof(τ)) + storage = TensorKit.similarstoragetype(storagetype(H), T) + O = PeriodicArray([cluster_initial_tensor(H[i], τ, storage) for i in 1:length(H)]) + # These cases are already exact products; avoid building channels from + # roundoff in differences of independently contracted exponentials. + (iszero(τ) || all(cluster_is_onsite, parent(H))) && return InfiniteMPO(O) + return InfiniteMPO(cluster_build(O, H, τ, n, tol)) +end + +cluster_is_onsite(h) = size(h, 1) == 2 && size(h, 4) == 2 + +function cluster_initial_tensor(h, τ, ::Type{A}) where {A} + P = physicalspace(h) + levels = ⊞(oneunit(P)) + blocktype = tensormaptype(spacetype(P), 2, 2, A) + O = SparseBlockTensorMap{blocktype}(undef, levels ⊗ P ← P ⊗ levels) + onsite = removeunit(removeunit(h[1, 1, 1, end], 4), 1) + O[1, 1, 1, 1] = add_util_leg(exp(τ * onsite)) + return O +end + +cluster_build(O, H, τ, ::Val{1}, tol) = O +function cluster_build(O, H, τ, n::Val{N}, tol) where {N} + previous = cluster_build(O, H, τ, Val(N - 1), tol) + nt = (N - 1) ÷ 2 + level = N ÷ 2 + 1 + # Every translated window must see the completed previous stage. Compute + # all corrections before inserting any of them, including across the cell seam. + if isodd(N) + centers = [cluster_correction(H, previous, s, τ, n, tol) for s in 1:length(H)] + for s in 1:length(H) + previous[s + nt][level, 1, 1, level] = centers[s] + end + return previous + else + pairs = [cluster_pair(H, previous, s, τ, n, tol) for s in 1:length(H)] + expanded = cluster_expand(previous, pairs, n) + for s in 1:length(H) + expanded[s + nt][level - 1, 1, 1, level] = pairs[s][1] + expanded[s + nt + 1][level, 1, 1, level - 1] = pairs[s][2] + end + return expanded + end +end + +function cluster_expand(O, pairs, ::Val{N}) where {N} + nt = (N - 1) ÷ 2 + bonds = PeriodicArray([right_virtualspace(pairs[mod1(i - nt, length(O))][1]) for i in 1:length(O)]) + return PeriodicArray([cluster_expand_tensor(O[i], bonds[i - 1], bonds[i]) for i in 1:length(O)]) +end + +function cluster_expand_tensor(O, left_bond, right_bond) + L, R = left_virtualspace(O), right_virtualspace(O) + left_levels = L ⊞ left_bond + right_levels = R ⊞ right_bond + P = physicalspace(O) + expanded = typeof(O)(undef, left_levels ⊗ P ← P ⊗ right_levels) + for (indices, block) in nonzero_pairs(O) + expanded[indices] = block + end + return expanded +end + +function cluster_residual(H, O, start::Int, τ, n::Val{N}) where {N} + boundary = size(H[start + N - 1], 4) + exact = exp(τ * cluster_contract(H, start, n, boundary)) + return add_util_leg(exact - cluster_contract(O, start, n, 1)) +end + +function cluster_correction(H, O, start, τ, n::Val, tol) + residual = cluster_residual(H, O, start, τ, n) + return cluster_center(O, start, residual, n, tol) +end + +function cluster_pair(H, O, start, τ, n::Val{2}, tol) + residual = cluster_residual(H, O, start, τ, n) + return cluster_split(residual, tol) +end + +function cluster_pair(H, O, start, τ, n::Val{N}, tol) where {N} + residual = cluster_residual(H, O, start, τ, n) + nt = (N - 1) ÷ 2 + left, right = cluster_environments(O, start, n) + left_map, right_map = cluster_environment_maps(left, right, Val(nt)) + center = cluster_center(left_map, right_map, residual, n, tol) + C = permute(center, ((1, 2, 4), (3, 5, 6))) + + # The pseudoinverses restrict the center to the active outer-environment + # subspaces. Complete only these directions, rather than the redundant + # virtual directions introduced by previous stages. + # Completed factors span each physical environment space. Untruncated + # LQ/QR bases retain all these directions, including weak symmetry sectors; + # a global cutoff here would disagree with the blockwise pseudoinverses. + _, Q_left = right_orth(left_map) + Q_right, _ = left_orth(right_map) + F_left = Q_left' ⊗ id(storagetype(C), codomain(C)[2] ⊗ codomain(C)[3]) + F_right = id(storagetype(C), domain(C)[1] ⊗ domain(C)[2]) ⊗ Q_right + A, B = cluster_complete_svd(F_left' * C * F_right, tol) + return permute(F_left * A, ((1, 2), (3, 4))), + permute(B * F_right', ((1, 2), (3, 4))) +end + +function cluster_split(center, tol) + A, B = cluster_complete_svd(permute(center, ((1, 2, 4), (3, 5, 6))), tol) + return permute(A, ((1, 2), (3, 4))), permute(B, ((1, 2), (3, 4))) +end + +""" + cluster_complete_svd(C, tol) + +Factor `C` through its shared SVD support and separate left/right complementary +spaces. For exact rank `r_c`, the virtual multiplicity in sector `c` is +`m_c + n_c - r_c`, the minimum that gives both factors full environment support. +For exact null directions, the factors have the form + + A = [U√S γU⊥ 0], B = [√S Vᴴ; 0; γV⊥ᴴ]. + +At a finite rank cutoff, retain the remaining block `E = U⊥' C V⊥` in the +complementary channels, so the product still reconstructs `C` up to roundoff. +The scale `γ` balances these channels against the retained singular values. +""" +function cluster_complete_svd(C::TensorMap, tol::Real) + cnorm = norm(C) + if iszero(cnorm) + # At zero rank, the two complementary spaces occupy disjoint channels. + L, R = fuse(codomain(C)), fuse(domain(C)) + storage = storagetype(C) + A = catdomain( + isomorphism(storage, codomain(C) ← L), zeros(storage, codomain(C) ← R) + ) + B = catcodomain( + zeros(storage, L ← domain(C)), isomorphism(storage, R ← domain(C)) + ) + return A, B + end + U, S, Vᴴ, _ = svd_trunc(C; trunc = trunctol(; rtol = tol, p = Inf)) + U_perp, V_perp = left_null(U), right_null(Vᴴ) + E = U_perp' * C * V_perp' + γ = sqrt(cnorm) + root = sqrt(S) + A = catdomain(catdomain(U * root, γ * U_perp), U_perp * E / (2γ)) + B = catcodomain(catcodomain(root * Vᴴ, E * V_perp / (2γ)), γ * V_perp) + return A, B +end + +function cluster_contract(O, start::Int, ::Val{1}, boundary::Int) + return removeunit(removeunit(TensorMap(O[start][1, :, :, boundary]), 4), 1) +end + +function cluster_contract(O, start::Int, n::Val{N}, boundary::Int) where {N} + left = removeunit(TensorMap(O[start][1, :, :, :]), 1) + right = removeunit(TensorMap(O[start + N - 1][:, :, :, boundary]), 4) + sites = cluster_bulk(O, start + 1, Val(N - 2)) + return cluster_contract(left, sites, right, n) +end + +cluster_bulk(O, start, ::Val{0}) = () +function cluster_bulk(O, start, ::Val{N}) where {N} + return (TensorMap(O[start]), cluster_bulk(O, start + 1, Val(N - 1))...) +end + +""" + cluster_contract(left, sites, right, ::Val{N}) + +Contract the open-boundary cluster, shown for `N = 4`. `L` and `R` are the +boundary tensors; `O₂` and `O₃` are the site-dependent bulk tensors. Negative labels are open +physical legs and positive labels are contracted virtual bonds. Domain legs +are above the tensors and codomain legs below. + +``` + -5 -6 -7 -8 + │ │ │ │ +┌─┴─┐ ┌─┴─┐ ┌─┴─┐ ┌─┴─┐ +│ L ├─1─┤O₂ ├─2─┤O₃ ├─3─┤ R │ +└─┬─┘ └─┬─┘ └─┬─┘ └─┬─┘ + │ │ │ │ + -1 -2 -3 -4 +``` +""" +@generated function cluster_contract( + left::AbstractTensorMap{<:Any, S, 1, 2}, sites::NTuple{K, AbstractTensorMap{<:Any, S, 2, 2}}, + right::AbstractTensorMap{<:Any, S, 2, 1}, ::Val{N} + ) where {S, N, K} + out = tensorexpr(:cluster, -(1:N), -((N + 1):(2N))) + first_site = tensorexpr(:left, -1, (-N - 1, 1)) + last_site = tensorexpr(:right, (N - 1, -N), -2N) + middle = [tensorexpr(:(sites[$(i - 1)]), (i - 1, -i), (-N - i, i)) for i in 2:(N - 1)] + return macroexpand(@__MODULE__, :(return @plansor $out := *($first_site, $(middle...), $last_site))) +end + +function cluster_environments(O, start, ::Val{N}) where {N} + nt = (N - 1) ÷ 2 + return cluster_left_environment(O, start, Val(nt)), + cluster_right_environment(O, start + N - nt, Val(nt)) +end + +cluster_left_environment(O, start, ::Val{1}) = O[start][1, 1, 1, 2] +function cluster_left_environment(O, start, ::Val{L}) where {L} + left = cluster_left_environment(O, start, Val(L - 1)) + return cluster_join_left(left, O[start + L - 1][L, 1, 1, L + 1]) +end + +cluster_right_environment(O, start, ::Val{1}) = O[start][2, 1, 1, 1] +function cluster_right_environment(O, start, ::Val{L}) where {L} + right = cluster_right_environment(O, start + 1, Val(L - 1)) + return cluster_join_right(O[start][L + 1, 1, 1, L], right) +end + +""" + cluster_join_left(left, site) + +Append a site to the left environment, shown for `K = numout(left) = 2`. +For general `K`, each vertical leg of `left` is a bundle of `K - 1` physical +legs. The outer virtual legs remain open and bond `1` is contracted. Physical +domain legs are above the tensors and codomain legs below. + +``` + -4 -5 + │ │ + ┌──┴──┐ ┌──┴──┐ +-1───┤left ├─1─┤site ├───-6 + └──┬──┘ └──┬──┘ + │ │ + -2 -3 +``` +""" +@generated function cluster_join_left(left::AbstractTensorMap{<:Any, <:Any, K, K}, site::AbstractTensorMap) where {K} + out = tensorexpr(:next_left, -(1:(K + 1)), -((K + 2):(2K + 2))) + env = tensorexpr(:left, -(1:K), (Tuple(-((K + 2):(2K)))..., 1)) + op = tensorexpr(:site, (1, -K - 1), (-2K - 1, -2K - 2)) + return macroexpand(@__MODULE__, :(return @plansor $out := $env * $op)) +end + +""" + cluster_join_right(site, right) + +Prepend a site to the right environment, shown for `K = numout(right) = 2`. +For general `K`, each vertical leg of `right` is a bundle of `K - 1` physical +legs. The outer virtual legs remain open and bond `1` is contracted. Physical +domain legs are above the tensors and codomain legs below. + +``` + -4 -5 + │ │ + ┌──┴──┐ ┌──┴──┐ +-1───┤site ├─1─┤right├───-6 + └──┬──┘ └──┬──┘ + │ │ + -2 -3 +``` +""" +@generated function cluster_join_right(site::AbstractTensorMap, right::AbstractTensorMap{<:Any, <:Any, K, K}) where {K} + out = tensorexpr(:next_right, -(1:(K + 1)), -((K + 2):(2K + 2))) + op = tensorexpr(:site, (-1, -2), (-K - 2, 1)) + env = tensorexpr(:right, (1, Tuple(-(3:(K + 1)))...), -((K + 3):(2K + 2))) + return macroexpand(@__MODULE__, :(return @plansor $out := $op * $env)) +end + +function cluster_environment_maps(left, right, ::Val{L}) where {L} + indices = ntuple(identity, Val(2L + 1)) + return permute(left, (indices, (2L + 2,))), + permute(right, ((1,), indices .+ 1)) +end + +function cluster_center(O, start::Int, residual::TensorMap, n::Val{N}, tol::Real) where {N} + left, right = cluster_environments(O, start, n) + left_map, right_map = cluster_environment_maps(left, right, Val((N - 1) ÷ 2)) + return cluster_center(left_map, right_map, residual, n, tol) +end + +function cluster_center(left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N}, tol::Real) where {N} + nt = (N - 1) ÷ 2 + # Blockwise pseudoinverses preserve symmetry sectors and dual spaces. + physical_indices = ntuple(identity, Val(nt)) + left_inverse = permute( + pinv(left_map; rtol = tol), + (((physical_indices .+ (nt + 2))..., 1), ntuple(identity, Val(nt + 1)) .+ 1) + ) + right_inverse = permute( + pinv(right_map; rtol = tol), + (ntuple(identity, Val(nt + 1)) .+ nt, (2nt + 2, physical_indices...)) + ) + return cluster_apply_inverses(left_inverse, residual, right_inverse, Val(N)) +end + +""" + cluster_apply_inverses(left_inverse, residual, right_inverse, ::Val{N}) + +Apply the environment pseudoinverses `L⁺` and `R⁺` to the residual, shown for +`N = 3`. Each tensor has its domain above and codomain below. Matching positive +labels are contracted; negative labels are the open legs of `center`. The tensors +are drawn separately to keep this orientation explicit. Bonds `1` and `2` are +the residual's trivial boundary legs. + +``` + 1 4 3 -3 5 2 -4 6 + │ │ │ │ │ │ │ │ + ┌──┴───┴──┐ ┌──┴────┴───┴───┴──┐ ┌──┴───┴──┐ + │ L⁺ │ │ residual │ │ R⁺ │ + └──┬───┬──┘ └──┬────┬───┬───┬──┘ └──┬───┬──┘ + │ │ │ │ │ │ │ │ + 3 -1 1 4 -2 6 5 2 +``` + +For larger clusters, each upper/lower environment connection is a bundle of +`nt = (N - 1) ÷ 2` physical legs. The center has one upper/lower physical leg +for odd `N`, or two for even `N`, in addition to its two virtual legs. +""" +@generated function cluster_apply_inverses( + left_inverse::AbstractTensorMap, residual::AbstractTensorMap, + right_inverse::AbstractTensorMap, ::Val{N} + ) where {N} + nt = (N - 1) ÷ 2 + width = iseven(N) ? 2 : 1 + # Positive labels connect the environments to the corresponding physical + # input/output legs of the exact cluster. Negative labels leave only the + # central site(s) and their virtual bonds open. + left_input = Tuple(3:(nt + 2)) + left_output = Tuple((nt + 3):(2nt + 2)) + right_input = Tuple((2nt + 3):(3nt + 2)) + right_output = Tuple((3nt + 3):(4nt + 2)) + center_output = Tuple(-(2:(width + 1))) + center_input = Tuple(-((width + 2):(2width + 1))) + right_bond = -2width - 2 + out = tensorexpr(:center, (-1, center_output...), (center_input..., right_bond)) + left = tensorexpr(:left_inverse, (left_input..., -1), (1, left_output...)) + rhs = tensorexpr( + :residual, (1, left_output..., center_output..., right_output...), + (left_input..., center_input..., right_input..., 2) + ) + right = tensorexpr(:right_inverse, (right_input..., 2), (right_bond, right_output...)) + return macroexpand(@__MODULE__, :(return @plansor $out := $left * $rhs * $right)) +end diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl new file mode 100644 index 000000000..ba069056b --- /dev/null +++ b/test/timeevolution/clusterexpansion.jl @@ -0,0 +1,234 @@ +using Test +using MPSKit, TensorKit, LinearAlgebra, Random +using TensorKit: ℙ + +# Construct finite windows independently of the cluster-contraction kernels, +# including lengths that are not multiples of the infinite unit cell. +function cluster_dense_window(U, start, L) + sites = [copy(U[start + i - 1]) for i in 1:L] + sites[1] = sites[1][1, :, :, :] + sites[end] = sites[end][:, :, :, 1] + return convert(TensorMap, MPSKit.DenseMPO(FiniteMPO(sites))) +end + +function cluster_exact_window(H, τ, start, L) + if L == 1 + onsite = H[start][1, 1, 1, end] + return exp(τ * MPSKit.removeunit(MPSKit.removeunit(onsite, 4), 1)) + end + sites = [copy(H[start + i - 1]) for i in 1:L] + sites[1] = sites[1][1, :, :, :] + sites[end] = sites[end][:, :, :, end] + return exp(τ * convert(TensorMap, FiniteMPOHamiltonian(sites))) +end + +@testset "Nonperturbative cluster expansion" begin + P = ℂ^2 + X = TensorMap(ComplexF64[0 1; 1 0], P ← P) + Z = TensorMap(ComplexF64[1 0; 0 -1], P ← P) + h = Z ⊗ Z + 0.3 * (X ⊗ X) + H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.7X, (1, 2) => h) + dense(U, L) = convert(TensorMap, MPSKit.DenseMPO(open_boundary_conditions(U, L))) + exact(H, dt, L) = L == 1 ? exp(-im * dt * 0.7X) : + exp(-im * dt * convert(TensorMap, open_boundary_conditions(H, L))) + + @testset "Exact clusters, N=$N" for N in 1:5 + U = make_time_mpo(H, 0.1, ClusterExpansion(N)) + @test U isa InfiniteMPO + @test dim(left_virtualspace(U, 1)) == sum(4^l for l in 0:(N ÷ 2)) + for L in 1:N + @test dense(U, L) ≈ exact(H, 0.1, L) atol = 1.0e-10 + end + end + + @testset "Time steps and rank-deficient environments" begin + onsite = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.7X) + commuting = InfiniteMPOHamiltonian(PeriodicArray([P]), (1, 2) => Z ⊗ Z) + for dt in (0.0, -0.1, -0.1im) + U = make_time_mpo(H, dt, ClusterExpansion(3)) + @test dense(U, 3) ≈ exact(H, dt, 3) atol = 1.0e-10 + end + for Hsimple in (onsite, commuting) + U = make_time_mpo(Hsimple, 0.1, ClusterExpansion(3)) + @test dense(U, 4) ≈ exp(-0.1im * convert(TensorMap, open_boundary_conditions(Hsimple, 4))) atol = 1.0e-10 + end + end + + @testset "Complex Hamiltonian and reflection symmetry" begin + Y = TensorMap(ComplexF64[0 -im; im 0], P ← P) + complex_h = h + 0.2 * (X ⊗ Y + Y ⊗ X) + complex_H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.4Y, (1, 2) => complex_h) + for N in (3, 4), dt in (0.1, -0.1im) + U = make_time_mpo(complex_H, dt, ClusterExpansion(N)) + @test dense(U, N) ≈ exp(-im * dt * convert(TensorMap, open_boundary_conditions(complex_H, N))) atol = 1.0e-10 + larger = dense(U, N + 1) + reflected = permute(larger, (Tuple((N + 1):-1:1), Tuple((2N + 2):-1:(N + 2)))) + @test larger ≈ reflected atol = 1.0e-10 + end + end + + @testset "Leading error on larger chains" begin + for N in 2:4 + errors = Float64[] + for dt in (0.04, 0.02) + U = make_time_mpo(H, dt, ClusterExpansion(N)) + push!(errors, norm(dense(U, N + 1) - exact(H, dt, N + 1))) + end + @test log2(errors[1] / errors[2]) ≥ N - 0.15 + end + end + + @testset "Symmetry sectors and dual physical spaces" begin + rng = MersenneTwister(1234) + spaces = ( + ℝ^2, ℙ^2, + Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1)), + Rep[SU₂](1 // 2 => 1), dual(Rep[SU₂](1 // 2 => 1)), + Vect[FermionParity](0 => 1, 1 => 1), + ) + for P in spaces + raw1 = randn!(rng, zeros(ComplexF64, P ← P)) + onsite = (raw1 + raw1') / 2 + raw2 = randn!(rng, zeros(ComplexF64, P ⊗ P ← P ⊗ P)) + interaction = (raw2 + raw2') / 2 + symmetric_H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => onsite, (1, 2) => interaction) + U = make_time_mpo(symmetric_H, 0.1, ClusterExpansion(1)) + @test dense(U, 1) ≈ exp(-0.1im * onsite) atol = 1.0e-10 + for N in 2:5 + U = @inferred MPSKit.make_cluster_mpo(symmetric_H, -0.1im, Val(N), 1.0e-12) + @test physicalspace(U, 1) == P + @test left_virtualspace(U, 1)[2] == fuse(P ⊗ dual(P)) + exact_cluster = exp(-0.1im * convert(TensorMap, open_boundary_conditions(symmetric_H, N))) + @test dense(U, N) ≈ exact_cluster atol = 1.0e-10 + end + residual = @inferred MPSKit.cluster_residual(symmetric_H, U, 1, -0.1im, Val(5)) + center = @inferred MPSKit.cluster_center(U, 1, residual, Val(5), 1.0e-12) + @test center isa TensorMap + end + end + + @testset "Inference after the cluster-size dispatch" begin + for n in (Val(4), Val(5)) + U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, n, 1.0e-12) + residual = @inferred MPSKit.cluster_residual(H, U, 1, -0.1im, n) + center = @inferred MPSKit.cluster_center(U, 1, residual, n, 1.0e-12) + @test center isa TensorMap + end + end + + @test_throws ArgumentError ClusterExpansion(0) + @test_throws ArgumentError ClusterExpansion(2; tol = 0) + @test_throws ArgumentError ClusterExpansion(2; tol = Inf) + @test ClusterExpansion(; N = 3).N == 3 + H2 = repeat(H, 2) + @testset "Repeated unit cell" for N in 2:5 + U = make_time_mpo(H, 0.1, ClusterExpansion(N)) + U2 = make_time_mpo(H2, 0.1, ClusterExpansion(N)) + @test length(U2) == 2 + @test cluster_dense_window(U2, 2, N + 1) ≈ dense(U, N + 1) atol = 1.0e-10 + end +end + +@testset "Periodic cluster expansion with mixed physical spaces" begin + rng = MersenneTwister(25) + lattices = ( + [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], + [ℝ^2, ℝ^3], [ℙ^2, ℙ^3], + [Rep[U₁](0 => 1, 1 => 1), Rep[U₁](0 => 1, 1 => 1, 2 => 1)], + [dual(Rep[U₁](0 => 1, 1 => 1)), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1))], + [Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1))], + [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1)], + [dual(Rep[SU₂](1 // 2 => 1)), dual(Rep[SU₂](1 => 1))], + [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1)], + ) + for lattice in lattices + Ps = PeriodicArray(lattice) + terms = Pair[] + for s in eachindex(lattice) + a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) + b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) + push!(terms, s => (a + a') / 2, (s, s + 1) => (b + b') / 2) + end + H = InfiniteMPOHamiltonian(Ps, terms...) + for N in 1:5 + U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), 1.0e-12) + @test length(U) == length(H) + for s in eachindex(lattice) + @test physicalspace(U, s) == Ps[s] + @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) + for L in 1:N + @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(H, -0.1im, s, L) atol = 1.0e-10 + end + end + end + end +end + +@testset "Completed SVD support" begin + rng = MersenneTwister(22) + spaces = ( + (ℂ^4, ℂ^9), + (Rep[U₁](0 => 4, 1 => 1), Rep[U₁](0 => 2, -1 => 2, 1 => 3)), + (Rep[SU₂](0 => 3, 1 // 2 => 2, 1 => 1), Rep[SU₂](0 => 1, 1 // 2 => 3, 2 => 1)), + (Vect[FermionParity](0 => 4, 1 => 1), Vect[FermionParity](0 => 2, 1 => 3)), + ) + for (L, R) in spaces, dualize in (false, true) + L, R = dualize ? (dual(L), dual(R)) : (L, R) + raw = randn!(rng, zeros(ComplexF64, L ← R)) + for scale in (1.0, 1.0e-28, 0.0), tol in (1.0e-12, 0.5) + C = scale * raw + A, B = @inferred MPSKit.cluster_complete_svd(C, tol) + @test norm(A * B - C) ≤ 2.0e-12 * max(norm(C), 1) + @test domain(A) == codomain(B) + @test A * pinv(A; rtol = 1.0e-12) ≈ id(storagetype(C), L) atol = 1.0e-10 + @test pinv(B; rtol = 1.0e-12) * B ≈ id(storagetype(C), R) atol = 1.0e-10 + if scale == 1 && tol == 1.0e-12 + V = only(domain(A)) + for c in union(sectors(L), sectors(R)) + @test dim(V, c) == max(dim(L, c), dim(R, c)) + end + end + end + end + # The SVD rank cutoff must not discard the tiny second singular value. + C = TensorMap(ComplexF64[1 0 0; 0 1.0e-14 0; 0 0 0], ℂ^3 ← ℂ^3) + A, B = MPSKit.cluster_complete_svd(C, 1.0e-12) + @test dim(only(domain(A))) == 5 + @test norm(A * B - C) < 1.0e-16 +end + +@testset "Mixed-cell edge cases and leading error" begin + rng = MersenneTwister(24) + for lattice in ([ℂ^2, ℂ^3], [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1)]) + Ps = PeriodicArray(lattice) + onsite_terms = Pair[] + full_terms = Pair[] + for s in eachindex(lattice) + a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) + b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) + push!(onsite_terms, s => (a + a') / 2) + push!(full_terms, s => (a + a') / 2, (s, s + 1) => (b + b') / 2) + end + onsite = InfiniteMPOHamiltonian(Ps, onsite_terms...) + H = InfiniteMPOHamiltonian(Ps, full_terms...) + # Remove the interaction across the cell seam to form independent dimers. + dimers = InfiniteMPOHamiltonian(Ps, filter(p -> p.first != (2, 3), full_terms)...) + for (model, dt) in ((onsite, 0.1), (dimers, 0.1), (H, 0.0), (H, -0.1im), (H, -0.1)) + U = make_time_mpo(model, dt, ClusterExpansion(5)) + longest = model === onsite || model === dimers || iszero(dt) ? 6 : 5 + for s in eachindex(lattice), L in 1:longest + @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(model, -im * dt, s, L) atol = 1.0e-10 + end + end + for N in 3:4, s in eachindex(lattice) + # Keep the SU(2) N=4 errors above the floating-point floor. + steps = spacetype(Ps[1]) == ComplexSpace || N == 3 ? (0.04, 0.02) : (0.16, 0.08) + errors = Float64[] + for dt in steps + U = make_time_mpo(H, dt, ClusterExpansion(N)) + push!(errors, norm(cluster_dense_window(U, s, N + 1) - cluster_exact_window(H, -im * dt, s, N + 1))) + end + @test log2(errors[1] / errors[2]) ≥ N - 0.15 + end + end +end From 4510336f40d8d5f09af4e6b4964eb594dfdaaa07 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 6 Oct 2026 18:05:03 -0400 Subject: [PATCH 2/8] Solve cluster environment equations without explicit pseudoinverses --- src/algorithms/timestep/clusterexpansion.jl | 127 ++++++++++---------- test/timeevolution/clusterexpansion.jl | 11 +- 2 files changed, 70 insertions(+), 68 deletions(-) diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index 4aeaf1db2..568bdb91d 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -15,8 +15,9 @@ For full-rank uniform physical spaces of dimension `d`, the virtual dimension gr as `1 + d^2 + ⋯ + d^(2 floor(N/2))`. Mixed physical spaces and rank-deficient residuals can require additional complementary virtual channels. Exact cluster exponentials also grow exponentially with `N`. `tol` is the relative singular-value cutoff used -to identify shared SVD support and compute environment pseudoinverses; components -below the SVD cutoff are retained in the complementary channels. +to identify shared SVD support; components below the cutoff are retained in the +complementary channels. Environment equations are solved directly without a +second singular-value cutoff. """ struct ClusterExpansion <: Algorithm N::Int @@ -65,7 +66,7 @@ function cluster_build(O, H, τ, n::Val{N}, tol) where {N} # Every translated window must see the completed previous stage. Compute # all corrections before inserting any of them, including across the cell seam. if isodd(N) - centers = [cluster_correction(H, previous, s, τ, n, tol) for s in 1:length(H)] + centers = [cluster_correction(H, previous, s, τ, n) for s in 1:length(H)] for s in 1:length(H) previous[s + nt][level, 1, 1, level] = centers[s] end @@ -105,9 +106,9 @@ function cluster_residual(H, O, start::Int, τ, n::Val{N}) where {N} return add_util_leg(exact - cluster_contract(O, start, n, 1)) end -function cluster_correction(H, O, start, τ, n::Val, tol) +function cluster_correction(H, O, start, τ, n::Val) residual = cluster_residual(H, O, start, τ, n) - return cluster_center(O, start, residual, n, tol) + return cluster_center(O, start, residual, n) end function cluster_pair(H, O, start, τ, n::Val{2}, tol) @@ -120,15 +121,13 @@ function cluster_pair(H, O, start, τ, n::Val{N}, tol) where {N} nt = (N - 1) ÷ 2 left, right = cluster_environments(O, start, n) left_map, right_map = cluster_environment_maps(left, right, Val(nt)) - center = cluster_center(left_map, right_map, residual, n, tol) + center = cluster_center(left_map, right_map, residual, n) C = permute(center, ((1, 2, 4), (3, 5, 6))) - # The pseudoinverses restrict the center to the active outer-environment + # The minimum-norm solves restrict the center to the active environment # subspaces. Complete only these directions, rather than the redundant - # virtual directions introduced by previous stages. - # Completed factors span each physical environment space. Untruncated - # LQ/QR bases retain all these directions, including weak symmetry sectors; - # a global cutoff here would disagree with the blockwise pseudoinverses. + # virtual directions introduced by previous stages. Untruncated LQ/QR bases + # retain the entire physical environment support, including weak sectors. _, Q_left = right_orth(left_map) Q_right, _ = left_orth(right_map) F_left = Q_left' ⊗ id(storagetype(C), codomain(C)[2] ⊗ codomain(C)[3]) @@ -300,72 +299,68 @@ function cluster_environment_maps(left, right, ::Val{L}) where {L} permute(right, ((1,), indices .+ 1)) end -function cluster_center(O, start::Int, residual::TensorMap, n::Val{N}, tol::Real) where {N} +function cluster_center(O, start::Int, residual::TensorMap, n::Val{N}) where {N} left, right = cluster_environments(O, start, n) left_map, right_map = cluster_environment_maps(left, right, Val((N - 1) ÷ 2)) - return cluster_center(left_map, right_map, residual, n, tol) -end - -function cluster_center(left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N}, tol::Real) where {N} - nt = (N - 1) ÷ 2 - # Blockwise pseudoinverses preserve symmetry sectors and dual spaces. - physical_indices = ntuple(identity, Val(nt)) - left_inverse = permute( - pinv(left_map; rtol = tol), - (((physical_indices .+ (nt + 2))..., 1), ntuple(identity, Val(nt + 1)) .+ 1) - ) - right_inverse = permute( - pinv(right_map; rtol = tol), - (ntuple(identity, Val(nt + 1)) .+ nt, (2nt + 2, physical_indices...)) - ) - return cluster_apply_inverses(left_inverse, residual, right_inverse, Val(N)) + return cluster_center(left_map, right_map, residual, n) end """ - cluster_apply_inverses(left_inverse, residual, right_inverse, ::Val{N}) + cluster_center(left_map, right_map, residual, ::Val{N}) + +Solve the environment equations from the left and right without forming inverse +maps. Completed environments have full row/column support in each symmetry +sector; rectangular solves choose the minimum-norm center. + +The grouped map layouts for `N = 3` are shown below. Domain bundles are above +and codomain bundles below. `pᵢ′` is the dual physical space, and `u_L`/`u_R` +are the trivial boundary spaces. -Apply the environment pseudoinverses `L⁺` and `R⁺` to the residual, shown for -`N = 3`. Each tensor has its domain above and codomain below. Matching positive -labels are contracted; negative labels are the open legs of `center`. The tensors -are drawn separately to keep this orientation explicit. Bonds `1` and `2` are -the residual's trivial boundary legs. +First group the left physical legs of the residual into `B_L` and solve `L \\ B_L`: ``` - 1 4 3 -3 5 2 -4 6 - │ │ │ │ │ │ │ │ - ┌──┴───┴──┐ ┌──┴────┴───┴───┴──┐ ┌──┴───┴──┐ - │ L⁺ │ │ residual │ │ R⁺ │ - └──┬───┬──┘ └──┬────┬───┬───┬──┘ └──┬───┬──┘ - │ │ │ │ │ │ │ │ - 3 -1 1 4 -2 6 5 2 + V_L (p₂′, p₃′, p₂, p₃, u_R) + │ │ + ┌─────┴─────┐ ┌─────┴─────┐ + │ L │ │ B_L │ + └─────┬─────┘ └─────┬─────┘ + │ │ + (u_L, p₁, p₁′) (u_L, p₁, p₁′) ``` -For larger clusters, each upper/lower environment connection is a bundle of -`nt = (N - 1) ÷ 2` physical legs. The center has one upper/lower physical leg -for odd `N`, or two for even `N`, in addition to its two virtual legs. +Regroup that solution into `B_R` and solve `B_R / R`: + +``` + (p₃′, p₃, u_R) (p₃′, p₃, u_R) + │ │ + ┌─────┴─────┐ ┌─────┴─────┐ + │ R │ │ B_R │ + └─────┬─────┘ └─────┬─────┘ + │ │ + V_R (V_L, p₂, p₂′) +``` + +Finally move `p₂′` back into the domain. The center has one physical leg per +side for odd `N`, or two for even `N`, and its two virtual legs. For larger +clusters, each outer physical space above is a bundle of `(N - 1) ÷ 2` legs. """ -@generated function cluster_apply_inverses( - left_inverse::AbstractTensorMap, residual::AbstractTensorMap, - right_inverse::AbstractTensorMap, ::Val{N} +@generated function cluster_center( + left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N} ) where {N} nt = (N - 1) ÷ 2 - width = iseven(N) ? 2 : 1 - # Positive labels connect the environments to the corresponding physical - # input/output legs of the exact cluster. Negative labels leave only the - # central site(s) and their virtual bonds open. - left_input = Tuple(3:(nt + 2)) - left_output = Tuple((nt + 3):(2nt + 2)) - right_input = Tuple((2nt + 3):(3nt + 2)) - right_output = Tuple((3nt + 3):(4nt + 2)) - center_output = Tuple(-(2:(width + 1))) - center_input = Tuple(-((width + 2):(2width + 1))) - right_bond = -2width - 2 - out = tensorexpr(:center, (-1, center_output...), (center_input..., right_bond)) - left = tensorexpr(:left_inverse, (left_input..., -1), (1, left_output...)) - rhs = tensorexpr( - :residual, (1, left_output..., center_output..., right_output...), - (left_input..., center_input..., right_input..., 2) - ) - right = tensorexpr(:right_inverse, (right_input..., 2), (right_bond, right_output...)) - return macroexpand(@__MODULE__, :(return @plansor $out := $left * $rhs * $right)) + width = N - 2nt + # Group all left environment legs into the codomain for the left solve. + left_indices = (1, Tuple(2:(nt + 1))..., Tuple((N + 2):(N + nt + 1))...) + remaining = (Tuple((nt + 2):(N + 1))..., Tuple((N + nt + 2):(2N + 2))...) + # After solving, group the right environment legs into the domain. + middle_indices = (1, Tuple(2:(width + 1))..., Tuple((nt + width + 2):(nt + 2width + 1))...) + right_indices = (Tuple((width + 2):(nt + width + 1))..., Tuple((nt + 2width + 2):(2nt + 2width + 2))...) + center_indices = (Tuple(1:(width + 1)), Tuple((width + 2):(2width + 2))) + return quote + left_rhs = permute(residual, $((left_indices, remaining))) + left_solved = left_map \ left_rhs + right_rhs = permute(left_solved, $((middle_indices, right_indices))) + solved = right_rhs / right_map + return permute(solved, $center_indices) + end end diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl index ba069056b..9cb915d51 100644 --- a/test/timeevolution/clusterexpansion.jl +++ b/test/timeevolution/clusterexpansion.jl @@ -102,7 +102,7 @@ end @test dense(U, N) ≈ exact_cluster atol = 1.0e-10 end residual = @inferred MPSKit.cluster_residual(symmetric_H, U, 1, -0.1im, Val(5)) - center = @inferred MPSKit.cluster_center(U, 1, residual, Val(5), 1.0e-12) + center = @inferred MPSKit.cluster_center(U, 1, residual, Val(5)) @test center isa TensorMap end end @@ -111,7 +111,7 @@ end for n in (Val(4), Val(5)) U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, n, 1.0e-12) residual = @inferred MPSKit.cluster_residual(H, U, 1, -0.1im, n) - center = @inferred MPSKit.cluster_center(U, 1, residual, n, 1.0e-12) + center = @inferred MPSKit.cluster_center(U, 1, residual, n) @test center isa TensorMap end end @@ -220,6 +220,13 @@ end @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(model, -im * dt, s, L) atol = 1.0e-10 end end + # A large cutoff changes the SVD channels, not the matched exponentials. + for N in (3, 5), tol in (0.1, 0.5) + U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), tol) + for s in eachindex(lattice), L in 1:N + @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(H, -0.1im, s, L) atol = 1.0e-10 + end + end for N in 3:4, s in eachindex(lattice) # Keep the SU(2) N=4 errors above the floating-point floor. steps = spacetype(Ps[1]) == ComplexSpace || N == 3 ? (0.04, 0.02) : (0.16, 0.08) From cdafc1226c52627df93347e09a4bbc6319969aa6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 6 Oct 2026 18:34:08 -0400 Subject: [PATCH 3/8] Support finite chains in nonperturbative cluster expansion Share the cluster builder and tensor kernels between finite and infinite Hamiltonians. Limit finite stages to windows contained in the chain, add empty virtual levels near the endpoints, and project and structurally prune the final MPO. Cap the cluster size at the finite chain length before the single dynamic size dispatch. Handle one-site finite Hamiltonian conversion by stripping both utility legs instead of contracting the site with itself. Update the algorithm documentation and cover mixed spaces, symmetries, boundary windows, inference, large cutoffs, time-step edge cases, finite error scaling, and agreement with projected infinite MPOs. All 2,030 cluster-expansion tests pass. --- src/algorithms/timestep/clusterexpansion.jl | 64 ++++++++++++---- src/operators/mpohamiltonian.jl | 2 + test/hamiltonian/finite.jl | 4 + test/timeevolution/clusterexpansion.jl | 84 ++++++++++++++++++++- 4 files changed, 134 insertions(+), 20 deletions(-) diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index 568bdb91d..3d2ddc06c 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -2,14 +2,17 @@ ClusterExpansion(N; tol = 1.0e-12) ClusterExpansion(; N = 2, tol = 1.0e-12) -Construct an infinite time evolution MPO for a nearest-neighbor Hamiltonian by +Construct a finite or infinite time evolution MPO for a nearest-neighbor Hamiltonian by matching exact exponentials on clusters of up to `N` sites. This is the nonperturbative construction of [Vanhecke, Vanderstraeten and Verstraete (2021)](https://doi.org/10.1103/PhysRevA.103.L020402). For nearest-neighbor Hamiltonians the error is `O(dt^N)`. -The unit cell may contain different physical spaces of the same TensorKit space type. -Symmetry sectors are preserved throughout the construction. Every translated cluster -is matched, including those crossing the unit-cell boundary. +Sites may have different physical spaces of the same TensorKit space type. +Symmetry sectors are preserved throughout the construction. Infinite Hamiltonians +match every translated cluster, including those crossing the unit-cell boundary. +Finite Hamiltonians match every cluster contained in the chain and return an +open-boundary MPO. For a chain of length `L`, the cluster size is capped at `L`; +when `N ≥ L`, the entire finite-chain exponential is matched up to roundoff. For full-rank uniform physical spaces of dimension `d`, the virtual dimension grows as `1 + d^2 + ⋯ + d^(2 floor(N/2))`. Mixed physical spaces and rank-deficient residuals @@ -30,23 +33,38 @@ struct ClusterExpansion <: Algorithm end ClusterExpansion(; N::Integer = 2, kwargs...) = ClusterExpansion(N; kwargs...) -function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::ClusterExpansion) +function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) # The only dispatch on a runtime cluster size. All subsequent stages and # contraction shapes are determined by this Val parameter. - return make_cluster_mpo(H, -im * dt, Val(alg.N), alg.tol) + isempty(H) && throw(ArgumentError("Hamiltonian must contain at least one site")) + N = isfinite(H) ? min(alg.N, length(H)) : alg.N + return make_cluster_mpo(H, -im * dt, Val(N), alg.tol) end -function make_cluster_mpo(H::InfiniteMPOHamiltonian, τ::Number, n::Val, tol::Real) +function make_cluster_mpo(H::MPOHamiltonian, τ::Number, n::Val, tol::Real) T = promote_type(scalartype(H), typeof(τ)) storage = TensorKit.similarstoragetype(storagetype(H), T) - O = PeriodicArray([cluster_initial_tensor(H[i], τ, storage) for i in 1:length(H)]) + O = map(parent(H)) do h + return cluster_initial_tensor(h, τ, storage) + end # These cases are already exact products; avoid building channels from # roundoff in differences of independently contracted exponentials. - (iszero(τ) || all(cluster_is_onsite, parent(H))) && return InfiniteMPO(O) - return InfiniteMPO(cluster_build(O, H, τ, n, tol)) + (iszero(τ) || all(cluster_is_onsite, parent(H))) && return cluster_finalize(O) + return cluster_finalize(cluster_build(O, H, τ, n, tol)) +end + +cluster_is_onsite(h) = size(h, 1) ≤ 2 && size(h, 4) ≤ 2 + +cluster_finalize(O::PeriodicVector) = InfiniteMPO(O) +function cluster_finalize(O::Vector) + O[1] = O[1][1, :, :, :] + O[end] = O[end][:, :, :, 1] + # Remove structurally unused levels only after all stages have finished. + return remove_orphans!(FiniteMPO(O); tol = 0) end -cluster_is_onsite(h) = size(h, 1) == 2 && size(h, 4) == 2 +cluster_starts(H::InfiniteMPOHamiltonian, ::Val) = 1:length(H) +cluster_starts(H::FiniteMPOHamiltonian, ::Val{N}) where {N} = 1:(length(H) - N + 1) function cluster_initial_tensor(h, τ, ::Type{A}) where {A} P = physicalspace(h) @@ -61,20 +79,21 @@ end cluster_build(O, H, τ, ::Val{1}, tol) = O function cluster_build(O, H, τ, n::Val{N}, tol) where {N} previous = cluster_build(O, H, τ, Val(N - 1), tol) + starts = cluster_starts(H, n) nt = (N - 1) ÷ 2 level = N ÷ 2 + 1 # Every translated window must see the completed previous stage. Compute # all corrections before inserting any of them, including across the cell seam. if isodd(N) - centers = [cluster_correction(H, previous, s, τ, n) for s in 1:length(H)] - for s in 1:length(H) + centers = [cluster_correction(H, previous, s, τ, n) for s in starts] + for s in starts previous[s + nt][level, 1, 1, level] = centers[s] end return previous else - pairs = [cluster_pair(H, previous, s, τ, n, tol) for s in 1:length(H)] + pairs = [cluster_pair(H, previous, s, τ, n, tol) for s in starts] expanded = cluster_expand(previous, pairs, n) - for s in 1:length(H) + for s in starts expanded[s + nt][level - 1, 1, 1, level] = pairs[s][1] expanded[s + nt + 1][level, 1, 1, level - 1] = pairs[s][2] end @@ -82,12 +101,25 @@ function cluster_build(O, H, τ, n::Val{N}, tol) where {N} end end -function cluster_expand(O, pairs, ::Val{N}) where {N} +function cluster_expand(O::PeriodicVector, pairs, ::Val{N}) where {N} nt = (N - 1) ÷ 2 bonds = PeriodicArray([right_virtualspace(pairs[mod1(i - nt, length(O))][1]) for i in 1:length(O)]) return PeriodicArray([cluster_expand_tensor(O[i], bonds[i - 1], bonds[i]) for i in 1:length(O)]) end +function cluster_expand(O::Vector, pairs, ::Val{N}) where {N} + nt = (N - 1) ÷ 2 + empty_bond = zero(left_virtualspace(O[1])[1]) + # Only the middle bonds can support an N-site cluster. Empty levels keep + # the same level indices everywhere until the final boundary projection. + bonds = vcat( + fill(empty_bond, nt + 1), + [right_virtualspace(pair[1]) for pair in pairs], + fill(empty_bond, nt + 1) + ) + return [cluster_expand_tensor(O[i], bonds[i], bonds[i + 1]) for i in 1:length(O)] +end + function cluster_expand_tensor(O, left_bond, right_bond) L, R = left_virtualspace(O), right_virtualspace(O) left_levels = L ⊞ left_bond diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 6a8c6e055..4b724a4eb 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -884,6 +884,8 @@ function isemptylevel(H::InfiniteMPOHamiltonian, i::Int) end function Base.convert(::Type{TensorMap}, H::FiniteMPOHamiltonian) + # A single site is both boundaries; strip its utility legs only once. + length(H) == 1 && return removeunit(removeunit(TensorMap(only(H)), 4), 1) L = removeunit(H[1], 1) R = removeunit(H[end], 4) M = Tuple(H[2:(end - 1)]) diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index e3c266249..14f7c9794 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -41,6 +41,10 @@ using TensorKit: ℙ end end + # A one-site Hamiltonian must not be contracted with itself. + Hsingle = FiniteMPOHamiltonian([V], 1 => O₁) + @test convert(TensorMap, Hsingle) ≈ O₁ + # check if constructor works by converting back to tensormap H1_tm = convert(TensorMap, H1) operators = vcat(fill(E, L - 1), O₁) diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl index 9cb915d51..664d3425e 100644 --- a/test/timeevolution/clusterexpansion.jl +++ b/test/timeevolution/clusterexpansion.jl @@ -12,10 +12,6 @@ function cluster_dense_window(U, start, L) end function cluster_exact_window(H, τ, start, L) - if L == 1 - onsite = H[start][1, 1, 1, end] - return exp(τ * MPSKit.removeunit(MPSKit.removeunit(onsite, 4), 1)) - end sites = [copy(H[start + i - 1]) for i in 1:L] sites[1] = sites[1][1, :, :, :] sites[end] = sites[end][:, :, :, end] @@ -239,3 +235,83 @@ end end end end + + +# Independent random couplings and onsite terms exercise nonuniform finite chains. +function cluster_finite_model(rng, Ps; onsite_only = false, dimers = false) + terms = Pair[] + for s in eachindex(Ps) + a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) + push!(terms, s => (a + a') / 2) + if s < length(Ps) && !onsite_only && (!dimers || isodd(s)) + b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) + push!(terms, (s, s + 1) => (b + b') / 2) + end + end + return FiniteMPOHamiltonian(Ps, terms...) +end + +@testset "Finite cluster expansion" begin + rng = MersenneTwister(723) + lattices = ( + [ℂ^2], [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], [ℂ^2, ℂ^3, ℂ^2, ℂ^3], + [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2], [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2, ℂ^3], + [ℝ^2, ℝ^3, ℝ^2, ℝ^3, ℝ^2], [ℙ^2, ℙ^3, ℙ^2, ℙ^3, ℙ^2], + [Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1)), Rep[U₁](0 => 1, 1 => 1)], + [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1), Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1), Rep[SU₂](1 // 2 => 1)], + [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1), Vect[FermionParity](0 => 1, 1 => 1)], + ) + for Ps in lattices + H = cluster_finite_model(rng, Ps) + L = length(H) + for N in 1:min(6, L) + U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), 1.0e-12) + @test U isa FiniteMPO + @test length(U) == L + @test dim(left_virtualspace(U, 1)) == dim(right_virtualspace(U, L)) == 1 + for s in 1:L + @test physicalspace(U, s) == Ps[s] + s < L && @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) + for width in 1:min(N, L - s + 1) + @test cluster_dense_window(U, s, width) ≈ cluster_exact_window(H, -0.1im, s, width) atol = 1.0e-10 + end + end + end + # The public entry point caps N at L, including one-site chains. + for N in (L, L + 3), tol in (1.0e-12, 0.5) + U = make_time_mpo(H, 0.1, ClusterExpansion(N; tol)) + @test convert(TensorMap, MPSKit.DenseMPO(U)) ≈ exp(-0.1im * convert(TensorMap, H)) atol = 1.0e-10 + end + end +end + +@testset "Finite edge cases and leading error" begin + rng = MersenneTwister(984) + Ps = [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2] + H = cluster_finite_model(rng, Ps) + onsite = cluster_finite_model(rng, Ps; onsite_only = true) + dimers = cluster_finite_model(rng, Ps; dimers = true) + for (model, dt) in ((H, 0.0), (H, -0.1), (H, -0.1im), (onsite, 0.1), (dimers, 0.1)) + U = make_time_mpo(model, dt, ClusterExpansion(5)) + @test convert(TensorMap, MPSKit.DenseMPO(U)) ≈ exp(-im * dt * convert(TensorMap, model)) atol = 1.0e-10 + end + @test dim(left_virtualspace(make_time_mpo(onsite, 0.1, ClusterExpansion(5)), 3)) == 1 + @test_throws ArgumentError make_time_mpo(FiniteMPOHamiltonian(similar(parent(H), 0)), 0.1, ClusterExpansion(2)) + for N in (3, 4) + errors = Float64[] + for dt in (0.1, 0.05) + U = make_time_mpo(H, dt, ClusterExpansion(N)) + push!(errors, norm(convert(TensorMap, MPSKit.DenseMPO(U)) - exp(-im * dt * convert(TensorMap, H)))) + end + @test log2(errors[1] / errors[2]) ≥ N - 0.15 + end + P = ℂ^2 + a = randn!(rng, zeros(ComplexF64, P ← P)) + b = randn!(rng, zeros(ComplexF64, P ⊗ P ← P ⊗ P)) + infinite = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => (a + a') / 2, (1, 2) => (b + b') / 2) + for L in (4, 6), N in 1:5 + finite = make_time_mpo(open_boundary_conditions(infinite, L), 0.1, ClusterExpansion(N)) + projected = open_boundary_conditions(make_time_mpo(infinite, 0.1, ClusterExpansion(N)), L) + @test convert(TensorMap, MPSKit.DenseMPO(finite)) ≈ convert(TensorMap, MPSKit.DenseMPO(projected)) atol = 1.0e-10 + end +end From 3849ff327b63d37d4b87ece1147e168620033015 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 09:39:01 -0400 Subject: [PATCH 4/8] Clarify cluster expansion helpers and algorithm documentation --- docs/src/assets/mpskit.bib | 14 ++ src/algorithms/timestep/clusterexpansion.jl | 221 ++++++++++---------- test/timeevolution/clusterexpansion.jl | 48 ++--- 3 files changed, 151 insertions(+), 132 deletions(-) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index fd1352050..0ac909bac 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -773,6 +773,20 @@ @article{vanderstraeten2019 abstract = {In these lecture notes we give a technical overview of tangent-space methods for matrix product states in the thermodynamic limit. We introduce the manifold of uniform matrix product states, show how to compute different types of observables, and discuss the concept of a tangent space. We explain how to variationally optimize ground-state approximations, implement real-time evolution and describe elementary excitations for a given model Hamiltonian. Also, we explain how matrix product states approximate fixed points of one-dimensional transfer matrices. We show how all these methods can be translated to the language of continuous matrix product states for one-dimensional field theories. We conclude with some extensions of the tangent-space formalism and with an outlook to new applications.} } +@article{vanhecke2021cluster, + title = {Symmetric Cluster Expansions with Tensor Networks}, + author = {Vanhecke, Bram and Vanderstraeten, Laurens and Verstraete, Frank}, + year = {2021}, + month = feb, + journal = {Physical Review A}, + volume = {103}, + number = {2}, + pages = {L020402}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevA.103.L020402}, + url = {https://link.aps.org/doi/10.1103/PhysRevA.103.L020402} +} + @article{vanhecke2021, title = {Tangent-Space Methods for Truncating Uniform {{MPS}}}, author = {Vanhecke, Bram and Van Damme, Maarten and Haegeman, Jutho and Vanderstraeten, Laurens and Verstraete, Frank}, diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index 3d2ddc06c..c17f34636 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -1,29 +1,52 @@ """ - ClusterExpansion(N; tol = 1.0e-12) - ClusterExpansion(; N = 2, tol = 1.0e-12) +$(TYPEDEF) -Construct a finite or infinite time evolution MPO for a nearest-neighbor Hamiltonian by -matching exact exponentials on clusters of up to `N` sites. This is the nonperturbative construction of -[Vanhecke, Vanderstraeten and Verstraete (2021)](https://doi.org/10.1103/PhysRevA.103.L020402). +Algorithm for constructing a finite or infinite time evolution MPO using the nonperturbative cluster expansion. +Exact exponentials are matched on clusters of up to `N` sites. For nearest-neighbor Hamiltonians the error is `O(dt^N)`. -Sites may have different physical spaces of the same TensorKit space type. -Symmetry sectors are preserved throughout the construction. Infinite Hamiltonians -match every translated cluster, including those crossing the unit-cell boundary. -Finite Hamiltonians match every cluster contained in the chain and return an -open-boundary MPO. For a chain of length `L`, the cluster size is capped at `L`; -when `N ≥ L`, the entire finite-chain exponential is matched up to roundoff. - -For full-rank uniform physical spaces of dimension `d`, the virtual dimension grows -as `1 + d^2 + ⋯ + d^(2 floor(N/2))`. Mixed physical spaces and rank-deficient residuals -can require additional complementary virtual channels. Exact cluster exponentials -also grow exponentially with `N`. `tol` is the relative singular-value cutoff used -to identify shared SVD support; components below the cutoff are retained in the -complementary channels. Environment equations are solved directly without a -second singular-value cutoff. +# Fields + +$(TYPEDFIELDS) + +# Algorithm + +Exact evolution within clusters of up to `N` sites is included to all orders in +`dt`, together with products of disjoint clusters. The remaining error comes from +connected contributions extending beyond `N` sites: for nearest-neighbor +Hamiltonians the one-step error is `O(dt^N)`. Infinite Hamiltonians include clusters +crossing the unit-cell boundary; finite chains become exact up to roundoff when +`N` reaches the chain length. + +Compared with [`TaylorCluster`](@ref), which matches a specified Taylor order, +`ClusterExpansion` controls the spatial extent of the exact evolution. Its `N` +is a cluster size, not a Taylor order. Resumming smaller clusters can reduce the +error prefactor and permit larger time steps. Compared with the compact, local +construction of [`WII`](@ref), increasing `N` systematically incorporates larger +clusters, at the cost of exponentially growing cluster exponentials and MPO bond +dimensions. For full-rank uniform physical spaces of dimension `d`, the latter +scale as `1 + d^2 + ⋯ + d^(2 floor(N/2))`; mixed spaces or rank-deficient residuals +can require extra channels. + +This is useful for evolution with fewer, larger steps and for growing the MPS bond +dimension from weakly entangled states. Larger steps can also increase compression +costs, so check convergence and runtime by varying `dt` and `N`. MPS compression +introduces a separate error, and the approximate MPO need not be unitary. `tol` +only identifies shared SVD support: smaller singular components are retained in +complementary channels rather than discarded. + +# See also + +Used as the `algorithm` argument of [`make_time_mpo`](@ref). + +# References + +* [Vanhecke et al. Phys. Rev. A 103 (2021)](@cite vanhecke2021cluster) """ struct ClusterExpansion <: Algorithm + "maximum cluster size (positive integer, default: `2`)" N::Int + "relative singular-value cutoff for shared SVD support (`0 < tol < 1`, default: `1.0e-12`)" tol::Float64 function ClusterExpansion(N::Integer; tol::Real = 1.0e-12) N ≥ 1 || throw(ArgumentError("cluster size must be positive")) @@ -38,35 +61,32 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) # contraction shapes are determined by this Val parameter. isempty(H) && throw(ArgumentError("Hamiltonian must contain at least one site")) N = isfinite(H) ? min(alg.N, length(H)) : alg.N - return make_cluster_mpo(H, -im * dt, Val(N), alg.tol) + return make_cluster_expansion_mpo(H, -im * dt, Val(N), alg.tol) end -function make_cluster_mpo(H::MPOHamiltonian, τ::Number, n::Val, tol::Real) +function make_cluster_expansion_mpo(H::MPOHamiltonian, τ::Number, n::Val, tol::Real) T = promote_type(scalartype(H), typeof(τ)) storage = TensorKit.similarstoragetype(storagetype(H), T) O = map(parent(H)) do h - return cluster_initial_tensor(h, τ, storage) + return exponentiate_onsite_tensor(h, τ, storage) end # These cases are already exact products; avoid building channels from # roundoff in differences of independently contracted exponentials. - (iszero(τ) || all(cluster_is_onsite, parent(H))) && return cluster_finalize(O) - return cluster_finalize(cluster_build(O, H, τ, n, tol)) + (iszero(τ) || all(has_only_onsite_terms, parent(H))) && return assemble_evolution_mpo(O) + return assemble_evolution_mpo(add_cluster_corrections(O, H, τ, n, tol)) end -cluster_is_onsite(h) = size(h, 1) ≤ 2 && size(h, 4) ≤ 2 +has_only_onsite_terms(h) = size(h, 1) ≤ 2 && size(h, 4) ≤ 2 -cluster_finalize(O::PeriodicVector) = InfiniteMPO(O) -function cluster_finalize(O::Vector) +assemble_evolution_mpo(O::PeriodicVector) = InfiniteMPO(O) +function assemble_evolution_mpo(O::Vector) O[1] = O[1][1, :, :, :] O[end] = O[end][:, :, :, 1] # Remove structurally unused levels only after all stages have finished. return remove_orphans!(FiniteMPO(O); tol = 0) end -cluster_starts(H::InfiniteMPOHamiltonian, ::Val) = 1:length(H) -cluster_starts(H::FiniteMPOHamiltonian, ::Val{N}) where {N} = 1:(length(H) - N + 1) - -function cluster_initial_tensor(h, τ, ::Type{A}) where {A} +function exponentiate_onsite_tensor(h, τ, ::Type{A}) where {A} P = physicalspace(h) levels = ⊞(oneunit(P)) blocktype = tensormaptype(spacetype(P), 2, 2, A) @@ -76,23 +96,23 @@ function cluster_initial_tensor(h, τ, ::Type{A}) where {A} return O end -cluster_build(O, H, τ, ::Val{1}, tol) = O -function cluster_build(O, H, τ, n::Val{N}, tol) where {N} - previous = cluster_build(O, H, τ, Val(N - 1), tol) - starts = cluster_starts(H, n) +add_cluster_corrections(O, H, τ, ::Val{1}, tol) = O +function add_cluster_corrections(O, H, τ, n::Val{N}, tol) where {N} + previous = add_cluster_corrections(O, H, τ, Val(N - 1), tol) + starts = isfinite(H) ? (1:(length(H) - N + 1)) : (1:length(H)) nt = (N - 1) ÷ 2 level = N ÷ 2 + 1 # Every translated window must see the completed previous stage. Compute # all corrections before inserting any of them, including across the cell seam. if isodd(N) - centers = [cluster_correction(H, previous, s, τ, n) for s in starts] + centers = [solve_center_correction(previous, s, evolution_cluster_residual(H, previous, s, τ, n), n) for s in starts] for s in starts previous[s + nt][level, 1, 1, level] = centers[s] end return previous else - pairs = [cluster_pair(H, previous, s, τ, n, tol) for s in starts] - expanded = cluster_expand(previous, pairs, n) + pairs = [factor_two_site_correction(H, previous, s, τ, n, tol) for s in starts] + expanded = expand_correction_bonds(previous, pairs, n) for s in starts expanded[s + nt][level - 1, 1, 1, level] = pairs[s][1] expanded[s + nt + 1][level, 1, 1, level - 1] = pairs[s][2] @@ -101,13 +121,13 @@ function cluster_build(O, H, τ, n::Val{N}, tol) where {N} end end -function cluster_expand(O::PeriodicVector, pairs, ::Val{N}) where {N} +function expand_correction_bonds(O::PeriodicVector, pairs, ::Val{N}) where {N} nt = (N - 1) ÷ 2 bonds = PeriodicArray([right_virtualspace(pairs[mod1(i - nt, length(O))][1]) for i in 1:length(O)]) - return PeriodicArray([cluster_expand_tensor(O[i], bonds[i - 1], bonds[i]) for i in 1:length(O)]) + return PeriodicArray([append_virtual_spaces(O[i], bonds[i - 1], bonds[i]) for i in 1:length(O)]) end -function cluster_expand(O::Vector, pairs, ::Val{N}) where {N} +function expand_correction_bonds(O::Vector, pairs, ::Val{N}) where {N} nt = (N - 1) ÷ 2 empty_bond = zero(left_virtualspace(O[1])[1]) # Only the middle bonds can support an N-site cluster. Empty levels keep @@ -117,10 +137,10 @@ function cluster_expand(O::Vector, pairs, ::Val{N}) where {N} [right_virtualspace(pair[1]) for pair in pairs], fill(empty_bond, nt + 1) ) - return [cluster_expand_tensor(O[i], bonds[i], bonds[i + 1]) for i in 1:length(O)] + return [append_virtual_spaces(O[i], bonds[i], bonds[i + 1]) for i in 1:length(O)] end -function cluster_expand_tensor(O, left_bond, right_bond) +function append_virtual_spaces(O, left_bond, right_bond) L, R = left_virtualspace(O), right_virtualspace(O) left_levels = L ⊞ left_bond right_levels = R ⊞ right_bond @@ -132,28 +152,22 @@ function cluster_expand_tensor(O, left_bond, right_bond) return expanded end -function cluster_residual(H, O, start::Int, τ, n::Val{N}) where {N} +function evolution_cluster_residual(H, O, start::Int, τ, n::Val{N}) where {N} boundary = size(H[start + N - 1], 4) - exact = exp(τ * cluster_contract(H, start, n, boundary)) - return add_util_leg(exact - cluster_contract(O, start, n, 1)) + exact = exp(τ * contract_mpo_window(H, start, n, boundary)) + return add_util_leg(exact - contract_mpo_window(O, start, n, 1)) end -function cluster_correction(H, O, start, τ, n::Val) - residual = cluster_residual(H, O, start, τ, n) - return cluster_center(O, start, residual, n) -end - -function cluster_pair(H, O, start, τ, n::Val{2}, tol) - residual = cluster_residual(H, O, start, τ, n) - return cluster_split(residual, tol) +function factor_two_site_correction(H, O, start, τ, n::Val{2}, tol) + residual = evolution_cluster_residual(H, O, start, τ, n) + A, B = factor_with_complementary_spaces(permute(residual, ((1, 2, 4), (3, 5, 6))), tol) + return permute(A, ((1, 2), (3, 4))), permute(B, ((1, 2), (3, 4))) end -function cluster_pair(H, O, start, τ, n::Val{N}, tol) where {N} - residual = cluster_residual(H, O, start, τ, n) - nt = (N - 1) ÷ 2 - left, right = cluster_environments(O, start, n) - left_map, right_map = cluster_environment_maps(left, right, Val(nt)) - center = cluster_center(left_map, right_map, residual, n) +function factor_two_site_correction(H, O, start, τ, n::Val{N}, tol) where {N} + residual = evolution_cluster_residual(H, O, start, τ, n) + left_map, right_map = correction_environment_maps(O, start, n) + center = solve_center_correction(left_map, right_map, residual, n) C = permute(center, ((1, 2, 4), (3, 5, 6))) # The minimum-norm solves restrict the center to the active environment @@ -164,18 +178,13 @@ function cluster_pair(H, O, start, τ, n::Val{N}, tol) where {N} Q_right, _ = left_orth(right_map) F_left = Q_left' ⊗ id(storagetype(C), codomain(C)[2] ⊗ codomain(C)[3]) F_right = id(storagetype(C), domain(C)[1] ⊗ domain(C)[2]) ⊗ Q_right - A, B = cluster_complete_svd(F_left' * C * F_right, tol) + A, B = factor_with_complementary_spaces(F_left' * C * F_right, tol) return permute(F_left * A, ((1, 2), (3, 4))), permute(B * F_right', ((1, 2), (3, 4))) end -function cluster_split(center, tol) - A, B = cluster_complete_svd(permute(center, ((1, 2, 4), (3, 5, 6))), tol) - return permute(A, ((1, 2), (3, 4))), permute(B, ((1, 2), (3, 4))) -end - """ - cluster_complete_svd(C, tol) + factor_with_complementary_spaces(C, tol) Factor `C` through its shared SVD support and separate left/right complementary spaces. For exact rank `r_c`, the virtual multiplicity in sector `c` is @@ -188,7 +197,7 @@ At a finite rank cutoff, retain the remaining block `E = U⊥' C V⊥` in the complementary channels, so the product still reconstructs `C` up to roundoff. The scale `γ` balances these channels against the retained singular values. """ -function cluster_complete_svd(C::TensorMap, tol::Real) +function factor_with_complementary_spaces(C::TensorMap, tol::Real) cnorm = norm(C) if iszero(cnorm) # At zero rank, the two complementary spaces occupy disjoint channels. @@ -212,24 +221,24 @@ function cluster_complete_svd(C::TensorMap, tol::Real) return A, B end -function cluster_contract(O, start::Int, ::Val{1}, boundary::Int) +function contract_mpo_window(O, start::Int, ::Val{1}, boundary::Int) return removeunit(removeunit(TensorMap(O[start][1, :, :, boundary]), 4), 1) end -function cluster_contract(O, start::Int, n::Val{N}, boundary::Int) where {N} +function contract_mpo_window(O, start::Int, n::Val{N}, boundary::Int) where {N} left = removeunit(TensorMap(O[start][1, :, :, :]), 1) right = removeunit(TensorMap(O[start + N - 1][:, :, :, boundary]), 4) - sites = cluster_bulk(O, start + 1, Val(N - 2)) - return cluster_contract(left, sites, right, n) + sites = collect_bulk_tensors(O, start + 1, Val(N - 2)) + return contract_mpo_window(left, sites, right, n) end -cluster_bulk(O, start, ::Val{0}) = () -function cluster_bulk(O, start, ::Val{N}) where {N} - return (TensorMap(O[start]), cluster_bulk(O, start + 1, Val(N - 1))...) +collect_bulk_tensors(O, start, ::Val{0}) = () +function collect_bulk_tensors(O, start, ::Val{N}) where {N} + return (TensorMap(O[start]), collect_bulk_tensors(O, start + 1, Val(N - 1))...) end """ - cluster_contract(left, sites, right, ::Val{N}) + contract_mpo_window(left, sites, right, ::Val{N}) Contract the open-boundary cluster, shown for `N = 4`. `L` and `R` are the boundary tensors; `O₂` and `O₃` are the site-dependent bulk tensors. Negative labels are open @@ -246,7 +255,7 @@ are above the tensors and codomain legs below. -1 -2 -3 -4 ``` """ -@generated function cluster_contract( +@generated function contract_mpo_window( left::AbstractTensorMap{<:Any, S, 1, 2}, sites::NTuple{K, AbstractTensorMap{<:Any, S, 2, 2}}, right::AbstractTensorMap{<:Any, S, 2, 1}, ::Val{N} ) where {S, N, K} @@ -257,26 +266,20 @@ are above the tensors and codomain legs below. return macroexpand(@__MODULE__, :(return @plansor $out := *($first_site, $(middle...), $last_site))) end -function cluster_environments(O, start, ::Val{N}) where {N} - nt = (N - 1) ÷ 2 - return cluster_left_environment(O, start, Val(nt)), - cluster_right_environment(O, start + N - nt, Val(nt)) -end - -cluster_left_environment(O, start, ::Val{1}) = O[start][1, 1, 1, 2] -function cluster_left_environment(O, start, ::Val{L}) where {L} - left = cluster_left_environment(O, start, Val(L - 1)) - return cluster_join_left(left, O[start + L - 1][L, 1, 1, L + 1]) +left_correction_environment(O, start, ::Val{1}) = O[start][1, 1, 1, 2] +function left_correction_environment(O, start, ::Val{L}) where {L} + left = left_correction_environment(O, start, Val(L - 1)) + return extend_left_environment(left, O[start + L - 1][L, 1, 1, L + 1]) end -cluster_right_environment(O, start, ::Val{1}) = O[start][2, 1, 1, 1] -function cluster_right_environment(O, start, ::Val{L}) where {L} - right = cluster_right_environment(O, start + 1, Val(L - 1)) - return cluster_join_right(O[start][L + 1, 1, 1, L], right) +right_correction_environment(O, start, ::Val{1}) = O[start][2, 1, 1, 1] +function right_correction_environment(O, start, ::Val{L}) where {L} + right = right_correction_environment(O, start + 1, Val(L - 1)) + return extend_right_environment(O[start][L + 1, 1, 1, L], right) end """ - cluster_join_left(left, site) + extend_left_environment(left, site) Append a site to the left environment, shown for `K = numout(left) = 2`. For general `K`, each vertical leg of `left` is a bundle of `K - 1` physical @@ -293,15 +296,15 @@ domain legs are above the tensors and codomain legs below. -2 -3 ``` """ -@generated function cluster_join_left(left::AbstractTensorMap{<:Any, <:Any, K, K}, site::AbstractTensorMap) where {K} +@generated function extend_left_environment(left::AbstractTensorMap{<:Any, <:Any, K, K}, site::AbstractTensorMap) where {K} out = tensorexpr(:next_left, -(1:(K + 1)), -((K + 2):(2K + 2))) - env = tensorexpr(:left, -(1:K), (Tuple(-((K + 2):(2K)))..., 1)) + env = tensorexpr(:left, -(1:K), (-((K + 2):(2K))..., 1)) op = tensorexpr(:site, (1, -K - 1), (-2K - 1, -2K - 2)) return macroexpand(@__MODULE__, :(return @plansor $out := $env * $op)) end """ - cluster_join_right(site, right) + extend_right_environment(site, right) Prepend a site to the right environment, shown for `K = numout(right) = 2`. For general `K`, each vertical leg of `right` is a bundle of `K - 1` physical @@ -318,27 +321,29 @@ domain legs are above the tensors and codomain legs below. -2 -3 ``` """ -@generated function cluster_join_right(site::AbstractTensorMap, right::AbstractTensorMap{<:Any, <:Any, K, K}) where {K} +@generated function extend_right_environment(site::AbstractTensorMap, right::AbstractTensorMap{<:Any, <:Any, K, K}) where {K} out = tensorexpr(:next_right, -(1:(K + 1)), -((K + 2):(2K + 2))) op = tensorexpr(:site, (-1, -2), (-K - 2, 1)) - env = tensorexpr(:right, (1, Tuple(-(3:(K + 1)))...), -((K + 3):(2K + 2))) + env = tensorexpr(:right, (1, -(3:(K + 1))...), -((K + 3):(2K + 2))) return macroexpand(@__MODULE__, :(return @plansor $out := $op * $env)) end -function cluster_environment_maps(left, right, ::Val{L}) where {L} +function correction_environment_maps(O, start, ::Val{N}) where {N} + L = (N - 1) ÷ 2 + left = left_correction_environment(O, start, Val(L)) + right = right_correction_environment(O, start + N - L, Val(L)) indices = ntuple(identity, Val(2L + 1)) return permute(left, (indices, (2L + 2,))), permute(right, ((1,), indices .+ 1)) end -function cluster_center(O, start::Int, residual::TensorMap, n::Val{N}) where {N} - left, right = cluster_environments(O, start, n) - left_map, right_map = cluster_environment_maps(left, right, Val((N - 1) ÷ 2)) - return cluster_center(left_map, right_map, residual, n) +function solve_center_correction(O, start::Int, residual::TensorMap, n::Val{N}) where {N} + left_map, right_map = correction_environment_maps(O, start, n) + return solve_center_correction(left_map, right_map, residual, n) end """ - cluster_center(left_map, right_map, residual, ::Val{N}) + solve_center_correction(left_map, right_map, residual, ::Val{N}) Solve the environment equations from the left and right without forming inverse maps. Completed environments have full row/column support in each symmetry @@ -376,17 +381,17 @@ Finally move `p₂′` back into the domain. The center has one physical leg per side for odd `N`, or two for even `N`, and its two virtual legs. For larger clusters, each outer physical space above is a bundle of `(N - 1) ÷ 2` legs. """ -@generated function cluster_center( +@generated function solve_center_correction( left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N} ) where {N} nt = (N - 1) ÷ 2 width = N - 2nt # Group all left environment legs into the codomain for the left solve. - left_indices = (1, Tuple(2:(nt + 1))..., Tuple((N + 2):(N + nt + 1))...) - remaining = (Tuple((nt + 2):(N + 1))..., Tuple((N + nt + 2):(2N + 2))...) + left_indices = (1, (2:(nt + 1))..., ((N + 2):(N + nt + 1))...) + remaining = (((nt + 2):(N + 1))..., ((N + nt + 2):(2N + 2))...) # After solving, group the right environment legs into the domain. - middle_indices = (1, Tuple(2:(width + 1))..., Tuple((nt + width + 2):(nt + 2width + 1))...) - right_indices = (Tuple((width + 2):(nt + width + 1))..., Tuple((nt + 2width + 2):(2nt + 2width + 2))...) + middle_indices = (1, (2:(width + 1))..., ((nt + width + 2):(nt + 2width + 1))...) + right_indices = (((width + 2):(nt + width + 1))..., ((nt + 2width + 2):(2nt + 2width + 2))...) center_indices = (Tuple(1:(width + 1)), Tuple((width + 2):(2width + 2))) return quote left_rhs = permute(residual, $((left_indices, remaining))) diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl index 664d3425e..f14f11cf9 100644 --- a/test/timeevolution/clusterexpansion.jl +++ b/test/timeevolution/clusterexpansion.jl @@ -4,14 +4,14 @@ using TensorKit: ℙ # Construct finite windows independently of the cluster-contraction kernels, # including lengths that are not multiples of the infinite unit cell. -function cluster_dense_window(U, start, L) +function dense_mpo_window(U, start, L) sites = [copy(U[start + i - 1]) for i in 1:L] sites[1] = sites[1][1, :, :, :] sites[end] = sites[end][:, :, :, 1] return convert(TensorMap, MPSKit.DenseMPO(FiniteMPO(sites))) end -function cluster_exact_window(H, τ, start, L) +function exact_evolution_window(H, τ, start, L) sites = [copy(H[start + i - 1]) for i in 1:L] sites[1] = sites[1][1, :, :, :] sites[end] = sites[end][:, :, :, end] @@ -91,23 +91,23 @@ end U = make_time_mpo(symmetric_H, 0.1, ClusterExpansion(1)) @test dense(U, 1) ≈ exp(-0.1im * onsite) atol = 1.0e-10 for N in 2:5 - U = @inferred MPSKit.make_cluster_mpo(symmetric_H, -0.1im, Val(N), 1.0e-12) + U = @inferred MPSKit.make_cluster_expansion_mpo(symmetric_H, -0.1im, Val(N), 1.0e-12) @test physicalspace(U, 1) == P @test left_virtualspace(U, 1)[2] == fuse(P ⊗ dual(P)) exact_cluster = exp(-0.1im * convert(TensorMap, open_boundary_conditions(symmetric_H, N))) @test dense(U, N) ≈ exact_cluster atol = 1.0e-10 end - residual = @inferred MPSKit.cluster_residual(symmetric_H, U, 1, -0.1im, Val(5)) - center = @inferred MPSKit.cluster_center(U, 1, residual, Val(5)) + residual = @inferred MPSKit.evolution_cluster_residual(symmetric_H, U, 1, -0.1im, Val(5)) + center = @inferred MPSKit.solve_center_correction(U, 1, residual, Val(5)) @test center isa TensorMap end end @testset "Inference after the cluster-size dispatch" begin for n in (Val(4), Val(5)) - U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, n, 1.0e-12) - residual = @inferred MPSKit.cluster_residual(H, U, 1, -0.1im, n) - center = @inferred MPSKit.cluster_center(U, 1, residual, n) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, n, 1.0e-12) + residual = @inferred MPSKit.evolution_cluster_residual(H, U, 1, -0.1im, n) + center = @inferred MPSKit.solve_center_correction(U, 1, residual, n) @test center isa TensorMap end end @@ -121,7 +121,7 @@ end U = make_time_mpo(H, 0.1, ClusterExpansion(N)) U2 = make_time_mpo(H2, 0.1, ClusterExpansion(N)) @test length(U2) == 2 - @test cluster_dense_window(U2, 2, N + 1) ≈ dense(U, N + 1) atol = 1.0e-10 + @test dense_mpo_window(U2, 2, N + 1) ≈ dense(U, N + 1) atol = 1.0e-10 end end @@ -147,13 +147,13 @@ end end H = InfiniteMPOHamiltonian(Ps, terms...) for N in 1:5 - U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), 1.0e-12) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) @test length(U) == length(H) for s in eachindex(lattice) @test physicalspace(U, s) == Ps[s] @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) for L in 1:N - @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(H, -0.1im, s, L) atol = 1.0e-10 + @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(H, -0.1im, s, L) atol = 1.0e-10 end end end @@ -173,7 +173,7 @@ end raw = randn!(rng, zeros(ComplexF64, L ← R)) for scale in (1.0, 1.0e-28, 0.0), tol in (1.0e-12, 0.5) C = scale * raw - A, B = @inferred MPSKit.cluster_complete_svd(C, tol) + A, B = @inferred MPSKit.factor_with_complementary_spaces(C, tol) @test norm(A * B - C) ≤ 2.0e-12 * max(norm(C), 1) @test domain(A) == codomain(B) @test A * pinv(A; rtol = 1.0e-12) ≈ id(storagetype(C), L) atol = 1.0e-10 @@ -188,7 +188,7 @@ end end # The SVD rank cutoff must not discard the tiny second singular value. C = TensorMap(ComplexF64[1 0 0; 0 1.0e-14 0; 0 0 0], ℂ^3 ← ℂ^3) - A, B = MPSKit.cluster_complete_svd(C, 1.0e-12) + A, B = MPSKit.factor_with_complementary_spaces(C, 1.0e-12) @test dim(only(domain(A))) == 5 @test norm(A * B - C) < 1.0e-16 end @@ -213,14 +213,14 @@ end U = make_time_mpo(model, dt, ClusterExpansion(5)) longest = model === onsite || model === dimers || iszero(dt) ? 6 : 5 for s in eachindex(lattice), L in 1:longest - @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(model, -im * dt, s, L) atol = 1.0e-10 + @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(model, -im * dt, s, L) atol = 1.0e-10 end end # A large cutoff changes the SVD channels, not the matched exponentials. for N in (3, 5), tol in (0.1, 0.5) - U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), tol) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), tol) for s in eachindex(lattice), L in 1:N - @test cluster_dense_window(U, s, L) ≈ cluster_exact_window(H, -0.1im, s, L) atol = 1.0e-10 + @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(H, -0.1im, s, L) atol = 1.0e-10 end end for N in 3:4, s in eachindex(lattice) @@ -229,7 +229,7 @@ end errors = Float64[] for dt in steps U = make_time_mpo(H, dt, ClusterExpansion(N)) - push!(errors, norm(cluster_dense_window(U, s, N + 1) - cluster_exact_window(H, -im * dt, s, N + 1))) + push!(errors, norm(dense_mpo_window(U, s, N + 1) - exact_evolution_window(H, -im * dt, s, N + 1))) end @test log2(errors[1] / errors[2]) ≥ N - 0.15 end @@ -238,7 +238,7 @@ end # Independent random couplings and onsite terms exercise nonuniform finite chains. -function cluster_finite_model(rng, Ps; onsite_only = false, dimers = false) +function finite_test_hamiltonian(rng, Ps; onsite_only = false, dimers = false) terms = Pair[] for s in eachindex(Ps) a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) @@ -262,10 +262,10 @@ end [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1), Vect[FermionParity](0 => 1, 1 => 1)], ) for Ps in lattices - H = cluster_finite_model(rng, Ps) + H = finite_test_hamiltonian(rng, Ps) L = length(H) for N in 1:min(6, L) - U = @inferred MPSKit.make_cluster_mpo(H, -0.1im, Val(N), 1.0e-12) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) @test U isa FiniteMPO @test length(U) == L @test dim(left_virtualspace(U, 1)) == dim(right_virtualspace(U, L)) == 1 @@ -273,7 +273,7 @@ end @test physicalspace(U, s) == Ps[s] s < L && @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) for width in 1:min(N, L - s + 1) - @test cluster_dense_window(U, s, width) ≈ cluster_exact_window(H, -0.1im, s, width) atol = 1.0e-10 + @test dense_mpo_window(U, s, width) ≈ exact_evolution_window(H, -0.1im, s, width) atol = 1.0e-10 end end end @@ -288,9 +288,9 @@ end @testset "Finite edge cases and leading error" begin rng = MersenneTwister(984) Ps = [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2] - H = cluster_finite_model(rng, Ps) - onsite = cluster_finite_model(rng, Ps; onsite_only = true) - dimers = cluster_finite_model(rng, Ps; dimers = true) + H = finite_test_hamiltonian(rng, Ps) + onsite = finite_test_hamiltonian(rng, Ps; onsite_only = true) + dimers = finite_test_hamiltonian(rng, Ps; dimers = true) for (model, dt) in ((H, 0.0), (H, -0.1), (H, -0.1im), (onsite, 0.1), (dimers, 0.1)) U = make_time_mpo(model, dt, ClusterExpansion(5)) @test convert(TensorMap, MPSKit.DenseMPO(U)) ≈ exp(-im * dt * convert(TensorMap, model)) atol = 1.0e-10 From 50ebce405b6d1473166517201a6ad5d636a27438 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 13:36:01 -0400 Subject: [PATCH 5/8] Streamline cluster expansion tests and document usage --- docs/src/man/algorithms.md | 20 ++ test/timeevolution/clusterexpansion.jl | 355 +++++++------------------ 2 files changed, 117 insertions(+), 258 deletions(-) diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 89c0d9ab2..5d4029c87 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -128,6 +128,26 @@ ClusterExpansion See also [Time evolution accuracy](@ref). +#### Nonperturbative cluster expansion + +[`ClusterExpansion`](@ref) constructs evolution MPOs for finite and infinite +nearest-neighbor Hamiltonians by matching exact exponentials on clusters of up to +`N` sites [vanhecke2021cluster](@cite). Unlike [`TaylorCluster`](@ref), `N` specifies +the cluster size rather than a Taylor order. For nearest-neighbor interactions, +the one-step error is `O(dt^N)`; a finite chain is matched up to roundoff when +`N` reaches its length. + +```julia +U = make_time_mpo(H, dt, ClusterExpansion(; N = 3)) +ψ, envs, info = approximate(ψ, (U, ψ); trunc = truncrank(32)) +``` + +Matching small clusters to all orders can permit larger time steps and help grow +the bond dimension from weakly entangled states. The construction costs more than +the compact [`WII`](@ref) approximation, and increasing `N` requires exponentially +larger cluster exponentials and MPO bonds. Check convergence by varying `dt` and +`N`, as well as the MPS compression settings. + ## Excitations It might also be desirable to obtain information beyond the lowest energy state of a given system, and study the dispersion relation. diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl index f14f11cf9..324044ce4 100644 --- a/test/timeevolution/clusterexpansion.jl +++ b/test/timeevolution/clusterexpansion.jl @@ -1,317 +1,156 @@ using Test using MPSKit, TensorKit, LinearAlgebra, Random using TensorKit: ℙ +using TensorKitTensors.SpinOperators: σˣ, σʸ, σᶻ -# Construct finite windows independently of the cluster-contraction kernels, -# including lengths that are not multiples of the infinite unit cell. -function dense_mpo_window(U, start, L) - sites = [copy(U[start + i - 1]) for i in 1:L] +# Independent dense windows also cover lengths that are not multiples of the unit cell. +function dense_mpo_window(O, start = 1, L = length(O)) + sites = [copy(O[start + i - 1]) for i in 1:L] sites[1] = sites[1][1, :, :, :] - sites[end] = sites[end][:, :, :, 1] - return convert(TensorMap, MPSKit.DenseMPO(FiniteMPO(sites))) + boundary = O isa MPOHamiltonian ? size(sites[end], 4) : 1 + sites[end] = sites[end][:, :, :, boundary] + window = O isa MPOHamiltonian ? FiniteMPOHamiltonian(sites) : MPSKit.DenseMPO(FiniteMPO(sites)) + return convert(TensorMap, window) end -function exact_evolution_window(H, τ, start, L) - sites = [copy(H[start + i - 1]) for i in 1:L] - sites[1] = sites[1][1, :, :, :] - sites[end] = sites[end][:, :, :, end] - return exp(τ * convert(TensorMap, FiniteMPOHamiltonian(sites))) +function random_test_hamiltonian(rng, Ps; onsite_only = false, dimers = false) + infinite = Ps isa MPSKit.PeriodicVector + terms = Pair[] + for s in eachindex(Ps) + a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) + push!(terms, s => (a + a') / 2) + if (infinite || s < length(Ps)) && !onsite_only && (!dimers || isodd(s)) + b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) + push!(terms, (s, s + 1) => (b + b') / 2) + end + end + return infinite ? InfiniteMPOHamiltonian(Ps, terms...) : FiniteMPOHamiltonian(Ps, terms...) end @testset "Nonperturbative cluster expansion" begin - P = ℂ^2 - X = TensorMap(ComplexF64[0 1; 1 0], P ← P) - Z = TensorMap(ComplexF64[1 0; 0 -1], P ← P) - h = Z ⊗ Z + 0.3 * (X ⊗ X) + X, Y, Z = σˣ(), σʸ(), σᶻ() + P = space(X, 1) + h = Z ⊗ Z + 0.3 * X ⊗ X + 0.2 * (X ⊗ Y + Y ⊗ X) H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.7X, (1, 2) => h) - dense(U, L) = convert(TensorMap, MPSKit.DenseMPO(open_boundary_conditions(U, L))) - exact(H, dt, L) = L == 1 ? exp(-im * dt * 0.7X) : - exp(-im * dt * convert(TensorMap, open_boundary_conditions(H, L))) - @testset "Exact clusters, N=$N" for N in 1:5 - U = make_time_mpo(H, 0.1, ClusterExpansion(N)) + @testset "Exact clusters and inference, N=$N" for N in 1:5 + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) @test U isa InfiniteMPO @test dim(left_virtualspace(U, 1)) == sum(4^l for l in 0:(N ÷ 2)) for L in 1:N - @test dense(U, L) ≈ exact(H, 0.1, L) atol = 1.0e-10 - end - end - - @testset "Time steps and rank-deficient environments" begin - onsite = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.7X) - commuting = InfiniteMPOHamiltonian(PeriodicArray([P]), (1, 2) => Z ⊗ Z) - for dt in (0.0, -0.1, -0.1im) - U = make_time_mpo(H, dt, ClusterExpansion(3)) - @test dense(U, 3) ≈ exact(H, dt, 3) atol = 1.0e-10 + @test dense_mpo_window(U, 1, L) ≈ exp(-0.1im * dense_mpo_window(H, 1, L)) atol = 1.0e-10 end - for Hsimple in (onsite, commuting) - U = make_time_mpo(Hsimple, 0.1, ClusterExpansion(3)) - @test dense(U, 4) ≈ exp(-0.1im * convert(TensorMap, open_boundary_conditions(Hsimple, 4))) atol = 1.0e-10 + if N ≥ 3 + residual = @inferred MPSKit.evolution_cluster_residual(H, U, 1, -0.1im, Val(N)) + @inferred MPSKit.solve_center_correction(U, 1, residual, Val(N)) end - end - - @testset "Complex Hamiltonian and reflection symmetry" begin - Y = TensorMap(ComplexF64[0 -im; im 0], P ← P) - complex_h = h + 0.2 * (X ⊗ Y + Y ⊗ X) - complex_H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => 0.4Y, (1, 2) => complex_h) - for N in (3, 4), dt in (0.1, -0.1im) - U = make_time_mpo(complex_H, dt, ClusterExpansion(N)) - @test dense(U, N) ≈ exp(-im * dt * convert(TensorMap, open_boundary_conditions(complex_H, N))) atol = 1.0e-10 - larger = dense(U, N + 1) + if N in (3, 4) + larger = dense_mpo_window(U, 1, N + 1) reflected = permute(larger, (Tuple((N + 1):-1:1), Tuple((2N + 2):-1:(N + 2)))) @test larger ≈ reflected atol = 1.0e-10 + U2 = make_time_mpo(repeat(H, 2), 0.1, ClusterExpansion(N)) + @test dense_mpo_window(U2, 2, N + 1) ≈ larger atol = 1.0e-10 end end - @testset "Leading error on larger chains" begin - for N in 2:4 - errors = Float64[] - for dt in (0.04, 0.02) - U = make_time_mpo(H, dt, ClusterExpansion(N)) - push!(errors, norm(dense(U, N + 1) - exact(H, dt, N + 1))) - end + @testset "Leading error, N=$N" for N in 2:4 + for model in (H, open_boundary_conditions(H, 5)) + errors = [ + norm( + dense_mpo_window(make_time_mpo(model, dt, ClusterExpansion(N)), 1, N + 1) - + exp(-im * dt * dense_mpo_window(model, 1, N + 1)) + ) for dt in (0.04, 0.02) + ] @test log2(errors[1] / errors[2]) ≥ N - 0.15 end end - @testset "Symmetry sectors and dual physical spaces" begin - rng = MersenneTwister(1234) - spaces = ( - ℝ^2, ℙ^2, - Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1)), - Rep[SU₂](1 // 2 => 1), dual(Rep[SU₂](1 // 2 => 1)), - Vect[FermionParity](0 => 1, 1 => 1), - ) - for P in spaces - raw1 = randn!(rng, zeros(ComplexF64, P ← P)) - onsite = (raw1 + raw1') / 2 - raw2 = randn!(rng, zeros(ComplexF64, P ⊗ P ← P ⊗ P)) - interaction = (raw2 + raw2') / 2 - symmetric_H = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => onsite, (1, 2) => interaction) - U = make_time_mpo(symmetric_H, 0.1, ClusterExpansion(1)) - @test dense(U, 1) ≈ exp(-0.1im * onsite) atol = 1.0e-10 - for N in 2:5 - U = @inferred MPSKit.make_cluster_expansion_mpo(symmetric_H, -0.1im, Val(N), 1.0e-12) - @test physicalspace(U, 1) == P - @test left_virtualspace(U, 1)[2] == fuse(P ⊗ dual(P)) - exact_cluster = exp(-0.1im * convert(TensorMap, open_boundary_conditions(symmetric_H, N))) - @test dense(U, N) ≈ exact_cluster atol = 1.0e-10 - end - residual = @inferred MPSKit.evolution_cluster_residual(symmetric_H, U, 1, -0.1im, Val(5)) - center = @inferred MPSKit.solve_center_correction(U, 1, residual, Val(5)) - @test center isa TensorMap - end + @testset "Finite boundaries and cluster-size cap, L=$L" for L in (1, 2, 3, 6) + finite = open_boundary_conditions(H, L) + U = make_time_mpo(finite, 0.1, ClusterExpansion(L + 3; tol = 0.5)) + @test U isa FiniteMPO + @test dim(left_virtualspace(U, 1)) == dim(right_virtualspace(U, L)) == 1 + @test dense_mpo_window(U) ≈ exp(-0.1im * convert(TensorMap, finite)) atol = 1.0e-10 end - - @testset "Inference after the cluster-size dispatch" begin - for n in (Val(4), Val(5)) - U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, n, 1.0e-12) - residual = @inferred MPSKit.evolution_cluster_residual(H, U, 1, -0.1im, n) - center = @inferred MPSKit.solve_center_correction(U, 1, residual, n) - @test center isa TensorMap - end + @testset "Finite/infinite consistency, N=$N" for N in (2, 3, 4, 5) + finite = make_time_mpo(open_boundary_conditions(H, 6), 0.1, ClusterExpansion(N)) + infinite = make_time_mpo(H, 0.1, ClusterExpansion(N)) + @test dense_mpo_window(finite) ≈ dense_mpo_window(infinite, 1, 6) atol = 1.0e-10 end @test_throws ArgumentError ClusterExpansion(0) @test_throws ArgumentError ClusterExpansion(2; tol = 0) @test_throws ArgumentError ClusterExpansion(2; tol = Inf) @test ClusterExpansion(; N = 3).N == 3 - H2 = repeat(H, 2) - @testset "Repeated unit cell" for N in 2:5 - U = make_time_mpo(H, 0.1, ClusterExpansion(N)) - U2 = make_time_mpo(H2, 0.1, ClusterExpansion(N)) - @test length(U2) == 2 - @test dense_mpo_window(U2, 2, N + 1) ≈ dense(U, N + 1) atol = 1.0e-10 - end + empty = FiniteMPOHamiltonian(similar(parent(open_boundary_conditions(H, 2)), 0)) + @test_throws ArgumentError make_time_mpo(empty, 0.1, ClusterExpansion(2)) end -@testset "Periodic cluster expansion with mixed physical spaces" begin +@testset "Mixed physical spaces and boundaries" begin rng = MersenneTwister(25) lattices = ( - [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], - [ℝ^2, ℝ^3], [ℙ^2, ℙ^3], - [Rep[U₁](0 => 1, 1 => 1), Rep[U₁](0 => 1, 1 => 1, 2 => 1)], - [dual(Rep[U₁](0 => 1, 1 => 1)), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1))], + [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], [ℝ^2, ℝ^3], [ℙ^2, ℙ^3], [Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1))], - [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1)], - [dual(Rep[SU₂](1 // 2 => 1)), dual(Rep[SU₂](1 => 1))], + [dual(Rep[SU₂](1 // 2 => 1)), Rep[SU₂](1 => 1)], [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1)], ) - for lattice in lattices - Ps = PeriodicArray(lattice) - terms = Pair[] - for s in eachindex(lattice) - a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) - b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) - push!(terms, s => (a + a') / 2, (s, s + 1) => (b + b') / 2) - end - H = InfiniteMPOHamiltonian(Ps, terms...) - for N in 1:5 + for lattice in lattices, infinite in (true, false) + Ps = infinite ? PeriodicArray(lattice) : repeat(lattice, 3)[1:5] + H = random_test_hamiltonian(rng, Ps) + for N in (4, 5) U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) @test length(U) == length(H) - for s in eachindex(lattice) + for s in eachindex(Ps) @test physicalspace(U, s) == Ps[s] - @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) - for L in 1:N - @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(H, -0.1im, s, L) atol = 1.0e-10 - end + (infinite || s < length(H)) && @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) + width = infinite ? N : min(N, length(H) - s + 1) + @test dense_mpo_window(U, s, width) ≈ exp(-0.1im * dense_mpo_window(H, s, width)) atol = 1.0e-10 end end end end -@testset "Completed SVD support" begin - rng = MersenneTwister(22) - spaces = ( - (ℂ^4, ℂ^9), - (Rep[U₁](0 => 4, 1 => 1), Rep[U₁](0 => 2, -1 => 2, 1 => 3)), - (Rep[SU₂](0 => 3, 1 // 2 => 2, 1 => 1), Rep[SU₂](0 => 1, 1 // 2 => 3, 2 => 1)), - (Vect[FermionParity](0 => 4, 1 => 1), Vect[FermionParity](0 => 2, 1 => 3)), - ) - for (L, R) in spaces, dualize in (false, true) - L, R = dualize ? (dual(L), dual(R)) : (L, R) - raw = randn!(rng, zeros(ComplexF64, L ← R)) - for scale in (1.0, 1.0e-28, 0.0), tol in (1.0e-12, 0.5) - C = scale * raw - A, B = @inferred MPSKit.factor_with_complementary_spaces(C, tol) - @test norm(A * B - C) ≤ 2.0e-12 * max(norm(C), 1) - @test domain(A) == codomain(B) - @test A * pinv(A; rtol = 1.0e-12) ≈ id(storagetype(C), L) atol = 1.0e-10 - @test pinv(B; rtol = 1.0e-12) * B ≈ id(storagetype(C), R) atol = 1.0e-10 - if scale == 1 && tol == 1.0e-12 - V = only(domain(A)) - for c in union(sectors(L), sectors(R)) - @test dim(V, c) == max(dim(L, c), dim(R, c)) - end - end - end - end - # The SVD rank cutoff must not discard the tiny second singular value. - C = TensorMap(ComplexF64[1 0 0; 0 1.0e-14 0; 0 0 0], ℂ^3 ← ℂ^3) - A, B = MPSKit.factor_with_complementary_spaces(C, 1.0e-12) - @test dim(only(domain(A))) == 5 - @test norm(A * B - C) < 1.0e-16 -end - -@testset "Mixed-cell edge cases and leading error" begin +@testset "Time steps and rank-deficient environments" begin rng = MersenneTwister(24) - for lattice in ([ℂ^2, ℂ^3], [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1)]) - Ps = PeriodicArray(lattice) - onsite_terms = Pair[] - full_terms = Pair[] - for s in eachindex(lattice) - a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) - b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) - push!(onsite_terms, s => (a + a') / 2) - push!(full_terms, s => (a + a') / 2, (s, s + 1) => (b + b') / 2) - end - onsite = InfiniteMPOHamiltonian(Ps, onsite_terms...) - H = InfiniteMPOHamiltonian(Ps, full_terms...) - # Remove the interaction across the cell seam to form independent dimers. - dimers = InfiniteMPOHamiltonian(Ps, filter(p -> p.first != (2, 3), full_terms)...) + Z = σᶻ() + P = space(Z, 1) + commuting = InfiniteMPOHamiltonian(PeriodicArray([P]), (1, 2) => Z ⊗ Z) + U = make_time_mpo(commuting, 0.1, ClusterExpansion(4)) + @test dense_mpo_window(U, 1, 5) ≈ exp(-0.1im * dense_mpo_window(commuting, 1, 5)) atol = 1.0e-10 + + for Ps in (PeriodicArray([ℂ^2, ℂ^3]), [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2]) + H = random_test_hamiltonian(rng, Ps) + onsite = random_test_hamiltonian(rng, Ps; onsite_only = true) + dimers = random_test_hamiltonian(rng, Ps; dimers = true) for (model, dt) in ((onsite, 0.1), (dimers, 0.1), (H, 0.0), (H, -0.1im), (H, -0.1)) U = make_time_mpo(model, dt, ClusterExpansion(5)) - longest = model === onsite || model === dimers || iszero(dt) ? 6 : 5 - for s in eachindex(lattice), L in 1:longest - @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(model, -im * dt, s, L) atol = 1.0e-10 - end - end - # A large cutoff changes the SVD channels, not the matched exponentials. - for N in (3, 5), tol in (0.1, 0.5) - U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), tol) - for s in eachindex(lattice), L in 1:N - @test dense_mpo_window(U, s, L) ≈ exact_evolution_window(H, -0.1im, s, L) atol = 1.0e-10 - end - end - for N in 3:4, s in eachindex(lattice) - # Keep the SU(2) N=4 errors above the floating-point floor. - steps = spacetype(Ps[1]) == ComplexSpace || N == 3 ? (0.04, 0.02) : (0.16, 0.08) - errors = Float64[] - for dt in steps - U = make_time_mpo(H, dt, ClusterExpansion(N)) - push!(errors, norm(dense_mpo_window(U, s, N + 1) - exact_evolution_window(H, -im * dt, s, N + 1))) + width = isfinite(model) || model === H && !iszero(dt) ? 5 : 6 + starts = isfinite(model) ? (1:1) : eachindex(Ps) + for s in starts + @test dense_mpo_window(U, s, width) ≈ exp(-im * dt * dense_mpo_window(model, s, width)) atol = 1.0e-10 end - @test log2(errors[1] / errors[2]) ≥ N - 0.15 - end - end -end - - -# Independent random couplings and onsite terms exercise nonuniform finite chains. -function finite_test_hamiltonian(rng, Ps; onsite_only = false, dimers = false) - terms = Pair[] - for s in eachindex(Ps) - a = randn!(rng, zeros(ComplexF64, Ps[s] ← Ps[s])) - push!(terms, s => (a + a') / 2) - if s < length(Ps) && !onsite_only && (!dimers || isodd(s)) - b = randn!(rng, zeros(ComplexF64, Ps[s] ⊗ Ps[s + 1] ← Ps[s] ⊗ Ps[s + 1])) - push!(terms, (s, s + 1) => (b + b') / 2) - end - end - return FiniteMPOHamiltonian(Ps, terms...) -end - -@testset "Finite cluster expansion" begin - rng = MersenneTwister(723) - lattices = ( - [ℂ^2], [ℂ^2, ℂ^3], [ℂ^2, ℂ^3, ℂ^2], [ℂ^2, ℂ^3, ℂ^2, ℂ^3], - [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2], [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2, ℂ^3], - [ℝ^2, ℝ^3, ℝ^2, ℝ^3, ℝ^2], [ℙ^2, ℙ^3, ℙ^2, ℙ^3, ℙ^2], - [Rep[U₁](0 => 1, 1 => 1), dual(Rep[U₁](0 => 1, 1 => 1, 2 => 1)), Rep[U₁](0 => 1, 1 => 1)], - [Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1), Rep[SU₂](1 // 2 => 1), Rep[SU₂](1 => 1), Rep[SU₂](1 // 2 => 1)], - [Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2, 1 => 1), Vect[FermionParity](0 => 1, 1 => 1)], - ) - for Ps in lattices - H = finite_test_hamiltonian(rng, Ps) - L = length(H) - for N in 1:min(6, L) - U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 1.0e-12) - @test U isa FiniteMPO - @test length(U) == L - @test dim(left_virtualspace(U, 1)) == dim(right_virtualspace(U, L)) == 1 - for s in 1:L - @test physicalspace(U, s) == Ps[s] - s < L && @test right_virtualspace(U, s) == left_virtualspace(U, s + 1) - for width in 1:min(N, L - s + 1) - @test dense_mpo_window(U, s, width) ≈ exact_evolution_window(H, -0.1im, s, width) atol = 1.0e-10 - end + model === onsite && @test dim(left_virtualspace(U, 1)) == 1 + end + # The SVD cutoff changes the channels, not the matched exponentials. + for N in (4, 5) + U = @inferred MPSKit.make_cluster_expansion_mpo(H, -0.1im, Val(N), 0.5) + starts = isfinite(H) ? (1:(length(H) - N + 1)) : eachindex(Ps) + for s in starts + @test dense_mpo_window(U, s, N) ≈ exp(-0.1im * dense_mpo_window(H, s, N)) atol = 1.0e-10 end end - # The public entry point caps N at L, including one-site chains. - for N in (L, L + 3), tol in (1.0e-12, 0.5) - U = make_time_mpo(H, 0.1, ClusterExpansion(N; tol)) - @test convert(TensorMap, MPSKit.DenseMPO(U)) ≈ exp(-0.1im * convert(TensorMap, H)) atol = 1.0e-10 - end end end -@testset "Finite edge cases and leading error" begin - rng = MersenneTwister(984) - Ps = [ℂ^2, ℂ^3, ℂ^2, ℂ^3, ℂ^2] - H = finite_test_hamiltonian(rng, Ps) - onsite = finite_test_hamiltonian(rng, Ps; onsite_only = true) - dimers = finite_test_hamiltonian(rng, Ps; dimers = true) - for (model, dt) in ((H, 0.0), (H, -0.1), (H, -0.1im), (onsite, 0.1), (dimers, 0.1)) - U = make_time_mpo(model, dt, ClusterExpansion(5)) - @test convert(TensorMap, MPSKit.DenseMPO(U)) ≈ exp(-im * dt * convert(TensorMap, model)) atol = 1.0e-10 - end - @test dim(left_virtualspace(make_time_mpo(onsite, 0.1, ClusterExpansion(5)), 3)) == 1 - @test_throws ArgumentError make_time_mpo(FiniteMPOHamiltonian(similar(parent(H), 0)), 0.1, ClusterExpansion(2)) - for N in (3, 4) - errors = Float64[] - for dt in (0.1, 0.05) - U = make_time_mpo(H, dt, ClusterExpansion(N)) - push!(errors, norm(convert(TensorMap, MPSKit.DenseMPO(U)) - exp(-im * dt * convert(TensorMap, H)))) - end - @test log2(errors[1] / errors[2]) ≥ N - 0.15 - end - P = ℂ^2 - a = randn!(rng, zeros(ComplexF64, P ← P)) - b = randn!(rng, zeros(ComplexF64, P ⊗ P ← P ⊗ P)) - infinite = InfiniteMPOHamiltonian(PeriodicArray([P]), 1 => (a + a') / 2, (1, 2) => (b + b') / 2) - for L in (4, 6), N in 1:5 - finite = make_time_mpo(open_boundary_conditions(infinite, L), 0.1, ClusterExpansion(N)) - projected = open_boundary_conditions(make_time_mpo(infinite, 0.1, ClusterExpansion(N)), L) - @test convert(TensorMap, MPSKit.DenseMPO(finite)) ≈ convert(TensorMap, MPSKit.DenseMPO(projected)) atol = 1.0e-10 +# Reconstruction alone does not ensure the factors retain all environment directions. +@testset "Completed SVD support" begin + L, R = ℂ^2, ℂ^3 + weak = TensorMap(ComplexF64[1 0 0; 0 1.0e-14 0], L ← R) + for C in (weak, zeros(ComplexF64, L ← R)) + A, B = @inferred MPSKit.factor_with_complementary_spaces(C, 1.0e-12) + @test norm(A * B - C) ≤ 1.0e-16 * norm(C) + @test A * pinv(A; rtol = 1.0e-12) ≈ id(storagetype(C), L) atol = 1.0e-10 + @test pinv(B; rtol = 1.0e-12) * B ≈ id(storagetype(C), R) atol = 1.0e-10 + @test dim(only(domain(A))) == (iszero(norm(C)) ? 5 : 4) end end From 577163cec03495ccb0a44c4685a5c0c3eac31438 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 15:06:54 -0400 Subject: [PATCH 6/8] formatter --- docs/overrides.jl | 6 +++--- test/timeevolution/clusterexpansion.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/overrides.jl b/docs/overrides.jl index 36e0afcce..49e0d3967 100644 --- a/docs/overrides.jl +++ b/docs/overrides.jl @@ -6,7 +6,7 @@ function DocumenterVitepress.vitepress_config_transform(::VitepressOverrides, co occursin(marker, config) || error("DocumenterVitepress config has no default export to extend") return replace(config, marker => "const upstreamConfig = "; count = 1) * """ - import { withOverrides } from './overrides.mts' - export default withOverrides(upstreamConfig) - """ + import { withOverrides } from './overrides.mts' + export default withOverrides(upstreamConfig) + """ end diff --git a/test/timeevolution/clusterexpansion.jl b/test/timeevolution/clusterexpansion.jl index 324044ce4..43fc85756 100644 --- a/test/timeevolution/clusterexpansion.jl +++ b/test/timeevolution/clusterexpansion.jl @@ -57,9 +57,9 @@ end for model in (H, open_boundary_conditions(H, 5)) errors = [ norm( - dense_mpo_window(make_time_mpo(model, dt, ClusterExpansion(N)), 1, N + 1) - + dense_mpo_window(make_time_mpo(model, dt, ClusterExpansion(N)), 1, N + 1) - exp(-im * dt * dense_mpo_window(model, 1, N + 1)) - ) for dt in (0.04, 0.02) + ) for dt in (0.04, 0.02) ] @test log2(errors[1] / errors[2]) ≥ N - 0.15 end From c4ca7a220b76177232943242413662c001fa2c84 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 08:40:19 -0400 Subject: [PATCH 7/8] Simplify cluster expansion helpers and fix Julia 1.10 inference --- src/algorithms/timestep/clusterexpansion.jl | 288 ++++++++++---------- 1 file changed, 143 insertions(+), 145 deletions(-) diff --git a/src/algorithms/timestep/clusterexpansion.jl b/src/algorithms/timestep/clusterexpansion.jl index c17f34636..44d38b083 100644 --- a/src/algorithms/timestep/clusterexpansion.jl +++ b/src/algorithms/timestep/clusterexpansion.jl @@ -64,177 +64,178 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::ClusterExpansion) return make_cluster_expansion_mpo(H, -im * dt, Val(N), alg.tol) end -function make_cluster_expansion_mpo(H::MPOHamiltonian, τ::Number, n::Val, tol::Real) - T = promote_type(scalartype(H), typeof(τ)) - storage = TensorKit.similarstoragetype(storagetype(H), T) - O = map(parent(H)) do h - return exponentiate_onsite_tensor(h, τ, storage) +function make_cluster_expansion_mpo(H::MPOHamiltonian, τ::Number, ::Val{N}, tol::Real) where {N} + storage_type = TensorKit.similarstoragetype(storagetype(H), promote_type(scalartype(H), typeof(τ))) + return make_cluster_expansion_mpo(H, τ, Val(N), tol, storage_type) +end +function make_cluster_expansion_mpo( + H::MPOHamiltonian, τ::Number, ::Val{N}, tol::Real, ::Type{A} + ) where {N, A} + tensors = map(parent(H)) do ham_tensor + return exponentiate_onsite_tensor(ham_tensor, τ, A) end # These cases are already exact products; avoid building channels from # roundoff in differences of independently contracted exponentials. - (iszero(τ) || all(has_only_onsite_terms, parent(H))) && return assemble_evolution_mpo(O) - return assemble_evolution_mpo(add_cluster_corrections(O, H, τ, n, tol)) + if !(iszero(τ) || all(has_only_onsite_terms, parent(H))) + tensors = add_cluster_corrections(tensors, H, τ, Val(N), tol) + end + + return assemble_evolution_mpo(tensors) end -has_only_onsite_terms(h) = size(h, 1) ≤ 2 && size(h, 4) ≤ 2 +has_only_onsite_terms(ham_tensor) = size(ham_tensor, 1) ≤ 2 && size(ham_tensor, 4) ≤ 2 -assemble_evolution_mpo(O::PeriodicVector) = InfiniteMPO(O) -function assemble_evolution_mpo(O::Vector) - O[1] = O[1][1, :, :, :] - O[end] = O[end][:, :, :, 1] +assemble_evolution_mpo(tensors::PeriodicVector) = InfiniteMPO(tensors) +function assemble_evolution_mpo(tensors::Vector) + tensors[1] = tensors[1][1, :, :, :] + tensors[end] = tensors[end][:, :, :, 1] # Remove structurally unused levels only after all stages have finished. - return remove_orphans!(FiniteMPO(O); tol = 0) + return remove_orphans!(FiniteMPO(tensors); tol = 0) end -function exponentiate_onsite_tensor(h, τ, ::Type{A}) where {A} - P = physicalspace(h) - levels = ⊞(oneunit(P)) - blocktype = tensormaptype(spacetype(P), 2, 2, A) - O = SparseBlockTensorMap{blocktype}(undef, levels ⊗ P ← P ⊗ levels) - onsite = removeunit(removeunit(h[1, 1, 1, end], 4), 1) - O[1, 1, 1, 1] = add_util_leg(exp(τ * onsite)) - return O +function exponentiate_onsite_tensor(ham_tensor, τ, ::Type{A}) where {A} + physical_space = physicalspace(ham_tensor) + levels = ⊞(oneunit(physical_space)) + blocktype = tensormaptype(spacetype(physical_space), 2, 2, A) + tensor = SparseBlockTensorMap{blocktype}(undef, levels ⊗ physical_space ← physical_space ⊗ levels) + onsite = removeunit(removeunit(ham_tensor[1, 1, 1, end], 4), 1) + tensor[1, 1, 1, 1] = add_util_leg(exp(τ * onsite)) + return tensor end -add_cluster_corrections(O, H, τ, ::Val{1}, tol) = O -function add_cluster_corrections(O, H, τ, n::Val{N}, tol) where {N} - previous = add_cluster_corrections(O, H, τ, Val(N - 1), tol) +function add_cluster_corrections(tensors, H, τ, ::Val{N}, tol) where {N} + N == 1 && return tensors + previous = add_cluster_corrections(tensors, H, τ, Val(N - 1), tol) starts = isfinite(H) ? (1:(length(H) - N + 1)) : (1:length(H)) - nt = (N - 1) ÷ 2 + environment_length = (N - 1) ÷ 2 level = N ÷ 2 + 1 # Every translated window must see the completed previous stage. Compute # all corrections before inserting any of them, including across the cell seam. if isodd(N) - centers = [solve_center_correction(previous, s, evolution_cluster_residual(H, previous, s, τ, n), n) for s in starts] - for s in starts - previous[s + nt][level, 1, 1, level] = centers[s] + centers = [solve_center_correction(previous, start, evolution_cluster_residual(H, previous, start, τ, Val(N)), Val(N)) for start in starts] + for start in starts + previous[start + environment_length][level, 1, 1, level] = centers[start] end return previous else - pairs = [factor_two_site_correction(H, previous, s, τ, n, tol) for s in starts] - expanded = expand_correction_bonds(previous, pairs, n) - for s in starts - expanded[s + nt][level - 1, 1, 1, level] = pairs[s][1] - expanded[s + nt + 1][level, 1, 1, level - 1] = pairs[s][2] + factors = [factor_two_site_correction(H, previous, start, τ, Val(N), tol) for start in starts] + expanded = expand_correction_bonds(previous, factors, Val(N)) + for start in starts + expanded[start + environment_length][level - 1, 1, 1, level] = factors[start][1] + expanded[start + environment_length + 1][level, 1, 1, level - 1] = factors[start][2] end return expanded end end -function expand_correction_bonds(O::PeriodicVector, pairs, ::Val{N}) where {N} - nt = (N - 1) ÷ 2 - bonds = PeriodicArray([right_virtualspace(pairs[mod1(i - nt, length(O))][1]) for i in 1:length(O)]) - return PeriodicArray([append_virtual_spaces(O[i], bonds[i - 1], bonds[i]) for i in 1:length(O)]) +function expand_correction_bonds(tensors::PeriodicVector, factors, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + bonds = PeriodicArray([right_virtualspace(factors[mod1(i - environment_length, length(tensors))][1]) for i in 1:length(tensors)]) + return PeriodicArray([append_virtual_spaces(tensors[i], bonds[i - 1], bonds[i]) for i in 1:length(tensors)]) end -function expand_correction_bonds(O::Vector, pairs, ::Val{N}) where {N} - nt = (N - 1) ÷ 2 - empty_bond = zero(left_virtualspace(O[1])[1]) +function expand_correction_bonds(tensors::Vector, factors, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + empty_bond = zero(left_virtualspace(tensors[1])[1]) # Only the middle bonds can support an N-site cluster. Empty levels keep # the same level indices everywhere until the final boundary projection. bonds = vcat( - fill(empty_bond, nt + 1), - [right_virtualspace(pair[1]) for pair in pairs], - fill(empty_bond, nt + 1) + fill(empty_bond, environment_length + 1), + [right_virtualspace(factor[1]) for factor in factors], + fill(empty_bond, environment_length + 1) ) - return [append_virtual_spaces(O[i], bonds[i], bonds[i + 1]) for i in 1:length(O)] + return [append_virtual_spaces(tensors[i], bonds[i], bonds[i + 1]) for i in 1:length(tensors)] end -function append_virtual_spaces(O, left_bond, right_bond) - L, R = left_virtualspace(O), right_virtualspace(O) - left_levels = L ⊞ left_bond - right_levels = R ⊞ right_bond - P = physicalspace(O) - expanded = typeof(O)(undef, left_levels ⊗ P ← P ⊗ right_levels) - for (indices, block) in nonzero_pairs(O) +function append_virtual_spaces(tensor, left_bond, right_bond) + left_levels = left_virtualspace(tensor) ⊞ left_bond + right_levels = right_virtualspace(tensor) ⊞ right_bond + physical_space = physicalspace(tensor) + expanded = typeof(tensor)(undef, left_levels ⊗ physical_space ← physical_space ⊗ right_levels) + for (indices, block) in nonzero_pairs(tensor) expanded[indices] = block end return expanded end -function evolution_cluster_residual(H, O, start::Int, τ, n::Val{N}) where {N} +function evolution_cluster_residual(H, tensors, start::Int, τ, ::Val{N}) where {N} boundary = size(H[start + N - 1], 4) - exact = exp(τ * contract_mpo_window(H, start, n, boundary)) - return add_util_leg(exact - contract_mpo_window(O, start, n, 1)) -end - -function factor_two_site_correction(H, O, start, τ, n::Val{2}, tol) - residual = evolution_cluster_residual(H, O, start, τ, n) - A, B = factor_with_complementary_spaces(permute(residual, ((1, 2, 4), (3, 5, 6))), tol) - return permute(A, ((1, 2), (3, 4))), permute(B, ((1, 2), (3, 4))) + exact = exp(τ * contract_mpo_window(H, start, Val(N), boundary)) + return add_util_leg(exact - contract_mpo_window(tensors, start, Val(N), 1)) end -function factor_two_site_correction(H, O, start, τ, n::Val{N}, tol) where {N} - residual = evolution_cluster_residual(H, O, start, τ, n) - left_map, right_map = correction_environment_maps(O, start, n) - center = solve_center_correction(left_map, right_map, residual, n) - C = permute(center, ((1, 2, 4), (3, 5, 6))) +function factor_two_site_correction(H, tensors, start, τ, ::Val{N}, tol) where {N} + residual = evolution_cluster_residual(H, tensors, start, τ, Val(N)) + if N == 2 + left_factor, right_factor = factor_with_complementary_spaces(permute(residual, ((1, 2, 4), (3, 5, 6))), tol) + return permute(left_factor, ((1, 2), (3, 4))), permute(right_factor, ((1, 2), (3, 4))) + end + left_map, right_map = correction_environment_maps(tensors, start, Val(N)) + center = solve_center_correction(left_map, right_map, residual, Val(N)) + center_map = permute(center, ((1, 2, 4), (3, 5, 6))) # The minimum-norm solves restrict the center to the active environment # subspaces. Complete only these directions, rather than the redundant # virtual directions introduced by previous stages. Untruncated LQ/QR bases # retain the entire physical environment support, including weak sectors. - _, Q_left = right_orth(left_map) - Q_right, _ = left_orth(right_map) - F_left = Q_left' ⊗ id(storagetype(C), codomain(C)[2] ⊗ codomain(C)[3]) - F_right = id(storagetype(C), domain(C)[1] ⊗ domain(C)[2]) ⊗ Q_right - A, B = factor_with_complementary_spaces(F_left' * C * F_right, tol) - return permute(F_left * A, ((1, 2), (3, 4))), - permute(B * F_right', ((1, 2), (3, 4))) + _, left_basis = right_orth(left_map) + right_basis, _ = left_orth(right_map) + left_embedding = left_basis' ⊗ id(storagetype(center_map), codomain(center_map)[2] ⊗ codomain(center_map)[3]) + right_embedding = id(storagetype(center_map), domain(center_map)[1] ⊗ domain(center_map)[2]) ⊗ right_basis + left_factor, right_factor = factor_with_complementary_spaces(left_embedding' * center_map * right_embedding, tol) + return permute(left_embedding * left_factor, ((1, 2), (3, 4))), + permute(right_factor * right_embedding', ((1, 2), (3, 4))) end """ - factor_with_complementary_spaces(C, tol) + factor_with_complementary_spaces(tensor, tol) -Factor `C` through its shared SVD support and separate left/right complementary +Factor `tensor` through its shared SVD support and separate left/right complementary spaces. For exact rank `r_c`, the virtual multiplicity in sector `c` is `m_c + n_c - r_c`, the minimum that gives both factors full environment support. For exact null directions, the factors have the form A = [U√S γU⊥ 0], B = [√S Vᴴ; 0; γV⊥ᴴ]. -At a finite rank cutoff, retain the remaining block `E = U⊥' C V⊥` in the -complementary channels, so the product still reconstructs `C` up to roundoff. +At a finite rank cutoff, retain the remaining block in the complementary +channels, so the product still reconstructs `tensor` up to roundoff. The scale `γ` balances these channels against the retained singular values. """ -function factor_with_complementary_spaces(C::TensorMap, tol::Real) - cnorm = norm(C) - if iszero(cnorm) +function factor_with_complementary_spaces(tensor::TensorMap, tol::Real) + tensor_norm = norm(tensor) + if iszero(tensor_norm) # At zero rank, the two complementary spaces occupy disjoint channels. - L, R = fuse(codomain(C)), fuse(domain(C)) - storage = storagetype(C) - A = catdomain( - isomorphism(storage, codomain(C) ← L), zeros(storage, codomain(C) ← R) + left_space, right_space = fuse(codomain(tensor)), fuse(domain(tensor)) + storage_type = storagetype(tensor) + left_factor = catdomain( + isomorphism(storage_type, codomain(tensor) ← left_space), zeros(storage_type, codomain(tensor) ← right_space) ) - B = catcodomain( - zeros(storage, L ← domain(C)), isomorphism(storage, R ← domain(C)) + right_factor = catcodomain( + zeros(storage_type, left_space ← domain(tensor)), isomorphism(storage_type, right_space ← domain(tensor)) ) - return A, B + return left_factor, right_factor end - U, S, Vᴴ, _ = svd_trunc(C; trunc = trunctol(; rtol = tol, p = Inf)) + U, S, Vᴴ, _ = svd_trunc(tensor; trunc = trunctol(; rtol = tol, p = Inf)) U_perp, V_perp = left_null(U), right_null(Vᴴ) - E = U_perp' * C * V_perp' - γ = sqrt(cnorm) - root = sqrt(S) - A = catdomain(catdomain(U * root, γ * U_perp), U_perp * E / (2γ)) - B = catcodomain(catcodomain(root * Vᴴ, E * V_perp / (2γ)), γ * V_perp) - return A, B -end - -function contract_mpo_window(O, start::Int, ::Val{1}, boundary::Int) - return removeunit(removeunit(TensorMap(O[start][1, :, :, boundary]), 4), 1) + complement = U_perp' * tensor * V_perp' + complement_scale = sqrt(tensor_norm) + sqrt_values = sqrt(S) + left_factor = catdomain(catdomain(U * sqrt_values, complement_scale * U_perp), U_perp * complement / (2complement_scale)) + right_factor = catcodomain(catcodomain(sqrt_values * Vᴴ, complement * V_perp / (2complement_scale)), complement_scale * V_perp) + return left_factor, right_factor end -function contract_mpo_window(O, start::Int, n::Val{N}, boundary::Int) where {N} - left = removeunit(TensorMap(O[start][1, :, :, :]), 1) - right = removeunit(TensorMap(O[start + N - 1][:, :, :, boundary]), 4) - sites = collect_bulk_tensors(O, start + 1, Val(N - 2)) - return contract_mpo_window(left, sites, right, n) +function contract_mpo_window(tensors, start::Int, ::Val{N}, boundary::Int) where {N} + N == 1 && return removeunit(removeunit(TensorMap(tensors[start][1, :, :, boundary]), 4), 1) + left = removeunit(TensorMap(tensors[start][1, :, :, :]), 1) + right = removeunit(TensorMap(tensors[start + N - 1][:, :, :, boundary]), 4) + sites = collect_bulk_tensors(tensors, start + 1, Val(N - 2)) + return contract_mpo_window(left, sites, right, Val(N)) end -collect_bulk_tensors(O, start, ::Val{0}) = () -function collect_bulk_tensors(O, start, ::Val{N}) where {N} - return (TensorMap(O[start]), collect_bulk_tensors(O, start + 1, Val(N - 1))...) +function collect_bulk_tensors(tensors, start, ::Val{N}) where {N} + N == 0 && return () + return (TensorMap(tensors[start]), collect_bulk_tensors(tensors, start + 1, Val(N - 1))...) end """ @@ -259,23 +260,23 @@ are above the tensors and codomain legs below. left::AbstractTensorMap{<:Any, S, 1, 2}, sites::NTuple{K, AbstractTensorMap{<:Any, S, 2, 2}}, right::AbstractTensorMap{<:Any, S, 2, 1}, ::Val{N} ) where {S, N, K} - out = tensorexpr(:cluster, -(1:N), -((N + 1):(2N))) + result = tensorexpr(:cluster, -(1:N), -((N + 1):(2N))) first_site = tensorexpr(:left, -1, (-N - 1, 1)) last_site = tensorexpr(:right, (N - 1, -N), -2N) - middle = [tensorexpr(:(sites[$(i - 1)]), (i - 1, -i), (-N - i, i)) for i in 2:(N - 1)] - return macroexpand(@__MODULE__, :(return @plansor $out := *($first_site, $(middle...), $last_site))) + bulk = [tensorexpr(:(sites[$(i - 1)]), (i - 1, -i), (-N - i, i)) for i in 2:(N - 1)] + return macroexpand(@__MODULE__, :(return @plansor $result := *($first_site, $(bulk...), $last_site))) end -left_correction_environment(O, start, ::Val{1}) = O[start][1, 1, 1, 2] -function left_correction_environment(O, start, ::Val{L}) where {L} - left = left_correction_environment(O, start, Val(L - 1)) - return extend_left_environment(left, O[start + L - 1][L, 1, 1, L + 1]) +function left_correction_environment(tensors, start, ::Val{N}) where {N} + N == 1 && return tensors[start][1, 1, 1, 2] + left = left_correction_environment(tensors, start, Val(N - 1)) + return extend_left_environment(left, tensors[start + N - 1][N, 1, 1, N + 1]) end -right_correction_environment(O, start, ::Val{1}) = O[start][2, 1, 1, 1] -function right_correction_environment(O, start, ::Val{L}) where {L} - right = right_correction_environment(O, start + 1, Val(L - 1)) - return extend_right_environment(O[start][L + 1, 1, 1, L], right) +function right_correction_environment(tensors, start, ::Val{N}) where {N} + N == 1 && return tensors[start][2, 1, 1, 1] + right = right_correction_environment(tensors, start + 1, Val(N - 1)) + return extend_right_environment(tensors[start][N + 1, 1, 1, N], right) end """ @@ -297,10 +298,10 @@ domain legs are above the tensors and codomain legs below. ``` """ @generated function extend_left_environment(left::AbstractTensorMap{<:Any, <:Any, K, K}, site::AbstractTensorMap) where {K} - out = tensorexpr(:next_left, -(1:(K + 1)), -((K + 2):(2K + 2))) - env = tensorexpr(:left, -(1:K), (-((K + 2):(2K))..., 1)) - op = tensorexpr(:site, (1, -K - 1), (-2K - 1, -2K - 2)) - return macroexpand(@__MODULE__, :(return @plansor $out := $env * $op)) + result = tensorexpr(:next_left, -(1:(K + 1)), -((K + 2):(2K + 2))) + environment = tensorexpr(:left, -(1:K), (-((K + 2):(2K))..., 1)) + site_expression = tensorexpr(:site, (1, -K - 1), (-2K - 1, -2K - 2)) + return macroexpand(@__MODULE__, :(return @plansor $result := $environment * $site_expression)) end """ @@ -322,24 +323,23 @@ domain legs are above the tensors and codomain legs below. ``` """ @generated function extend_right_environment(site::AbstractTensorMap, right::AbstractTensorMap{<:Any, <:Any, K, K}) where {K} - out = tensorexpr(:next_right, -(1:(K + 1)), -((K + 2):(2K + 2))) - op = tensorexpr(:site, (-1, -2), (-K - 2, 1)) - env = tensorexpr(:right, (1, -(3:(K + 1))...), -((K + 3):(2K + 2))) - return macroexpand(@__MODULE__, :(return @plansor $out := $op * $env)) + result = tensorexpr(:next_right, -(1:(K + 1)), -((K + 2):(2K + 2))) + site_expression = tensorexpr(:site, (-1, -2), (-K - 2, 1)) + environment = tensorexpr(:right, (1, -(3:(K + 1))...), -((K + 3):(2K + 2))) + return macroexpand(@__MODULE__, :(return @plansor $result := $site_expression * $environment)) end -function correction_environment_maps(O, start, ::Val{N}) where {N} - L = (N - 1) ÷ 2 - left = left_correction_environment(O, start, Val(L)) - right = right_correction_environment(O, start + N - L, Val(L)) - indices = ntuple(identity, Val(2L + 1)) - return permute(left, (indices, (2L + 2,))), - permute(right, ((1,), indices .+ 1)) +function correction_environment_maps(tensors, start, ::Val{N}) where {N} + environment_length = (N - 1) ÷ 2 + left = left_correction_environment(tensors, start, Val(environment_length)) + right = right_correction_environment(tensors, start + N - environment_length, Val(environment_length)) + environment_indices = ntuple(identity, Val(2environment_length + 1)) + return permute(left, (environment_indices, (2environment_length + 2,))), permute(right, ((1,), environment_indices .+ 1)) end -function solve_center_correction(O, start::Int, residual::TensorMap, n::Val{N}) where {N} - left_map, right_map = correction_environment_maps(O, start, n) - return solve_center_correction(left_map, right_map, residual, n) +function solve_center_correction(tensors, start::Int, residual::TensorMap, ::Val{N}) where {N} + left_map, right_map = correction_environment_maps(tensors, start, Val(N)) + return solve_center_correction(left_map, right_map, residual, Val(N)) end """ @@ -384,20 +384,18 @@ clusters, each outer physical space above is a bundle of `(N - 1) ÷ 2` legs. @generated function solve_center_correction( left_map::TensorMap, right_map::TensorMap, residual::TensorMap, ::Val{N} ) where {N} - nt = (N - 1) ÷ 2 - width = N - 2nt + environment_length = (N - 1) ÷ 2 + center_width = N - 2environment_length # Group all left environment legs into the codomain for the left solve. - left_indices = (1, (2:(nt + 1))..., ((N + 2):(N + nt + 1))...) - remaining = (((nt + 2):(N + 1))..., ((N + nt + 2):(2N + 2))...) + left_indices = (1, (2:(environment_length + 1))..., ((N + 2):(N + environment_length + 1))...) + remaining_indices = (((environment_length + 2):(N + 1))..., ((N + environment_length + 2):(2N + 2))...) # After solving, group the right environment legs into the domain. - middle_indices = (1, (2:(width + 1))..., ((nt + width + 2):(nt + 2width + 1))...) - right_indices = (((width + 2):(nt + width + 1))..., ((nt + 2width + 2):(2nt + 2width + 2))...) - center_indices = (Tuple(1:(width + 1)), Tuple((width + 2):(2width + 2))) + middle_indices = (1, (2:(center_width + 1))..., ((environment_length + center_width + 2):(environment_length + 2center_width + 1))...) + right_indices = (((center_width + 2):(environment_length + center_width + 1))..., ((environment_length + 2center_width + 2):(2environment_length + 2center_width + 2))...) + center_indices = (Tuple(1:(center_width + 1)), Tuple((center_width + 2):(2center_width + 2))) return quote - left_rhs = permute(residual, $((left_indices, remaining))) - left_solved = left_map \ left_rhs - right_rhs = permute(left_solved, $((middle_indices, right_indices))) - solved = right_rhs / right_map + left_solved = left_map \ permute(residual, $((left_indices, remaining_indices))) + solved = permute(left_solved, $((middle_indices, right_indices))) / right_map return permute(solved, $center_indices) end end From 608c08d828213ceb7a231ba8684d54d322657609 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 11:56:01 -0400 Subject: [PATCH 8/8] update changelog --- docs/src/changelog.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index c0fcf73b6..7efbec564 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -21,6 +21,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs + for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion + that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) + ### Changed ### Deprecated @@ -29,6 +33,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Fixed +- `convert(TensorMap, H)` for single-site `FiniteMPOHamiltonian`s no longer contracts the site tensor with itself. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) + ### Performance ## [0.14.0](https://github.com/QuantumKitHub/MPSKit.jl/compare/v0.13.13...v0.14.0) - 2026-10-05