From 63d0152c96333b2c03613203fdb206730be9f5a2 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 29 May 2026 11:16:26 -0400 Subject: [PATCH 01/20] initial CBE implementation --- src/MPSKit.jl | 2 +- src/algorithms/groundstate/dmrg.jl | 151 ++++++++++++++++++++++++++++- 2 files changed, 150 insertions(+), 3 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index a4e1fe605..03d54c215 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -31,7 +31,7 @@ export leftenv, rightenv export find_groundstate, find_groundstate! export leading_boundary export approximate, approximate! -export VUMPS, VOMPS, DMRG, DMRG2, IDMRG, IDMRG2, GradientGrassmann +export VUMPS, VOMPS, DMRG, DMRG2, CBEDMRG, IDMRG, IDMRG2, GradientGrassmann export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index ca8b09d18..a4f177b88 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -78,6 +78,154 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme return ψ, envs, ϵ end +struct CBEDMRG{A, F} <: Algorithm + "tolerance for convergence criterium" + tol::Float64 + + "maximal amount of iterations" + maxiter::Int + + "setting for how much information is displayed" + verbosity::Int + + "algorithm used for the eigenvalue solvers" + alg_eigsolve::A + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F + + selection::TruncationStrategy + + "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update" + trscheme::TruncationStrategy +end + +function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg::CBEDMRG, envs = environments(ψ, H)) + ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) + ϵ = maximum(ϵs) + # extra debug error measures (do not drive convergence; ϵ remains the two-site fidelity) + ϵs_galerkin = zero(ϵs) # local (one-site) Galerkin error + ϵs_trunc = zero(ϵs) # bond-truncation error (discarded weight) + ϵs_2site = zero(ϵs) # two-site Galerkin error (complement-projected two-site gradient) + log = IterLog("DMRG") + + LoggingExtras.withlevel(; alg.verbosity) do + @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs)) + for iter in 1:(alg.maxiter) + alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) + + zerovector!(ϵs) + zerovector!(ϵs_galerkin) + zerovector!(ϵs_trunc) + zerovector!(ϵs_2site) + for pos in 1:(length(ψ) - 1) + @plansor ac2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3] + + # exact two-site update at bond (pos, pos+1), projected onto the complement + # of the current state (this is the exact analog of the randomized SVD) + AC2 = AC2_projection(pos, ψ, H, ψ, envs) + NL = left_null(ψ.AC[pos]) + NR = right_null!(_transpose_tail(ψ.AR[pos + 1]; copy = true)) + g2 = adjoint(NL) * AC2 * adjoint(NR) + # two-site Galerkin: norm of the two-site gradient in the complement, normalized + ϵs_2site[pos] = max(ϵs_2site[pos], norm(g2) / norm(AC2)) + intermediate = normalize!(g2) + _, _, V = svd_trunc!(intermediate; trunc = alg.selection, alg = Defaults.alg_svd()) + + # expand the bond: optimal vectors at pos+1, zero weight at pos + ar_re = V * NR + ar_le = zerovector!(similar(ψ.AC[pos], codomain(ψ.AC[pos]) ← space(V, 1))) + Ql, C = qr_compact!(catdomain(ψ.AC[pos], ar_le)) + AR_exp = _transpose_front(catcodomain(_transpose_tail(ψ.AR[pos + 1]), ar_re)) + ψ.AC[pos] = (Ql, normalize!(C)) + ψ.AC[pos + 1] = (C, AR_exp) + + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular + # values, computed directly from the spectrum (no 1 - ‖C‖² cancellation) + U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) + ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) + AL′ = U + C = S * Vᴴ + normalize!(C) + + # DMRG2-style error: 1 - fidelity between the two-site state before the + # eigensolver and the truncated, optimized two-site state after it. The + # truncated bond is internal, so this -> 0 at convergence. + @plansor ac2′[-1 -2; -3 -4] := AL′[-1 -2; 1] * C[1; 2] * AR_exp[2 -4; -3] + ϵs[pos] = max(ϵs[pos], abs(1 - abs(dot(ac2, ac2′)))) + + ψ.AC[pos] = (AL′, C) + ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) + # @debug "CBEDMRG L→R" pos ϵ = ϵs[pos] + end + + for pos in length(ψ):-1:2 + @plansor ac2[-1 -2; -3 -4] := ψ.AL[pos - 1][-1 -2; 1] * ψ.AC[pos][1 -4; -3] + + # exact two-site update at bond (pos-1, pos), projected onto the complement + AC2 = AC2_projection(pos - 1, ψ, H, ψ, envs) + NL = left_null(ψ.AL[pos - 1]) + NR = right_null!(_transpose_tail(ψ.AC[pos]; copy = true)) + g2 = adjoint(NL) * AC2 * adjoint(NR) + # two-site Galerkin: norm of the two-site gradient in the complement, normalized + ϵs_2site[pos - 1] = max(ϵs_2site[pos - 1], norm(g2) / norm(AC2)) + intermediate = normalize!(g2) + U, _, _ = svd_trunc!(intermediate; trunc = alg.selection, alg = Defaults.alg_svd()) + + # optimal new left vectors for AL[pos-1]; zero padding for AC[pos]'s left virtual + Q = NL * U + AL = ψ.AL[pos - 1] + al_le = zerovector!( + similar(ψ.AC[pos], space(Q, 3)' ⊗ physicalspace(ψ.AC[pos]) ← right_virtualspace(ψ.AC[pos])) + ) + C, Qr = lq_compact!(catcodomain(_transpose_tail(ψ.AC[pos]), _transpose_tail(al_le))) + AL_exp = catdomain(AL, Q) + ψ.AC[pos] = (normalize!(C), _transpose_front(Qr)) + ψ.AC[pos - 1] = (AL_exp, C) + + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular + # values, computed directly from the spectrum (no 1 - ‖C‖² cancellation) + U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) + ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) + C = U * S + normalize!(C) + AR_mps = _transpose_front(AR′) + + # DMRG2-style error: 1 - fidelity between the two-site state before the + # eigensolver and the truncated, optimized two-site state after it. The + # truncated bond is internal, so this -> 0 at convergence. AR_mps is in MPS + # form (physical in codomain) so the overlap with ac2 stays planar. + @plansor ac2′[-1 -2; -3 -4] := AL_exp[-1 -2; 1] * C[1; 2] * AR_mps[2 -4; -3] + ϵs[pos] = max(ϵs[pos], abs(1 - abs(dot(ac2, ac2′)))) + # @debug "CBEDMRG R→L" pos ϵ = ϵs[pos] + ψ.AC[pos] = (C, AR_mps) + ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) + end + + ϵ = maximum(ϵs) + @infov 2 "CBEDMRG errors" iter fidelity = ϵ local_galerkin = maximum(ϵs_galerkin) twosite_galerkin = maximum(ϵs_2site) truncation = maximum(ϵs_trunc) + + ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} + + if ϵ <= alg.tol + @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs)) + break + end + if iter == alg.maxiter + @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs)) + else + @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs)) + end + end + end + return ψ, envs, ϵ +end + + """ $(TYPEDEF) @@ -172,7 +320,6 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm ψ.AC[pos] = (al, complex(c)) end end - end ϵ = maximum(ϵs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( @@ -195,6 +342,6 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm return ψ, envs, ϵ end -function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) +function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2, CBEDMRG}, envs...; kwargs...) return find_groundstate!(copy(ψ), H, alg, envs...; kwargs...) end From b37d24e50124c222c3abd9eb33b5b0146aa87f88 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Jun 2026 11:23:33 -0400 Subject: [PATCH 02/20] add `changebond!` implementation --- src/algorithms/changebonds/changebonds.jl | 34 ++++++-- src/algorithms/changebonds/optimalexpand.jl | 88 ++++++++++++++------- src/algorithms/changebonds/randexpand.jl | 75 +++++++++++++----- 3 files changed, 143 insertions(+), 54 deletions(-) diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl index de549ca24..56de20f1c 100644 --- a/src/algorithms/changebonds/changebonds.jl +++ b/src/algorithms/changebonds/changebonds.jl @@ -12,27 +12,41 @@ See also: [`SvdCut`](@ref), [`RandExpand`](@ref), [`VUMPSSvdCut`](@ref), [`Optim function changebonds end function changebonds! end +@doc """ + changebond(site, dir, ψ, [H], alg, [envs]) -> ψ, info + changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ, info + +Expand a single bond of `ψ` in place by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the enriched bond. +The sweep direction `dir` is a `Val(:right)` or `Val(:left)` used for dispatch. +For `Val(:right)` the bond `(site, site + 1)` is enriched on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions; +for `Val(:left)` the mirror is applied to bond `(site - 1, site)`. + +See also [`changebonds`](@ref), [`changebonds!`](@ref). +""" changebond, changebond! +function changebond end +function changebond! end + _expand(ψ, AL′, AR′) = _expand!(copy(ψ), AL′, AR′) function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector) for i in 1:length(ψ) # update AL: add vectors, make room for new vectors: # AL -> [AL expansion; 0 0] al′ = _transpose_tail(catdomain(ψ.AL[i], AL′[i])) - lz = zerovector!(similar(al′, _lastspace(AL′[i - 1])' ← domain(al′))) - ψ.AL[i] = _transpose_front(catcodomain(al′, lz)) + al_space = (codomain(al′)[1] ⊕ _lastspace(AL′[i - 1])') ← domain(al′) + ψ.AL[i] = _transpose_front(absorb!(zerovector!(similar(al′, al_space)), al′)) # update AR: add vectors, make room for new vectors: # AR -> [AR 0; expansion 0] ar′ = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), AR′[i + 1])) - rz = zerovector!(similar(ar′, codomain(ar′) ← _firstspace(AR′[i + 2]))) - ψ.AR[i + 1] = catdomain(ar′, rz) + ar_space = codomain(ar′) ← (domain(ar′)[1] ⊕ _firstspace(AR′[i + 2])) + ψ.AR[i + 1] = absorb!(zerovector!(similar(ar′, ar_space)), ar′) # update C: add vectors, make room for new vectors: # C -> [C 0; 0 expansion] - l = zerovector!(similar(ψ.C[i], codomain(ψ.C[i]) ← _firstspace(AR′[i + 1]))) - ψ.C[i] = catdomain(ψ.C[i], l) - r = zerovector!(similar(ψ.C[i], _lastspace(AL′[i])' ← domain(ψ.C[i]))) - ψ.C[i] = catcodomain(ψ.C[i], r) + c_dom = codomain(ψ.C[i]) ← (domain(ψ.C[i])[1] ⊕ _firstspace(AR′[i + 1])) + ψ.C[i] = absorb!(zerovector!(similar(ψ.C[i], c_dom)), ψ.C[i]) + c_cod = (codomain(ψ.C[i])[1] ⊕ _lastspace(AL′[i])') ← domain(ψ.C[i]) + ψ.C[i] = absorb!(zerovector!(similar(ψ.C[i], c_cod)), ψ.C[i]) # update AC: recalculate ψ.AC[i] = ψ.AL[i] * ψ.C[i] @@ -45,3 +59,7 @@ function _expand!(ψ::MultilineMPS, AL′::PeriodicMatrix, AR′::PeriodicMatrix end return ψ end + +function changebond(site::Int, dir::Val, ψ::AbstractFiniteMPS, H, alg::Algorithm, envs) + return changebond!(site, dir, copy(ψ), H, alg, envs) +end diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 194808704..28b46f467 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -81,37 +81,69 @@ function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environment return newψ, envs end -function changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) - return changebonds!(copy(ψ), H, alg, envs) -end -function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) - #inspired by the infinite mps algorithm, alternative is to use https://arxiv.org/pdf/1501.05504.pdf +# Finite system +# ------------- +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs) + bond = site + left = ψ.AC[site] + right = ψ.AR[site + 1] + NL = left_null(left) + NR = right_null!(_transpose_tail(right; copy = true)) + + # two-site update from the projected effective Hamiltonian + AC2 = AC2_projection(bond, ψ, H, ψ, envs) + + # select the dominant directions in the complement of the current state + g2 = adjoint(NL) * AC2 * adjoint(NR) + ϵ_2site = norm(g2) / norm(AC2) + _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + + # optimal vectors at site+1, zero weight at site + ar_re = Vᴴ * NR + # embed `left` into the enlarged domain (zero weight in the new directions) + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) + + ψ.AC[site] = (nal, normalize!(nc)) + ψ.AC[site + 1] = (nc, nar) + return ψ, (; ϵ_select, ϵ_2site) +end +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs) + bond = site - 1 + left = ψ.AL[site - 1] + right = ψ.AC[site] + NL = left_null(left) + NR = right_null!(_transpose_tail(right; copy = true)) + + # two-site update from the projected effective Hamiltonian + AC2 = AC2_projection(bond, ψ, H, ψ, envs) + + # select the dominant directions in the complement of the current state + g2 = adjoint(NL) * AC2 * adjoint(NR) + ϵ_2site = norm(g2) / norm(AC2) + U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + + # optimal vectors at site-1, zero weight at site + Q = NL * U + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) + right_tail = _transpose_tail(right) + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) + AL_exp = catdomain(left, Q) + + ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) + ψ.AC[site - 1] = (AL_exp, nc) + return ψ, (; ϵ_select, ϵ_2site) +end - #the idea is that we always want to expand the state in such a way that there are zeros at site i - #but "optimal vectors" at site i+1 - #so during optimization of site i, you have access to these optimal vectors :) +changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) = + changebonds!(copy(ψ), H, alg, envs) +function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) for i in 1:(length(ψ) - 1) - AC2 = AC2_projection(i, ψ, H, ψ, envs) - - #Calculate nullspaces for left and right - NL = left_null(ψ.AC[i]) - NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) - - #Use this nullspaces and SVD decomposition to determine the optimal expansion space - intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR)) - _, _, V, = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd) - - ar_re = V * NR - ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(V, 1))) - - nal, nc = qr_compact!(catdomain(ψ.AC[i], ar_le)) - nar = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), ar_re)) - - ψ.AC[i] = (nal, nc) - ψ.AC[i + 1] = (nc, nar) + changebond!(i, Val(:right), ψ, H, alg, envs) end - - return (ψ, envs) + return ψ, envs end diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 346755e4c..af849c28c 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -68,30 +68,69 @@ end function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) + # the expansion directions are sampled from a randomized two-site update, so no Hamiltonian + # or environments are required for i in 1:(length(ψ) - 1) - AC2 = randomize!(_transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])) - - #Calculate nullspaces for left and right - NL = left_null(ψ.AC[i]) - NR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true)) - - #Use this nullspaces and SVD decomposition to determine the optimal expansion space - intermediate = normalize!(adjoint(NL) * AC2 * adjoint(NR)) - _, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd) - - ar_re = Vᴴ * NR - ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(Vᴴ, 1))) - - nal, nc = qr_compact!(catdomain(ψ.AC[i], ar_le)) - nar = _transpose_front(catcodomain(_transpose_tail(ψ.AR[i + 1]), ar_re)) - - ψ.AC[i] = (nal, nc) - ψ.AC[i + 1] = (nc, nar) + changebond!(i, Val(:right), ψ, nothing, alg, nothing) end return normalize!(ψ) end +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs) + bond = site + left = ψ.AC[site] + right = ψ.AR[site + 1] + NL = left_null(left) + NR = right_null!(_transpose_tail(right; copy = true)) + + # randomized two-site update; H and envs are unused + ac2 = randomize!(AC2(ψ, bond)) + + # select the dominant directions in the complement of the current state + g2 = adjoint(NL) * ac2 * adjoint(NR) + ϵ_2site = norm(g2) / norm(ac2) + _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + + # optimal vectors at site+1, zero weight at site + ar_re = Vᴴ * NR + # embed `left` into the enlarged domain (zero weight in the new directions) + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) + + ψ.AC[site] = (nal, normalize!(nc)) + ψ.AC[site + 1] = (nc, nar) + return ψ, (; ϵ_select, ϵ_2site) +end +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs) + bond = site - 1 + left = ψ.AL[site - 1] + right = ψ.AC[site] + NL = left_null(left) + NR = right_null!(_transpose_tail(right; copy = true)) + + # randomized two-site update; H and envs are unused + ac2 = randomize!(AC2(ψ, bond)) + + # select the dominant directions in the complement of the current state + g2 = adjoint(NL) * ac2 * adjoint(NR) + ϵ_2site = norm(g2) / norm(ac2) + U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + + # optimal vectors at site-1, zero weight at site + Q = NL * U + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) + right_tail = _transpose_tail(right) + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) + AL_exp = catdomain(left, Q) + + ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) + ψ.AC[site - 1] = (AL_exp, nc) + return ψ, (; ϵ_select, ϵ_2site) +end + """ sample_space(V, strategy) From ff0e162d9b7029fc65e3b2f0c19db8bc2354effa Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Jun 2026 11:24:02 -0400 Subject: [PATCH 03/20] update CBE DMRG implementation --- src/algorithms/groundstate/dmrg.jl | 177 +++++++++++++++-------------- 1 file changed, 89 insertions(+), 88 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index a4f177b88..017ab3e51 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -78,7 +78,21 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme return ψ, envs, ϵ end -struct CBEDMRG{A, F} <: Algorithm +""" +$(TYPEDEF) + +Single-site DMRG algorithm with Controlled Bond Expansion (CBE): before each single-site +optimization the bond is enriched with directions orthogonal to the current state, so the +optimization can grow the bond dimension while only ever diagonalizing a single-site effective +Hamiltonian. The expansion is delegated to a pluggable bond-expansion sub-algorithm, +allowing different selection strategies (e.g. [`OptimalExpand`](@ref), +[`RandExpand`](@ref)) to be swapped in. + +## Fields + +$(TYPEDFIELDS) +""" +struct CBEDMRG{A, F, E <: Algorithm} <: Algorithm "tolerance for convergence criterium" tol::Float64 @@ -94,11 +108,82 @@ struct CBEDMRG{A, F} <: Algorithm "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" finalize::F - selection::TruncationStrategy + "algorithm used to expand the bond dimension ahead of each single-site optimization" + expand::E "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update" trscheme::TruncationStrategy end +function CBEDMRG(; + tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, + alg_eigsolve = (;), finalize = Defaults._finalize, expand, trscheme + ) + alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : + alg_eigsolve + return CBEDMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, expand, trscheme) +end + +# One CBEDMRG half-sweep: for each bond in the sweep direction, expand the bond with directions +# orthogonal to ψ, optimize the single-site tensor in the enlarged space, then truncate back down +# and move the center. The expansion is shared with `changebonds` via `changebond!`; the +# eigensolve + truncation + fidelity error are CBEDMRG-specific. The extra error measures +# (ϵs_galerkin/ϵs_trunc/ϵs_2site) are diagnostics and do not drive convergence; ϵs is the +# two-site fidelity that does. +function cbe_sweep!( + ψ, H, envs, alg::CBEDMRG, alg_eigsolve, dir::Val{D}, + ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site + ) where {D} + positions = D === :right ? (1:(length(ψ) - 1)) : (length(ψ):-1:2) + for pos in positions + # the two MPS tensors of the pre-update two-site state, kept by reference (the expansion + # builds new tensors rather than mutating these) for the DMRG2-style fidelity error + # below, so the dense two-site tensor is never formed + if D === :right + bra_left, bra_right = ψ.AC[pos], ψ.AR[pos + 1] + bond = pos + else + bra_left, bra_right = ψ.AL[pos - 1], ψ.AC[pos] + bond = pos - 1 + end + + # enrich the bond ahead of the optimization with directions orthogonal to ψ + _, info = changebond!(pos, dir, ψ, H, alg.expand, envs) + ϵs_2site[bond] = max(ϵs_2site[bond], info.ϵ_2site) + # the expanded neighbour tensor, kept for the fidelity contraction below + neighbor = D === :right ? ψ.AR[pos + 1] : ψ.AL[pos - 1] + + # single-site optimization in the enlarged space + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + + # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular values, + # computed directly from the spectrum (no 1 - ‖C‖² cancellation). `v` is the overlap of + # the pre-update two-site state with the truncated, optimized one, contracted directly + # from the MPS factors in zipper order (peak intermediate is a single 3-leg tensor) so + # neither dense two-site tensor is formed; the truncated bond is internal, so the + # fidelity 1 - |v| -> 0 at convergence. + if D === :right + U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) + AL′ = U + C = normalize!(S * Vᴴ) + v = @plansor bra_left[1 2; 7] * conj(AL′[1 2; 5]) * conj(C[5; 6]) * + bra_right[7 4; 3] * conj(neighbor[6 4; 3]) + ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) + ψ.AC[pos] = (AL′, C) + else + U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) + C = normalize!(U * S) + AR_mps = _transpose_front(AR′) + v = @plansor bra_left[1 2; 7] * conj(neighbor[1 2; 5]) * conj(C[5; 6]) * + bra_right[7 4; 3] * conj(AR_mps[6 4; 3]) + ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) + ψ.AC[pos] = (C, AR_mps) + end + ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) + ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) + end + return ψ +end function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg::CBEDMRG, envs = environments(ψ, H)) ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) @@ -118,93 +203,9 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg:: zerovector!(ϵs_galerkin) zerovector!(ϵs_trunc) zerovector!(ϵs_2site) - for pos in 1:(length(ψ) - 1) - @plansor ac2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3] - - # exact two-site update at bond (pos, pos+1), projected onto the complement - # of the current state (this is the exact analog of the randomized SVD) - AC2 = AC2_projection(pos, ψ, H, ψ, envs) - NL = left_null(ψ.AC[pos]) - NR = right_null!(_transpose_tail(ψ.AR[pos + 1]; copy = true)) - g2 = adjoint(NL) * AC2 * adjoint(NR) - # two-site Galerkin: norm of the two-site gradient in the complement, normalized - ϵs_2site[pos] = max(ϵs_2site[pos], norm(g2) / norm(AC2)) - intermediate = normalize!(g2) - _, _, V = svd_trunc!(intermediate; trunc = alg.selection, alg = Defaults.alg_svd()) - - # expand the bond: optimal vectors at pos+1, zero weight at pos - ar_re = V * NR - ar_le = zerovector!(similar(ψ.AC[pos], codomain(ψ.AC[pos]) ← space(V, 1))) - Ql, C = qr_compact!(catdomain(ψ.AC[pos], ar_le)) - AR_exp = _transpose_front(catcodomain(_transpose_tail(ψ.AR[pos + 1]), ar_re)) - ψ.AC[pos] = (Ql, normalize!(C)) - ψ.AC[pos + 1] = (C, AR_exp) - - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) - # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular - # values, computed directly from the spectrum (no 1 - ‖C‖² cancellation) - U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) - ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) - AL′ = U - C = S * Vᴴ - normalize!(C) - - # DMRG2-style error: 1 - fidelity between the two-site state before the - # eigensolver and the truncated, optimized two-site state after it. The - # truncated bond is internal, so this -> 0 at convergence. - @plansor ac2′[-1 -2; -3 -4] := AL′[-1 -2; 1] * C[1; 2] * AR_exp[2 -4; -3] - ϵs[pos] = max(ϵs[pos], abs(1 - abs(dot(ac2, ac2′)))) - - ψ.AC[pos] = (AL′, C) - ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) - # @debug "CBEDMRG L→R" pos ϵ = ϵs[pos] - end - for pos in length(ψ):-1:2 - @plansor ac2[-1 -2; -3 -4] := ψ.AL[pos - 1][-1 -2; 1] * ψ.AC[pos][1 -4; -3] - - # exact two-site update at bond (pos-1, pos), projected onto the complement - AC2 = AC2_projection(pos - 1, ψ, H, ψ, envs) - NL = left_null(ψ.AL[pos - 1]) - NR = right_null!(_transpose_tail(ψ.AC[pos]; copy = true)) - g2 = adjoint(NL) * AC2 * adjoint(NR) - # two-site Galerkin: norm of the two-site gradient in the complement, normalized - ϵs_2site[pos - 1] = max(ϵs_2site[pos - 1], norm(g2) / norm(AC2)) - intermediate = normalize!(g2) - U, _, _ = svd_trunc!(intermediate; trunc = alg.selection, alg = Defaults.alg_svd()) - - # optimal new left vectors for AL[pos-1]; zero padding for AC[pos]'s left virtual - Q = NL * U - AL = ψ.AL[pos - 1] - al_le = zerovector!( - similar(ψ.AC[pos], space(Q, 3)' ⊗ physicalspace(ψ.AC[pos]) ← right_virtualspace(ψ.AC[pos])) - ) - C, Qr = lq_compact!(catcodomain(_transpose_tail(ψ.AC[pos]), _transpose_tail(al_le))) - AL_exp = catdomain(AL, Q) - ψ.AC[pos] = (normalize!(C), _transpose_front(Qr)) - ψ.AC[pos - 1] = (AL_exp, C) - - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) - # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular - # values, computed directly from the spectrum (no 1 - ‖C‖² cancellation) - U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) - ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) - C = U * S - normalize!(C) - AR_mps = _transpose_front(AR′) - - # DMRG2-style error: 1 - fidelity between the two-site state before the - # eigensolver and the truncated, optimized two-site state after it. The - # truncated bond is internal, so this -> 0 at convergence. AR_mps is in MPS - # form (physical in codomain) so the overlap with ac2 stays planar. - @plansor ac2′[-1 -2; -3 -4] := AL_exp[-1 -2; 1] * C[1; 2] * AR_mps[2 -4; -3] - ϵs[pos] = max(ϵs[pos], abs(1 - abs(dot(ac2, ac2′)))) - # @debug "CBEDMRG R→L" pos ϵ = ϵs[pos] - ψ.AC[pos] = (C, AR_mps) - ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) - end + cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:right), ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site) + cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:left), ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site) ϵ = maximum(ϵs) @infov 2 "CBEDMRG errors" iter fidelity = ϵ local_galerkin = maximum(ϵs_galerkin) twosite_galerkin = maximum(ϵs_2site) truncation = maximum(ϵs_trunc) From a076e572676c3e9f3ff006448f3b6d02d04006af Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Jun 2026 11:24:12 -0400 Subject: [PATCH 04/20] add CBEDMRG test --- test/algorithms/groundstate.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index cb1faa261..766bfc3a7 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -61,6 +61,31 @@ verbosity_conv = 1 @test v < 1.0e-2 end + @testset "CBEDMRG" begin + # start from a small bond so the bond expansion is exercised + ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) + v₀ = variance(ψ₀, H) + expand = OptimalExpand(; trscheme = truncrank(D ÷ 2)) + trscheme = truncrank(D) + + # test logging + ψ, envs, δ = find_groundstate( + ψ₀, H, CBEDMRG(; verbosity = verbosity_full, maxiter = 2, expand, trscheme) + ) + + ψ, envs, δ = find_groundstate( + ψ, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 10, expand, trscheme), envs + ) + v = variance(ψ, H) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1.0e-3 + @test v < v₀ + @test v < 1.0e-2 + # the bond should have grown to the truncation target + @test dim(left_virtualspace(ψ, L ÷ 2)) == D + end + @testset "GradientGrassmann" begin ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D) v₀ = variance(ψ₀, H) From 5ea3d2f1455685b3dcb49505d76194b1af076276 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 19 Jun 2026 17:19:28 -0400 Subject: [PATCH 05/20] variety of updates MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Reconcile CBE with main after rebase: - migrate CBEDMRG/SketchedExpand to the 3-arg environments(ψ, H, ψ) API (#436 removed the defaulted operator-environments form) - instrument CBEDMRG/cbe_sweep! with TimerOutputs (expand / AC_eigsolve / svd_trunc / sweep / finalize), matching the #431 DMRG/DMRG2 style Co-Authored-By: Claude Opus 4.8 (1M context) --- src/MPSKit.jl | 3 +- src/algorithms/changebonds/changebonds.jl | 5 +- src/algorithms/changebonds/optimalexpand.jl | 63 ++++-- src/algorithms/changebonds/randexpand.jl | 19 +- src/algorithms/changebonds/sketchedexpand.jl | 201 +++++++++++++++++++ src/algorithms/groundstate/dmrg.jl | 86 ++++---- test/algorithms/changebonds.jl | 16 ++ test/algorithms/groundstate.jl | 42 ++++ 8 files changed, 368 insertions(+), 67 deletions(-) create mode 100644 src/algorithms/changebonds/sketchedexpand.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 03d54c215..4c999cf58 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -37,7 +37,7 @@ export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo export TDVP, TDVP2, WI, WII, TaylorCluster export changebonds, changebonds! -export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand +export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann export exact_diagonalization, fidelity_susceptibility @@ -156,6 +156,7 @@ include("algorithms/changebonds/optimalexpand.jl") include("algorithms/changebonds/vumpssvd.jl") include("algorithms/changebonds/svdcut.jl") include("algorithms/changebonds/randexpand.jl") +include("algorithms/changebonds/sketchedexpand.jl") include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/taylorcluster.jl") diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl index 56de20f1c..c8b8bd17a 100644 --- a/src/algorithms/changebonds/changebonds.jl +++ b/src/algorithms/changebonds/changebonds.jl @@ -13,10 +13,11 @@ function changebonds end function changebonds! end @doc """ - changebond(site, dir, ψ, [H], alg, [envs]) -> ψ, info - changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ, info + changebond(site, dir, ψ, [H], alg, [envs]) -> ψ + changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ Expand a single bond of `ψ` in place by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the enriched bond. +Diagnostic error measures of the expansion (e.g. the discarded selection weight) are reported through the logging system at `@infov 4` rather than returned. The sweep direction `dir` is a `Val(:right)` or `Val(:left)` used for dispatch. For `Val(:right)` the bond `(site, site + 1)` is enriched on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions; for `Val(:left)` the mirror is applied to bond `(site - 1, site)`. diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 28b46f467..9b7bd0ca3 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -5,6 +5,15 @@ An algorithm that expands the given mps as described in [Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018), by selecting the dominant contributions of a two-site updated MPS tensor, orthogonal to the original ψ. +By default the expansion does not alter the state: the added directions are connected through a +zero block, so that the expanded state is identical to the original one (as required for e.g. +TDVP). When `warmstart = true`, the new directions are instead seeded with the (physically +scaled) two-site gradient itself, so the expanded state already moves toward the optimal +two-site update. This changes the state and is therefore only useful for ground-state search +(e.g. as the `expand` strategy of [`CBEDMRG`](@ref)), where it warm-starts the subsequent +single-site optimization; the injected amplitude scales with the gradient norm and so vanishes +automatically at convergence. + !!! note [`changebonds!`](@ref) is only defined for `FiniteMPS`, and modifies both the state and its environment. @@ -18,6 +27,12 @@ $(TYPEDFIELDS) "algorithm used for truncating the expanded space" trscheme::TruncationStrategy + + "whether to seed the new directions with the two-site gradient (warm start, alters the state) instead of a zero block" + warmstart::Bool = false + + "setting for how much information is displayed" + verbosity::Int = Defaults.verbosity end # Simple wrapper to convert between diffrent type of InifniteMPS. @@ -96,19 +111,26 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) - ϵ_2site = norm(g2) / norm(AC2) - _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + gradnorm = norm(g2) + U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + @infov 4 "bond expansion" site dir = :right ϵ_select ϵ_2site = gradnorm / norm(AC2) - # optimal vectors at site+1, zero weight at site + # optimal vectors at site+1 ar_re = Vᴴ * NR - # embed `left` into the enlarged domain (zero weight in the new directions) - nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) - nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + if alg.warmstart + # seed the new left directions with the (physically scaled) gradient instead of a zero + # block, warm-starting the subsequent optimization (alters the state) + nal, nc = qr_compact!(catdomain(left, gradnorm * (NL * U * S))) + else + # embed `left` into the enlarged domain (zero weight in the new directions) + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) ψ.AC[site] = (nal, normalize!(nc)) ψ.AC[site + 1] = (nc, nar) - return ψ, (; ϵ_select, ϵ_2site) + return ψ end function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs) bond = site - 1 @@ -122,28 +144,37 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Opt # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) - ϵ_2site = norm(g2) / norm(AC2) - U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + gradnorm = norm(g2) + U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + @infov 4 "bond expansion" site dir = :left ϵ_select ϵ_2site = gradnorm / norm(AC2) - # optimal vectors at site-1, zero weight at site + # optimal vectors at site-1 Q = NL * U - # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) right_tail = _transpose_tail(right) - nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) - nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) + if alg.warmstart + # seed the new right directions with the (physically scaled) gradient instead of a zero + # block, warm-starting the subsequent optimization (alters the state) + nc, Qr = lq_compact!(catcodomain(right_tail, gradnorm * (S * Vᴴ * NR))) + else + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) + end AL_exp = catdomain(left, Q) ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) ψ.AC[site - 1] = (AL_exp, nc) - return ψ, (; ϵ_select, ϵ_2site) + return ψ end changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) = changebonds!(copy(ψ), H, alg, envs) function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) - for i in 1:(length(ψ) - 1) - changebond!(i, Val(:right), ψ, H, alg, envs) + LoggingExtras.withlevel(; alg.verbosity) do + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, H, alg, envs) + end end return ψ, envs end diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index af849c28c..ad84a66af 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -23,6 +23,9 @@ $(TYPEDFIELDS) "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" trscheme::TruncationStrategy + + "setting for how much information is displayed" + verbosity::Int = Defaults.verbosity end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) @@ -70,8 +73,10 @@ end function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) # the expansion directions are sampled from a randomized two-site update, so no Hamiltonian # or environments are required - for i in 1:(length(ψ) - 1) - changebond!(i, Val(:right), ψ, nothing, alg, nothing) + LoggingExtras.withlevel(; alg.verbosity) do + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, nothing, alg, nothing) + end end return normalize!(ψ) @@ -89,8 +94,9 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra # select the dominant directions in the complement of the current state g2 = adjoint(NL) * ac2 * adjoint(NR) - ϵ_2site = norm(g2) / norm(ac2) + gradnorm = norm(g2) _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + @infov 4 "bond expansion" site dir = :right ϵ_select ϵ_2site = gradnorm / norm(ac2) # optimal vectors at site+1, zero weight at site ar_re = Vᴴ * NR @@ -101,7 +107,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra ψ.AC[site] = (nal, normalize!(nc)) ψ.AC[site + 1] = (nc, nar) - return ψ, (; ϵ_select, ϵ_2site) + return ψ end function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs) bond = site - 1 @@ -115,8 +121,9 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ran # select the dominant directions in the complement of the current state g2 = adjoint(NL) * ac2 * adjoint(NR) - ϵ_2site = norm(g2) / norm(ac2) + gradnorm = norm(g2) U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + @infov 4 "bond expansion" site dir = :left ϵ_select ϵ_2site = gradnorm / norm(ac2) # optimal vectors at site-1, zero weight at site Q = NL * U @@ -128,7 +135,7 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ran ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) ψ.AC[site - 1] = (AL_exp, nc) - return ψ, (; ϵ_select, ϵ_2site) + return ψ end """ diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl new file mode 100644 index 000000000..cf98ccd20 --- /dev/null +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -0,0 +1,201 @@ +""" +$(TYPEDEF) + +An algorithm that expands the bond dimension by selecting the dominant directions of the +projected two-site update orthogonal to the current state, just like [`OptimalExpand`](@ref), +but at single-site cost using a randomized range-finder (the "shrewd selection" of Controlled +Bond Expansion, Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)). + +Rather than forming the full projected two-site update (which requires applying the two-site +effective Hamiltonian), a random sketch of width `rank + oversampling` of the orthogonal +complement is folded into the effective *environment* first, collapsing the large bond before +the two-site contraction is ever materialized. The dominant directions are then recovered from +a small singular value decomposition of the sketched gradient. + +!!! note + This strategy is only defined for `FiniteMPS` (through [`changebond!`](@ref)), so that it + can be used both standalone and as the `expand` strategy of [`CBEDMRG`](@ref). + +Like [`OptimalExpand`](@ref), the expansion is by default state-preserving (the added +directions are connected through a zero block, as required for e.g. TDVP). When +`warmstart = true`, the new directions are instead seeded with the (physically scaled) sketched +two-site gradient, warm-starting the subsequent single-site optimization in ground-state search +(e.g. as the `expand` strategy of [`CBEDMRG`](@ref)); this alters the state, and the injected +amplitude scales with the gradient norm so that it vanishes automatically at convergence. + +!!! note + `ϵ_2site` is a randomized estimate of the two-site complement-gradient norm (a diagnostic + only), not the exact ratio `‖g2‖ / ‖AC2‖` reported by [`OptimalExpand`](@ref). The folded + application routes the `MPOHamiltonian` tensors through the generic transfer/one-site + machinery and therefore does not exploit `JordanMPO` sparsity in the MPO bond dimension. + +## Fields + +$(TYPEDFIELDS) +""" +@kwdef struct SketchedExpand{S} <: Algorithm + "algorithm used for the (small) singular value decomposition" + alg_svd::S = Defaults.alg_svd() + + "algorithm used for truncating the expanded space" + trscheme::TruncationStrategy + + "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)" + oversampling::Int = 5 + + "whether to seed the new directions with the sketched two-site gradient (warm start, alters the state) instead of a zero block" + warmstart::Bool = false + + "setting for how much information is displayed" + verbosity::Int = Defaults.verbosity +end + +""" + sketch_space(V, alg::SketchedExpand) + +Determine the space of the random sketch drawn from the complement space `V`: the target +subspace selected by `alg.trscheme`, enlarged by `alg.oversampling` extra directions (capped +by `V`). Reuses [`sample_space`](@ref) so the result is sector-correct for graded spaces. +""" +function sketch_space(V, alg::SketchedExpand) + Vk = sample_space(V, alg.trscheme) + alg.oversampling == 0 && return Vk + Vp = sample_space(V, MatrixAlgebraKit.truncrank(alg.oversampling)) + return infimum(V, Vk ⊕ Vp) +end + +""" + complement_space(Vfull::ElementarySpace, image::ElementarySpace) + +The orthogonal complement `Vfull ⊖ image` of the range `image` inside the fused space `Vfull`, +computed sector-wise from the spaces alone (no factorization). For an isometry `A : Vfull ← image` +this is the space carried by `left_null(A)`; here `Vfull` is the fused codomain/domain of the +isometric MPS tensor and `image` its small virtual leg. Used to size and cap the random sketch. +""" +function complement_space(Vfull::ElementarySpace, image::ElementarySpace) + @assert !isdual(Vfull) + pairs = [c => (dim(Vfull, c) - dim(image, c)) for c in sectors(Vfull)] + filter!(p -> last(p) > 0, pairs) + return typeof(Vfull)(pairs) +end + +# Finite system +# ------------- +# The "shrewd selection" needs the orthogonal complement of the current state on either side of +# the bond. Rather than materializing an explicit complement basis with `left_null`/`right_null`, +# we use the complement projector of the *isometric* MPS tensor `A`: `I - A A'` (left-iso) and +# `I - A' A` (right-iso). The isometric tensor is obtained from the center tensor by a local +# `left_orth`/`right_orth!`, which leaves the MPS gauge (and any environments a caller maintains +# incrementally, e.g. `CBEDMRG`) untouched. On the sketch side this turns the range-finder into a +# narrow random draw projected into the complement; on the selection side the gradient is projected +# and its dominant directions read off directly from the SVD (its singular vectors land in the +# complement automatically). The kept rank is capped at the complement dimension: the projected +# object has only that many genuine singular values, and at edge bonds the complement may be empty, +# where an uncapped `svd_trunc` would otherwise pad the selection with non-orthogonal noise. +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs) + left = ψ.AC[site] + right = ψ.AR[site + 1] + AL, _ = left_orth(left) # local left-isometric form (range(AL) = range(left)) + ARtt = _transpose_tail(right) # right-isometric: ARtt * ARtt' = I + + # nothing to add when either complement is empty (e.g. edge bonds) + compL = complement_space(fuse(codomain(AL)), only(domain(AL))) + compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) + (dim(compL) == 0 || dim(compR) == 0) && return ψ + + # sketch the left complement: Q = (I - AL AL') Ω lies in the complement of `left`, re-orthonormalized + # by a narrow QR into a proper randomized range-finder basis (no full complement basis is formed) + Vℓ = sketch_space(compL, alg) + Ω = randisometry(scalartype(left), codomain(AL) ← Vℓ) + Q, _ = qr_compact!(Ω - AL * (AL' * Ω)) + + # fold the sketch into the left environment (bra-leg becomes the sketch leg), then apply the + # existing single-site effective Hamiltonian on site+1: this reconstructs `Q† AC2_projection` + # without ever forming the full two-site update + GL = leftenv(envs, site, ψ) * TransferMatrix(left, H[site], Q) + Hac = MPO_AC_Hamiltonian(GL, H[site + 1], rightenv(envs, site + 1, ψ)) + Y = Hac * right + + # project onto the right complement with (I - AR' AR): the SVD right-vectors lie in the complement + Ytt = _transpose_tail(Y) + B = Ytt - (Ytt * ARtt') * ARtt + nrm = norm(B) + trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compR)) + U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) + @infov 4 "bond expansion (sketched)" site dir = :right ϵ_select ϵ_2site = nrm + + # optimal vectors at site+1 (identical bookkeeping to OptimalExpand) + ar_re = Vᴴ + if alg.warmstart + # seed the new left directions with the (physically scaled) sketched gradient instead of + # a zero block; `Q * U` maps the sketch-basis left vectors back to the full space + nal, nc = qr_compact!(catdomain(left, nrm * (Q * U * S))) + else + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + end + nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) + + ψ.AC[site] = (nal, normalize!(nc)) + ψ.AC[site + 1] = (nc, nar) + return ψ +end +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs) + left = ψ.AL[site - 1] # left-isometric, left of center (valid) + right = ψ.AC[site] + _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc()) # local right-iso + + # nothing to add when either complement is empty (e.g. edge bonds) + compL = complement_space(fuse(codomain(left)), only(domain(left))) + compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) + (dim(compL) == 0 || dim(compR) == 0) && return ψ + + # sketch the right complement: Qr = (I - AR' AR) Ω lies in the complement of `right`, with + # orthonormal rows recovered by a narrow LQ (the randomized range-finder basis) + Vℓ = sketch_space(compR, alg) + Ω = adjoint(randisometry(scalartype(right), domain(ARtt) ← Vℓ)) + _, Qr_o = lq_compact!(Ω - (Ω * ARtt') * ARtt) + Qr = _transpose_front(Qr_o) + + # fold the sketch into the right environment, then apply the single-site effective + # Hamiltonian on site-1: reconstructs `AC2_projection NR†` (sketched) at single-site cost + GR = TransferMatrix(right, H[site], Qr) * rightenv(envs, site, ψ) + Hac = MPO_AC_Hamiltonian(leftenv(envs, site - 1, ψ), H[site - 1], GR) + Y = Hac * left + + # project onto the left complement with (I - AL AL'): the SVD left-vectors lie in the complement + B = Y - left * (left' * Y) + nrm = norm(B) + trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compL)) + U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) + @infov 4 "bond expansion (sketched)" site dir = :left ϵ_select ϵ_2site = nrm + + # optimal vectors at site-1 (identical bookkeeping to OptimalExpand) + Q = U + right_tail = _transpose_tail(right) + if alg.warmstart + # seed the new right directions with the (physically scaled) sketched gradient instead of + # a zero block; `Vᴴ * Qr_o` maps the sketch-basis right vectors back to the full space + nc, Qr2 = lq_compact!(catcodomain(right_tail, nrm * (S * Vᴴ * Qr_o))) + else + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) + end + AL_exp = catdomain(left, Q) + + ψ.AC[site] = (normalize!(nc), _transpose_front(Qr2)) + ψ.AC[site - 1] = (AL_exp, nc) + return ψ +end + +changebonds(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments(ψ, H, ψ)) = + changebonds!(copy(ψ), H, alg, envs) + +function changebonds!(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments(ψ, H, ψ)) + LoggingExtras.withlevel(; alg.verbosity) do + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, H, alg, envs) + end + end + return ψ, envs +end diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 017ab3e51..35c1abc84 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -126,35 +126,35 @@ end # One CBEDMRG half-sweep: for each bond in the sweep direction, expand the bond with directions # orthogonal to ψ, optimize the single-site tensor in the enlarged space, then truncate back down # and move the center. The expansion is shared with `changebonds` via `changebond!`; the -# eigensolve + truncation + fidelity error are CBEDMRG-specific. The extra error measures -# (ϵs_galerkin/ϵs_trunc/ϵs_2site) are diagnostics and do not drive convergence; ϵs is the -# two-site fidelity that does. +# eigensolve + truncation + fidelity error are CBEDMRG-specific. Only `ϵs` (the two-site +# fidelity) is threaded back to drive convergence; the per-bond diagnostics (truncation, +# local Galerkin, and the expansion's own ϵ_select/ϵ_2site) are emitted via `@infov 4`. function cbe_sweep!( - ψ, H, envs, alg::CBEDMRG, alg_eigsolve, dir::Val{D}, - ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site + ψ, H, envs, alg::CBEDMRG, alg_eigsolve, dir::Val{D}, ϵs, timeroutput ) where {D} positions = D === :right ? (1:(length(ψ) - 1)) : (length(ψ):-1:2) for pos in positions + local AC′, ϵ_trunc, v # the two MPS tensors of the pre-update two-site state, kept by reference (the expansion # builds new tensors rather than mutating these) for the DMRG2-style fidelity error # below, so the dense two-site tensor is never formed if D === :right bra_left, bra_right = ψ.AC[pos], ψ.AR[pos + 1] - bond = pos else bra_left, bra_right = ψ.AL[pos - 1], ψ.AC[pos] - bond = pos - 1 end # enrich the bond ahead of the optimization with directions orthogonal to ψ - _, info = changebond!(pos, dir, ψ, H, alg.expand, envs) - ϵs_2site[bond] = max(ϵs_2site[bond], info.ϵ_2site) + # (the expansion logs its own ϵ_select/ϵ_2site diagnostics at @infov 4) + @timeit timeroutput "expand" changebond!(pos, dir, ψ, H, alg.expand, envs) # the expanded neighbour tensor, kept for the fidelity contraction below neighbor = D === :right ? ψ.AR[pos + 1] : ψ.AL[pos - 1] # single-site optimization in the enlarged space - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + @timeit timeroutput "AC_eigsolve" begin + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + end # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular values, # computed directly from the spectrum (no 1 - ‖C‖² cancellation). `v` is the overlap of @@ -162,37 +162,36 @@ function cbe_sweep!( # from the MPS factors in zipper order (peak intermediate is a single 3-leg tensor) so # neither dense two-site tensor is formed; the truncated bond is internal, so the # fidelity 1 - |v| -> 0 at convergence. - if D === :right - U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) - AL′ = U - C = normalize!(S * Vᴴ) - v = @plansor bra_left[1 2; 7] * conj(AL′[1 2; 5]) * conj(C[5; 6]) * - bra_right[7 4; 3] * conj(neighbor[6 4; 3]) - ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) - ψ.AC[pos] = (AL′, C) - else - U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) - C = normalize!(U * S) - AR_mps = _transpose_front(AR′) - v = @plansor bra_left[1 2; 7] * conj(neighbor[1 2; 5]) * conj(C[5; 6]) * - bra_right[7 4; 3] * conj(AR_mps[6 4; 3]) - ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) - ψ.AC[pos] = (C, AR_mps) + @timeit timeroutput "svd_trunc" begin + if D === :right + U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) + AL′ = U + C = normalize!(S * Vᴴ) + v = @plansor bra_left[1 2; 7] * conj(AL′[1 2; 5]) * conj(C[5; 6]) * + bra_right[7 4; 3] * conj(neighbor[6 4; 3]) + ψ.AC[pos] = (AL′, C) + else + U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) + C = normalize!(U * S) + AR_mps = _transpose_front(AR′) + v = @plansor bra_left[1 2; 7] * conj(neighbor[1 2; 5]) * conj(C[5; 6]) * + bra_right[7 4; 3] * conj(AR_mps[6 4; 3]) + ψ.AC[pos] = (C, AR_mps) + end end - ϵs_trunc[pos] = max(ϵs_trunc[pos], ϵ_trunc) - ϵs_galerkin[pos] = max(ϵs_galerkin[pos], calc_galerkin(pos, ψ, H, ψ, envs)) + fid = abs(1 - abs(v)) + ϵs[pos] = max(ϵs[pos], fid) + @infov 4 "CBEDMRG bond" pos fidelity = fid truncation = ϵ_trunc galerkin = calc_galerkin(pos, ψ, H, ψ, envs) end return ψ end -function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg::CBEDMRG, envs = environments(ψ, H)) +function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg::CBEDMRG, envs = environments(ψ, H, ψ)) ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) ϵ = maximum(ϵs) - # extra debug error measures (do not drive convergence; ϵ remains the two-site fidelity) - ϵs_galerkin = zero(ϵs) # local (one-site) Galerkin error - ϵs_trunc = zero(ϵs) # bond-truncation error (discarded weight) - ϵs_2site = zero(ϵs) # two-site Galerkin error (complement-projected two-site gradient) - log = IterLog("DMRG") + log = IterLog("CBEDMRG") + timeroutput = TimerOutput("CBEDMRG") + alg.verbosity > 3 || disable_timer!(timeroutput) LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs)) @@ -200,23 +199,25 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg:: alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) zerovector!(ϵs) - zerovector!(ϵs_galerkin) - zerovector!(ϵs_trunc) - zerovector!(ϵs_2site) - cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:right), ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site) - cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:left), ϵs, ϵs_galerkin, ϵs_trunc, ϵs_2site) + @timeit timeroutput "sweep" begin + cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:right), ϵs, timeroutput) + cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:left), ϵs, timeroutput) + end ϵ = maximum(ϵs) - @infov 2 "CBEDMRG errors" iter fidelity = ϵ local_galerkin = maximum(ϵs_galerkin) twosite_galerkin = maximum(ϵs_2site) truncation = maximum(ϵs_trunc) - ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} + ψ, envs = @timeit timeroutput "finalize" alg.finalize( + iter, ψ, H, envs + )::Tuple{typeof(ψ), typeof(envs)} if ϵ <= alg.tol + @infov 4 timeroutput @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs)) break end if iter == alg.maxiter + @infov 4 timeroutput @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs)) else @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs)) @@ -321,6 +322,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm ψ.AC[pos] = (al, complex(c)) end end + end ϵ = maximum(ϵs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl index 0bf95b7cf..958df806d 100644 --- a/test/algorithms/changebonds.jl +++ b/test/algorithms/changebonds.jl @@ -84,6 +84,22 @@ end ) @test dot(state, state_oe) ≈ 1 atol = 1.0e-8 + state_se, _ = changebonds( + state, H, + SketchedExpand(; trscheme = truncrank(dim(Dspace) * dim(Dspace)), oversampling = 4) + ) + @test dot(state, state_se) ≈ 1 atol = 1.0e-8 + + # `warmstart` seeds the new directions with the two-site gradient: the expansion is no longer + # norm-preserving, but the resulting state is still normalized + for (Exp, kw) in ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) + state_ws, _ = changebonds( + state, H, Exp(; trscheme = truncrank(dim(Dspace) * dim(Dspace)), warmstart = true, kw...) + ) + @test dot(state_ws, state_ws) ≈ 1 atol = 1.0e-8 + @test !isapprox(abs(dot(state, state_ws)), 1; atol = 1.0e-4) + end + state_tr = changebonds(state_oe, SvdCut(; trscheme = truncrank(dim(Dspace)))) @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5)) diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index 766bfc3a7..05bb693d5 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -9,6 +9,7 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ +using Random verbosity_full = 5 verbosity_conv = 1 @@ -86,6 +87,47 @@ verbosity_conv = 1 @test dim(left_virtualspace(ψ, L ÷ 2)) == D end + @testset "CBEDMRG (SketchedExpand)" begin + # randomized bond expansion at single-site cost + Random.seed!(1234) + ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) + v₀ = variance(ψ₀, H) + expand = SketchedExpand(; trscheme = truncrank(D ÷ 2), oversampling = 4) + trscheme = truncrank(D) + + # test logging + ψ, envs, δ = find_groundstate( + ψ₀, H, CBEDMRG(; verbosity = verbosity_full, maxiter = 2, expand, trscheme) + ) + + ψ, envs, δ = find_groundstate( + ψ, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 15, expand, trscheme), envs + ) + v = variance(ψ, H) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1.0e-3 + @test v < v₀ + @test v < 1.0e-2 + # the bond should have grown to the truncation target + @test dim(left_virtualspace(ψ, L ÷ 2)) == D + end + + @testset "CBEDMRG warmstart $(nameof(Exp))" for (Exp, kw) in + ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) + # warmstart seeds the expansion with the two-site gradient (alters the state); the + # ground-state search should still converge to the correct low-variance state + Random.seed!(2025) + ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) + expand = Exp(; trscheme = truncrank(D ÷ 2), warmstart = true, kw...) + ψ, envs, δ = find_groundstate( + ψ₀, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 15, expand, trscheme = truncrank(D)) + ) + @test sum(δ) ≈ 0 atol = 1.0e-3 + @test variance(ψ, H) < 1.0e-2 + @test dim(left_virtualspace(ψ, L ÷ 2)) == D + end + @testset "GradientGrassmann" begin ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D) v₀ = variance(ψ₀, H) From a1c391ffaccb49b90db0445e011214935e1e8aae Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 14:38:29 -0400 Subject: [PATCH 06/20] clean up CBEDMRG --- src/MPSKit.jl | 2 +- src/algorithms/changebonds/optimalexpand.jl | 2 +- src/algorithms/changebonds/sketchedexpand.jl | 6 +- src/algorithms/groundstate/dmrg.jl | 221 ++++++------------- test/algorithms/groundstate.jl | 10 +- 5 files changed, 77 insertions(+), 164 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 4c999cf58..b29174216 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -31,7 +31,7 @@ export leftenv, rightenv export find_groundstate, find_groundstate! export leading_boundary export approximate, approximate! -export VUMPS, VOMPS, DMRG, DMRG2, CBEDMRG, IDMRG, IDMRG2, GradientGrassmann +export VUMPS, VOMPS, DMRG, DMRG2, IDMRG, IDMRG2, GradientGrassmann export excitations export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2 export time_evolve, timestep, timestep!, make_time_mpo diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 9b7bd0ca3..3d2d15db1 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -10,7 +10,7 @@ zero block, so that the expanded state is identical to the original one (as requ TDVP). When `warmstart = true`, the new directions are instead seeded with the (physically scaled) two-site gradient itself, so the expanded state already moves toward the optimal two-site update. This changes the state and is therefore only useful for ground-state search -(e.g. as the `expand` strategy of [`CBEDMRG`](@ref)), where it warm-starts the subsequent +(e.g. as the `alg_expand` strategy of [`DMRG`](@ref)), where it warm-starts the subsequent single-site optimization; the injected amplitude scales with the gradient norm and so vanishes automatically at convergence. diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index cf98ccd20..3f33f58ad 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -14,13 +14,13 @@ a small singular value decomposition of the sketched gradient. !!! note This strategy is only defined for `FiniteMPS` (through [`changebond!`](@ref)), so that it - can be used both standalone and as the `expand` strategy of [`CBEDMRG`](@ref). + can be used both standalone and as the `alg_expand` strategy of [`DMRG`](@ref). Like [`OptimalExpand`](@ref), the expansion is by default state-preserving (the added directions are connected through a zero block, as required for e.g. TDVP). When `warmstart = true`, the new directions are instead seeded with the (physically scaled) sketched two-site gradient, warm-starting the subsequent single-site optimization in ground-state search -(e.g. as the `expand` strategy of [`CBEDMRG`](@ref)); this alters the state, and the injected +(e.g. as the `alg_expand` strategy of [`DMRG`](@ref)); this alters the state, and the injected amplitude scales with the gradient norm so that it vanishes automatically at convergence. !!! note @@ -86,7 +86,7 @@ end # we use the complement projector of the *isometric* MPS tensor `A`: `I - A A'` (left-iso) and # `I - A' A` (right-iso). The isometric tensor is obtained from the center tensor by a local # `left_orth`/`right_orth!`, which leaves the MPS gauge (and any environments a caller maintains -# incrementally, e.g. `CBEDMRG`) untouched. On the sketch side this turns the range-finder into a +# incrementally, e.g. CBE `DMRG`) untouched. On the sketch side this turns the range-finder into a # narrow random draw projected into the complement; on the selection side the gradient is projected # and its dominant directions read off directly from the SVD (its singular vectors land in the # complement automatically). The kept rank is capped at the complement dimension: the projected diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 35c1abc84..bf8c8451d 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -1,13 +1,31 @@ +# whether the gauge algorithm truncates the bond (SVD-based) or preserves it (QR-based, a plain +# center-move as in textbook single-site DMRG) +_truncates(::MatrixAlgebraKit.AbstractAlgorithm) = false +_truncates(::MatrixAlgebraKit.TruncatedAlgorithm) = true + """ $(TYPEDEF) Single-site DMRG algorithm for finding the dominant eigenvector. +Each site update is, in order: (1) an optional bond expansion (`alg_expand`), (2) a single-site +eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and a +non-truncating QR `alg_gauge`) this is textbook single-site DMRG, which cannot change the bond +dimension. Setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref), +[`RandExpand`](@ref), [`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to +the current state ahead of each eigensolve, recovering Controlled Bond Expansion (CBE) DMRG; a +truncating `alg_gauge` is then desirable to cut the enlarged bond back down. + +The gauge algorithm is selected in the keyword constructor from the `trscheme` argument: when it +is `notrunc()` the gauge is a QR decomposition (`alg_orth`, [`Householder`](@extref +MatrixAlgebraKit.Householder) by default), otherwise it is a truncated SVD (`alg_svd` with the +given `trscheme`). + ## Fields $(TYPEDFIELDS) """ -struct DMRG{A, F} <: Algorithm +struct DMRG{A, F, E, G} <: Algorithm "tolerance for convergence criterium" tol::Float64 @@ -22,14 +40,28 @@ struct DMRG{A, F} <: Algorithm "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" finalize::F + + "algorithm used to expand the bond ahead of each local update, or `nothing` for none" + alg_expand::E + + "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD" + alg_gauge::G end function DMRG(; tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;), - verbosity = Defaults.verbosity, finalize = Defaults._finalize + verbosity = Defaults.verbosity, finalize = Defaults._finalize, + alg_expand = nothing, trscheme = notrunc(), + alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth() ) alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : alg_eigsolve - return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize) + # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD + alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth : + MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme) + if !isnothing(alg_expand) && !_truncates(alg_gauge) + @warn "DMRG with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep." + end + return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge) end function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ)) @@ -46,165 +78,46 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environme zerovector!(ϵs) @timeit timeroutput "sweep" begin - for pos in [1:(length(ψ) - 1); length(ψ):-1:2] - local vec - @timeit timeroutput "AC_eigsolve" begin - h = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, vec = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) - end + # left-to-right + for pos in 1:(length(ψ) - 1) + local AC′ + # convergence: pre-expansion single-site Galerkin error ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs)) - @timeit timeroutput "AC_update" ψ.AC[pos] = vec - end - end - ϵ = maximum(ϵs) - - ψ, envs = @timeit timeroutput "finalize" alg.finalize( - iter, ψ, H, envs - )::Tuple{typeof(ψ), typeof(envs)} - - if ϵ <= alg.tol - @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs)) - break - end - if iter == alg.maxiter - @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs)) - else - @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs)) - end - end - end - return ψ, envs, ϵ -end - -""" -$(TYPEDEF) - -Single-site DMRG algorithm with Controlled Bond Expansion (CBE): before each single-site -optimization the bond is enriched with directions orthogonal to the current state, so the -optimization can grow the bond dimension while only ever diagonalizing a single-site effective -Hamiltonian. The expansion is delegated to a pluggable bond-expansion sub-algorithm, -allowing different selection strategies (e.g. [`OptimalExpand`](@ref), -[`RandExpand`](@ref)) to be swapped in. - -## Fields - -$(TYPEDFIELDS) -""" -struct CBEDMRG{A, F, E <: Algorithm} <: Algorithm - "tolerance for convergence criterium" - tol::Float64 - - "maximal amount of iterations" - maxiter::Int - - "setting for how much information is displayed" - verbosity::Int - - "algorithm used for the eigenvalue solvers" - alg_eigsolve::A - - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" - finalize::F - "algorithm used to expand the bond dimension ahead of each single-site optimization" - expand::E + # 1. expand + isnothing(alg.alg_expand) || + @timeit timeroutput "expand" changebond!(pos, Val(:right), ψ, H, alg.alg_expand, envs) - "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update" - trscheme::TruncationStrategy -end -function CBEDMRG(; - tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, - alg_eigsolve = (;), finalize = Defaults._finalize, expand, trscheme - ) - alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : - alg_eigsolve - return CBEDMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, expand, trscheme) -end - -# One CBEDMRG half-sweep: for each bond in the sweep direction, expand the bond with directions -# orthogonal to ψ, optimize the single-site tensor in the enlarged space, then truncate back down -# and move the center. The expansion is shared with `changebonds` via `changebond!`; the -# eigensolve + truncation + fidelity error are CBEDMRG-specific. Only `ϵs` (the two-site -# fidelity) is threaded back to drive convergence; the per-bond diagnostics (truncation, -# local Galerkin, and the expansion's own ϵ_select/ϵ_2site) are emitted via `@infov 4`. -function cbe_sweep!( - ψ, H, envs, alg::CBEDMRG, alg_eigsolve, dir::Val{D}, ϵs, timeroutput - ) where {D} - positions = D === :right ? (1:(length(ψ) - 1)) : (length(ψ):-1:2) - for pos in positions - local AC′, ϵ_trunc, v - # the two MPS tensors of the pre-update two-site state, kept by reference (the expansion - # builds new tensors rather than mutating these) for the DMRG2-style fidelity error - # below, so the dense two-site tensor is never formed - if D === :right - bra_left, bra_right = ψ.AC[pos], ψ.AR[pos + 1] - else - bra_left, bra_right = ψ.AL[pos - 1], ψ.AC[pos] - end - - # enrich the bond ahead of the optimization with directions orthogonal to ψ - # (the expansion logs its own ϵ_select/ϵ_2site diagnostics at @infov 4) - @timeit timeroutput "expand" changebond!(pos, dir, ψ, H, alg.expand, envs) - # the expanded neighbour tensor, kept for the fidelity contraction below - neighbor = D === :right ? ψ.AR[pos + 1] : ψ.AL[pos - 1] - - # single-site optimization in the enlarged space - @timeit timeroutput "AC_eigsolve" begin - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) - end + # 2. local update + @timeit timeroutput "AC_eigsolve" begin + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + end - # explicit truncated SVD: ϵ_trunc is the 2-norm of the discarded singular values, - # computed directly from the spectrum (no 1 - ‖C‖² cancellation). `v` is the overlap of - # the pre-update two-site state with the truncated, optimized one, contracted directly - # from the MPS factors in zipper order (peak intermediate is a single 3-leg tensor) so - # neither dense two-site tensor is formed; the truncated bond is internal, so the - # fidelity 1 - |v| -> 0 at convergence. - @timeit timeroutput "svd_trunc" begin - if D === :right - U, S, Vᴴ, ϵ_trunc = svd_trunc!(AC′; trunc = alg.trscheme, alg = Defaults.alg_svd()) - AL′ = U - C = normalize!(S * Vᴴ) - v = @plansor bra_left[1 2; 7] * conj(AL′[1 2; 5]) * conj(C[5; 6]) * - bra_right[7 4; 3] * conj(neighbor[6 4; 3]) - ψ.AC[pos] = (AL′, C) - else - U, S, AR′, ϵ_trunc = svd_trunc!(_transpose_tail(AC′); trunc = alg.trscheme, alg = Defaults.alg_svd()) - C = normalize!(U * S) - AR_mps = _transpose_front(AR′) - v = @plansor bra_left[1 2; 7] * conj(neighbor[1 2; 5]) * conj(C[5; 6]) * - bra_right[7 4; 3] * conj(AR_mps[6 4; 3]) - ψ.AC[pos] = (C, AR_mps) - end - end - fid = abs(1 - abs(v)) - ϵs[pos] = max(ϵs[pos], fid) - @infov 4 "CBEDMRG bond" pos fidelity = fid truncation = ϵ_trunc galerkin = calc_galerkin(pos, ψ, H, ψ, envs) - end - return ψ -end + # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`) + @timeit timeroutput "gauge" left_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true) + end -function find_groundstate!(ψ::AbstractFiniteMPS, H::FiniteMPOHamiltonian, alg::CBEDMRG, envs = environments(ψ, H, ψ)) - ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) - ϵ = maximum(ϵs) - log = IterLog("CBEDMRG") - timeroutput = TimerOutput("CBEDMRG") - alg.verbosity > 3 || disable_timer!(timeroutput) + # right-to-left + for pos in length(ψ):-1:2 + local AC′ + # convergence: pre-expansion single-site Galerkin error + ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs)) - LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs)) - for iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) + # 1. expand + isnothing(alg.alg_expand) || + @timeit timeroutput "expand" changebond!(pos, Val(:left), ψ, H, alg.alg_expand, envs) - zerovector!(ϵs) + # 2. local update + @timeit timeroutput "AC_eigsolve" begin + Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) + end - @timeit timeroutput "sweep" begin - cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:right), ϵs, timeroutput) - cbe_sweep!(ψ, H, envs, alg, alg_eigsolve, Val(:left), ϵs, timeroutput) + # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`) + @timeit timeroutput "gauge" right_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true) + end end - ϵ = maximum(ϵs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( @@ -345,6 +258,6 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm return ψ, envs, ϵ end -function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2, CBEDMRG}, envs...; kwargs...) +function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) return find_groundstate!(copy(ψ), H, alg, envs...; kwargs...) end diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index 05bb693d5..66bbb7bae 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -71,11 +71,11 @@ verbosity_conv = 1 # test logging ψ, envs, δ = find_groundstate( - ψ₀, H, CBEDMRG(; verbosity = verbosity_full, maxiter = 2, expand, trscheme) + ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trscheme) ) ψ, envs, δ = find_groundstate( - ψ, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 10, expand, trscheme), envs + ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_expand = expand, trscheme), envs ) v = variance(ψ, H) @@ -97,11 +97,11 @@ verbosity_conv = 1 # test logging ψ, envs, δ = find_groundstate( - ψ₀, H, CBEDMRG(; verbosity = verbosity_full, maxiter = 2, expand, trscheme) + ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trscheme) ) ψ, envs, δ = find_groundstate( - ψ, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 15, expand, trscheme), envs + ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trscheme), envs ) v = variance(ψ, H) @@ -121,7 +121,7 @@ verbosity_conv = 1 ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) expand = Exp(; trscheme = truncrank(D ÷ 2), warmstart = true, kw...) ψ, envs, δ = find_groundstate( - ψ₀, H, CBEDMRG(; verbosity = verbosity_conv, maxiter = 15, expand, trscheme = truncrank(D)) + ψ₀, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trscheme = truncrank(D)) ) @test sum(δ) ≈ 0 atol = 1.0e-3 @test variance(ψ, H) < 1.0e-2 From 3c5af8ed0dad8a38ab05130a043c25c0c06c57ff Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 14:38:35 -0400 Subject: [PATCH 07/20] more gauging utilities --- src/states/orthoview.jl | 79 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 73 insertions(+), 6 deletions(-) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index facc91a9c..072b99e0f 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -60,8 +60,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E} end for j in Iterators.reverse((i + 1):center) - v.parent.Cs[j], tmp = lq_compact!(_transpose_tail(v.parent.ACs[j]; copy = true); positive = true) - v.parent.ARs[j] = _transpose_front(tmp) + v.parent.Cs[j], v.parent.ARs[j] = right_gauge(v.parent.ACs[j]) if j != i + 1 # last AC not needed v.parent.ACs[j - 1] = _mul_tail(v.parent.ALs[j - 1], v.parent.Cs[j]) end @@ -76,7 +75,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E} end for j in center:i - v.parent.ALs[j], v.parent.Cs[j + 1] = qr_compact(v.parent.ACs[j]; positive = true) + v.parent.ALs[j], v.parent.Cs[j + 1] = left_gauge(v.parent.ACs[j]) if j != i # last AC not needed v.parent.ACs[j + 1] = _mul_front(v.parent.Cs[j + 1], v.parent.ARs[j + 1]) end @@ -89,10 +88,9 @@ end function Base.setindex!(v::CView{<:FiniteMPS}, vec, i::Int) if ismissing(v.parent.Cs[i + 1]) if !ismissing(v.parent.ALs[i]) - v.parent.Cs[i + 1], temp = lq_compact!(_transpose_tail(v.parent.AC[i + 1]; copy = true); positive = true) - v.parent.ARs[i + 1] = _transpose_front(temp) + v.parent.Cs[i + 1], v.parent.ARs[i + 1] = right_gauge(v.parent.AC[i + 1]) else - v.parent.ALs[i], v.parent.Cs[i + 1] = qr_compact(v.parent.AC[i]; positive = true) + v.parent.ALs[i], v.parent.Cs[i + 1] = left_gauge(v.parent.AC[i]) end end @@ -221,3 +219,72 @@ function Base.checkbounds( checkbounds(Bool, CView(first(psi.parent.data)), b) end end + +# Gauging routines +# ---------------- +@doc """ + left_gauge(AC, [alg]) -> AL, C + right_gauge(AC, [alg]) -> C, AR + +Factor an updated center MPS tensor `AC` into left- or right-canonical form, `AC ≈ AL * C` +(left, with `AL` left-isometric) or `AC ≈ C * AR` (right, with `AR` right-isometric), without +modifying `AC`. `right_gauge` handles the MPS leg permutation internally, so `AR` is returned in +standard MPS-tensor form. + +`alg` selects the factorization and defaults to a (positive) QR/LQ center-move that preserves +the virtual bond ([`Defaults.alg_orth`](@ref)). Passing a +[`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated +SVD that shrinks the bond, e.g. to cut an expanded bond back down in controlled bond expansion +[`DMRG`](@ref). + +!!! note + A `TruncatedAlgorithm` built from the default SVD placeholder (`Defaults.alg_svd()`) has no + registered `left_orth_alg`/`right_orth_alg`, so it is wrapped explicitly in + `LeftOrthViaSVD`/`RightOrthViaSVD`; the inner default SVD is resolved by `svd_trunc!`. +""" +left_gauge +@doc (@doc left_gauge) right_gauge + +left_gauge(AC) = left_gauge(AC, Defaults.alg_orth()) +left_gauge(AC, alg) = left_orth(AC; alg) +left_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) = + left_orth(AC; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg)) + +right_gauge(AC) = right_gauge(AC, Defaults.alg_orth()) +function right_gauge(AC, alg) + C, AR = right_orth(_transpose_tail(AC); alg) + return C, _transpose_front(AR) +end +function right_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) + C, AR = right_orth(_transpose_tail(AC); alg = MatrixAlgebraKit.RightOrthViaSVD(alg)) + return C, _transpose_front(AR) +end + +@doc """ + left_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ + right_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ + +Gauge an updated center tensor `AC` at site `pos` and install it into `ψ` in one step: factor +`AC` with [`left_gauge`](@ref) / [`right_gauge`](@ref) and write the canonical tensors back, +shifting the gauge center past `pos` (to the right for `left_gauge!`, to the left for +`right_gauge!`). `alg` is forwarded to [`left_gauge`](@ref) / [`right_gauge`](@ref) and hence may +be a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) to truncate the bond. + +By default the factors are installed as-is. Pass `normalize = true` to renormalize the bond +tensor, so `ψ` stays normalized after a local update that changed its norm. +""" +left_gauge! +@doc (@doc left_gauge!) right_gauge! + +function left_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false) + AL, C = left_gauge(AC, alg) + normalize && normalize!(C) + ψ.AC[pos] = (AL, C) + return ψ +end +function right_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false) + C, AR = right_gauge(AC, alg) + normalize && normalize!(C) + ψ.AC[pos] = (C, AR) + return ψ +end From c0be28a431f0e344faafdad9e3a42316ef30e95a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 14:38:43 -0400 Subject: [PATCH 08/20] fix test environment --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 562dd4339..b3c91b9ab 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,6 +12,7 @@ MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From d226f44976c8815fce9463f478a9923b8952cfbd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 14:39:38 -0400 Subject: [PATCH 09/20] format --- test/algorithms/groundstate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index 66bbb7bae..1338669ce 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -114,7 +114,7 @@ verbosity_conv = 1 end @testset "CBEDMRG warmstart $(nameof(Exp))" for (Exp, kw) in - ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) + ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) # warmstart seeds the expansion with the two-site gradient (alters the state); the # ground-state search should still converge to the correct low-variance state Random.seed!(2025) From c7e2c0f6675e81aaa238a97d846c9c5291458278 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 20:23:35 -0400 Subject: [PATCH 10/20] docstring fix --- src/states/orthoview.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 072b99e0f..bb2338340 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -231,16 +231,8 @@ Factor an updated center MPS tensor `AC` into left- or right-canonical form, `AC modifying `AC`. `right_gauge` handles the MPS leg permutation internally, so `AR` is returned in standard MPS-tensor form. -`alg` selects the factorization and defaults to a (positive) QR/LQ center-move that preserves -the virtual bond ([`Defaults.alg_orth`](@ref)). Passing a -[`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated -SVD that shrinks the bond, e.g. to cut an expanded bond back down in controlled bond expansion -[`DMRG`](@ref). - -!!! note - A `TruncatedAlgorithm` built from the default SVD placeholder (`Defaults.alg_svd()`) has no - registered `left_orth_alg`/`right_orth_alg`, so it is wrapped explicitly in - `LeftOrthViaSVD`/`RightOrthViaSVD`; the inner default SVD is resolved by `svd_trunc!`. +`alg` selects the factorization and defaults to a (positive) QR/LQ center-move that preserves the virtual bond. +Passing a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated SVD may shrink the bond. """ left_gauge @doc (@doc left_gauge) right_gauge From a4ed9fac4eaa0d8312e88d6656e9a86e4ebb722d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 20:23:48 -0400 Subject: [PATCH 11/20] CBE-TDVP --- src/algorithms/timestep/tdvp.jl | 70 +++++++++++++++++++++++++++++---- test/algorithms/timestep.jl | 46 ++++++++++++++++++++++ 2 files changed, 109 insertions(+), 7 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index a6fc56c2f..c835163bd 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -3,6 +3,19 @@ $(TYPEDEF) Single site MPS time-evolution algorithm based on the Time-Dependent Variational Principle. +For finite MPS, setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref), +[`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to the current state +ahead of each local integration, recovering Controlled Bond Expansion (CBE) TDVP and lifting the +fixed-bond limitation of plain single-site TDVP. A truncating `trscheme` is then required to cut +the enlarged bond back down (selecting the truncated-SVD gauge). The expansion must be +state-preserving, so leave the expander's `warmstart = false` (the default); `warmstart = true` +injects a gradient and corrupts the evolution. + +!!! note + The expansion renormalizes the bond tensor, which is a no-op for the norm-1 states of + real-time evolution but renormalizes imaginary-time evolution at each expansion. CBE is only + available for finite MPS. + ## Fields $(TYPEDFIELDS) @@ -11,18 +24,38 @@ $(TYPEDFIELDS) * [Haegeman et al. Phys. Rev. Lett. 107 (2011)](@cite haegeman2011) """ -@kwdef struct TDVP{A, F} <: Algorithm +struct TDVP{A, E, G, F} <: Algorithm "algorithm used in the exponential solvers" - integrator::A = Defaults.alg_expsolve() + integrator::A "tolerance for gauging algorithm" - tolgauge::Float64 = Defaults.tolgauge + tolgauge::Float64 "maximal amount of iterations for gauging algorithm" - gaugemaxiter::Int = Defaults.maxiter + gaugemaxiter::Int + + "algorithm used to expand the bond ahead of each local update, or `nothing` for none (finite CBE-TDVP)" + alg_expand::E + + "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD" + alg_gauge::G "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" - finalize::F = Defaults._finalize + finalize::F +end +function TDVP(; + integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, + gaugemaxiter = Defaults.maxiter, finalize = Defaults._finalize, + alg_expand = nothing, trscheme = notrunc(), + alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth() + ) + # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD + alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth : + MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme) + if !isnothing(alg_expand) && !_truncates(alg_gauge) + @warn "TDVP with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep." + end + return TDVP(integrator, tolgauge, gaugemaxiter, alg_expand, alg_gauge, finalize) end function timestep( @@ -92,9 +125,21 @@ function timestep!( # sweep left to right for i in 1:(length(ψ) - 1) + # 1. optionally expand the bond ahead of the local update (CBE) + isnothing(alg.alg_expand) || + LoggingExtras.withlevel(; alg.alg_expand.verbosity) do + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs) + end + + # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) - ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) + AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) + # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the + # enlarged bond back down) and move the center to i+1 + left_gauge!(ψ, i, AC, alg.alg_gauge) + + # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs) ψ.C[i] = integrate( Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator; @@ -108,12 +153,23 @@ function timestep!( # sweep right to left for i in length(ψ):-1:2 + # 1. optionally expand the bond ahead of the local update (CBE) + isnothing(alg.alg_expand) || + LoggingExtras.withlevel(; alg.alg_expand.verbosity) do + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs) + end + + # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) - ψ.AC[i] = integrate( + AC = integrate( Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator; imaginary_evolution ) + # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 + right_gauge!(ψ, i, AC, alg.alg_gauge) + + # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs) ψ.C[i - 1] = integrate( Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator; diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index b6a3eb8dc..224507a44 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -9,6 +9,8 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ +using LinearAlgebra: norm +using Random verbosity_full = 5 verbosity_conv = 1 @@ -91,6 +93,50 @@ verbosity_conv = 1 end end +@testset "Finite CBE-TDVP" verbose = true begin + L = 10 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + Dstart, Dcap, dt = 2, 16, 0.05 + + # controlled bond expansion lets single-site TDVP grow the bond; the evolution should stay + # unitary (norm-preserving) and energy-conserving while tracking the bond-adaptive TDVP2 + # reference better than fixed-bond single-site TDVP + @testset "$(nameof(Exp))" for (Exp, kw) in + ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) + Random.seed!(4) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) + alg = TDVP(; alg_expand = Exp(; trscheme = truncrank(Dstart), kw...), trscheme = truncrank(Dcap)) + E₀ = real(expectation_value(ψ₀, H)) + + ref, cbe, plain = ψ₀, ψ₀, ψ₀ + for _ in 1:6 + ref, = timestep(ref, H, 0.0, dt, TDVP2(; trscheme = truncrank(Dcap))) + cbe, = timestep(cbe, H, 0.0, dt, alg) + plain, = timestep(plain, H, 0.0, dt, TDVP()) # stuck at Dstart + end + + @test norm(cbe) ≈ 1 atol = 1.0e-8 + @test real(expectation_value(cbe, H)) ≈ E₀ atol = 1.0e-2 + @test dim(left_virtualspace(cbe, L ÷ 2)) > Dstart + @test dim(left_virtualspace(plain, L ÷ 2)) == Dstart + @test abs(dot(ref, cbe)) > abs(dot(ref, plain)) + end + + @testset "imaginary-time lowers energy" begin + Random.seed!(5) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) + alg = TDVP(; alg_expand = OptimalExpand(; trscheme = truncrank(Dstart)), trscheme = truncrank(Dcap)) + E₀ = real(expectation_value(ψ₀, H)) + ψ = ψ₀ + for _ in 1:8 + ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true) + normalize!(ψ) + end + @test real(expectation_value(ψ, H)) < E₀ + @test dim(left_virtualspace(ψ, L ÷ 2)) > Dstart + end +end + @testset "time_evolve" verbose = true begin t_span = 0:0.1:0.1 algs = [TDVP(), TDVP2(; trscheme = truncrank(10))] From f1f7c46df3f7aebaafc7074b4ac50602a7c228be Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 20:24:15 -0400 Subject: [PATCH 12/20] project_complement --- src/algorithms/changebonds/sketchedexpand.jl | 4 ++-- src/algorithms/toolbox.jl | 3 +-- src/utility/utility.jl | 8 ++++++++ 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index 3f33f58ad..4ab04cfd1 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -107,7 +107,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk # by a narrow QR into a proper randomized range-finder basis (no full complement basis is formed) Vℓ = sketch_space(compL, alg) Ω = randisometry(scalartype(left), codomain(AL) ← Vℓ) - Q, _ = qr_compact!(Ω - AL * (AL' * Ω)) + Q, _ = qr_compact!(project_complement!(Ω, AL)) # fold the sketch into the left environment (bra-leg becomes the sketch leg), then apply the # existing single-site effective Hamiltonian on site+1: this reconstructs `Q† AC2_projection` @@ -164,7 +164,7 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske Y = Hac * left # project onto the left complement with (I - AL AL'): the SVD left-vectors lie in the complement - B = Y - left * (left' * Y) + B = project_complement!(Y, left) nrm = norm(B) trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compL)) U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index a319baede..4b9df6e58 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -56,8 +56,7 @@ Concretely, this is the overlap of the current state with the single-site deriva function calc_galerkin(pos::Int, below, operator, above, envs) AC´ = AC_projection(pos, below, operator, above, envs) normalize!(AC´) - out = mul!(AC´, below.AL[pos], below.AL[pos]' * AC´, -1, +1) - return norm(out) + return norm(project_complement!(AC´, below.AL[pos])) end function calc_galerkin(pos::CartesianIndex{2}, below, operator, above, envs) row, col = Tuple(pos) diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 670439523..6c925a151 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -11,6 +11,14 @@ end _mul_front(C, A) = mul_front(C, A) # _transpose_front(C * _transpose_tail(A)) _mul_tail(A, C) = mul_tail(A, C) # A * C +""" + project_complement!(Y, X) -> Y + +In-place projection of `Y` onto the orthogonal complement of the range of the left-isometry `X` +(`X' X = I`): `Y ← (I - X X') Y = Y - X (X' Y)`. `Y` is overwritten and returned. +""" +project_complement!(Y, X) = mul!(Y, X, X' * Y, -1, +1) + function _similar_tail(A::AbstractTensorMap) cod = _firstspace(A) dom = ⊗(dual(_lastspace(A)), dual.(space.(Ref(A), reverse(2:(numind(A) - 1))))...) From 97e00b0f5cdf247c4fc4b8738c2a7eb9dbbb22d7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 21:20:24 -0400 Subject: [PATCH 13/20] more updates and normalization --- src/algorithms/changebonds/optimalexpand.jl | 10 ++++--- src/algorithms/changebonds/randexpand.jl | 10 ++++--- src/algorithms/changebonds/sketchedexpand.jl | 16 ++++++----- src/algorithms/timestep/tdvp.jl | 19 +++++++------ src/utility/utility.jl | 12 +++++++++ test/algorithms/changebonds.jl | 28 ++++++++++++++++++++ test/algorithms/timestep.jl | 26 +++++++++++++++--- 7 files changed, 95 insertions(+), 26 deletions(-) diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 3d2d15db1..43f713e2c 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -99,7 +99,7 @@ end # Finite system # ------------- -function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs) +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs; normalize::Bool = true) bond = site left = ψ.AC[site] right = ψ.AR[site + 1] @@ -128,11 +128,12 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) - ψ.AC[site] = (nal, normalize!(nc)) + normalize && normalize!(nc) + ψ.AC[site] = (nal, nc) ψ.AC[site + 1] = (nc, nar) return ψ end -function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs) +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs; normalize::Bool = true) bond = site - 1 left = ψ.AL[site - 1] right = ψ.AC[site] @@ -162,7 +163,8 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Opt end AL_exp = catdomain(left, Q) - ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) + normalize && normalize!(nc) + ψ.AC[site] = (nc, _transpose_front(Qr)) ψ.AC[site - 1] = (AL_exp, nc) return ψ end diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index ad84a66af..6eacafbe8 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -82,7 +82,7 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) return normalize!(ψ) end -function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs) +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs; normalize::Bool = true) bond = site left = ψ.AC[site] right = ψ.AR[site + 1] @@ -105,11 +105,12 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) - ψ.AC[site] = (nal, normalize!(nc)) + normalize && normalize!(nc) + ψ.AC[site] = (nal, nc) ψ.AC[site + 1] = (nc, nar) return ψ end -function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs) +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::RandExpand, envs; normalize::Bool = true) bond = site - 1 left = ψ.AL[site - 1] right = ψ.AC[site] @@ -133,7 +134,8 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ran nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) AL_exp = catdomain(left, Q) - ψ.AC[site] = (normalize!(nc), _transpose_front(Qr)) + normalize && normalize!(nc) + ψ.AC[site] = (nc, _transpose_front(Qr)) ψ.AC[site - 1] = (AL_exp, nc) return ψ end diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index 4ab04cfd1..284296c20 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -60,7 +60,7 @@ by `V`). Reuses [`sample_space`](@ref) so the result is sector-correct for grade function sketch_space(V, alg::SketchedExpand) Vk = sample_space(V, alg.trscheme) alg.oversampling == 0 && return Vk - Vp = sample_space(V, MatrixAlgebraKit.truncrank(alg.oversampling)) + Vp = typeof(V)(c => min(dim(V, c), alg.oversampling) for c in sectors(V)) return infimum(V, Vk ⊕ Vp) end @@ -92,7 +92,7 @@ end # complement automatically). The kept rank is capped at the complement dimension: the projected # object has only that many genuine singular values, and at edge bonds the complement may be empty, # where an uncapped `svd_trunc` would otherwise pad the selection with non-orthogonal noise. -function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs) +function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true) left = ψ.AC[site] right = ψ.AR[site + 1] AL, _ = left_orth(left) # local left-isometric form (range(AL) = range(left)) @@ -118,7 +118,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk # project onto the right complement with (I - AR' AR): the SVD right-vectors lie in the complement Ytt = _transpose_tail(Y) - B = Ytt - (Ytt * ARtt') * ARtt + B = project_complement_right!(Ytt, ARtt) nrm = norm(B) trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compR)) U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) @@ -136,11 +136,12 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) - ψ.AC[site] = (nal, normalize!(nc)) + normalize && normalize!(nc) + ψ.AC[site] = (nal, nc) ψ.AC[site + 1] = (nc, nar) return ψ end -function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs) +function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true) left = ψ.AL[site - 1] # left-isometric, left of center (valid) right = ψ.AC[site] _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc()) # local right-iso @@ -154,7 +155,7 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske # orthonormal rows recovered by a narrow LQ (the randomized range-finder basis) Vℓ = sketch_space(compR, alg) Ω = adjoint(randisometry(scalartype(right), domain(ARtt) ← Vℓ)) - _, Qr_o = lq_compact!(Ω - (Ω * ARtt') * ARtt) + _, Qr_o = lq_compact!(project_complement_right!(Ω, ARtt)) Qr = _transpose_front(Qr_o) # fold the sketch into the right environment, then apply the single-site effective @@ -183,7 +184,8 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske end AL_exp = catdomain(left, Q) - ψ.AC[site] = (normalize!(nc), _transpose_front(Qr2)) + normalize && normalize!(nc) + ψ.AC[site] = (nc, _transpose_front(Qr2)) ψ.AC[site - 1] = (AL_exp, nc) return ψ end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index c835163bd..396a68482 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -12,8 +12,9 @@ state-preserving, so leave the expander's `warmstart = false` (the default); `wa injects a gradient and corrupts the evolution. !!! note - The expansion renormalizes the bond tensor, which is a no-op for the norm-1 states of - real-time evolution but renormalizes imaginary-time evolution at each expansion. CBE is only + Real-time evolution preserves the norm: neither the bond expansion nor the truncation + renormalizes, so the state norm reflects the accumulated truncation error. Imaginary-time + evolution instead renormalizes at every step, like a ground-state search. CBE is only available for finite MPS. ## Fields @@ -128,7 +129,7 @@ function timestep!( # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || LoggingExtras.withlevel(; alg.alg_expand.verbosity) do - changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs) + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) end # 2. evolve the (possibly expanded) center tensor forward @@ -136,8 +137,9 @@ function timestep!( AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the - # enlarged bond back down) and move the center to i+1 - left_gauge!(ψ, i, AC, alg.alg_gauge) + # enlarged bond back down) and move the center to i+1. Real-time evolution preserves + # the norm; imaginary-time evolution renormalizes. + left_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs) @@ -156,7 +158,7 @@ function timestep!( # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || LoggingExtras.withlevel(; alg.alg_expand.verbosity) do - changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs) + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) end # 2. evolve the (possibly expanded) center tensor forward @@ -166,8 +168,9 @@ function timestep!( imaginary_evolution ) - # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 - right_gauge!(ψ, i, AC, alg.alg_gauge) + # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (real-time preserves the + # norm; imaginary-time renormalizes) + right_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs) diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 6c925a151..83f67aadc 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -16,9 +16,21 @@ _mul_tail(A, C) = mul_tail(A, C) # A * C In-place projection of `Y` onto the orthogonal complement of the range of the left-isometry `X` (`X' X = I`): `Y ← (I - X X') Y = Y - X (X' Y)`. `Y` is overwritten and returned. + +See also [`project_complement_right!`](@ref). """ project_complement!(Y, X) = mul!(Y, X, X' * Y, -1, +1) +""" + project_complement_right!(Y, X) -> Y + +In-place projection of `Y` onto the orthogonal complement of the co-range of the right-isometry +`X` (`X X' = I`): `Y ← Y (I - X' X) = Y - (Y X') X`. `Y` is overwritten and returned. + +See also [`project_complement!`](@ref). +""" +project_complement_right!(Y, X) = mul!(Y, Y * X', X, -1, +1) + function _similar_tail(A::AbstractTensorMap) cod = _firstspace(A) dom = ⊗(dual(_lastspace(A)), dual.(space.(Ref(A), reverse(2:(numind(A) - 1))))...) diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl index 958df806d..a085f587e 100644 --- a/test/algorithms/changebonds.jl +++ b/test/algorithms/changebonds.jl @@ -9,6 +9,7 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ +using Random spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))] @@ -105,6 +106,33 @@ end @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5)) end +# density-matrix-style MPS: each site carries two physical legs (ket ⊗ bra). The operator-free +# bond-change algorithms (`RandExpand` expansion, `SvdCut` truncation) must handle the extra +# physical leg. Operator-based expanders (`OptimalExpand`/`SketchedExpand`) are not covered here +# because `FiniteMPOHamiltonian`/`FiniteMPO` only support a single physical leg per site. +@testset "Density-matrix FiniteMPS $(spacetype(pcomp))" for (pcomp, Dspace) in [ + (ℙ^2 ⊗ (ℙ^2)', ℙ^6), + (Rep[SU₂](1 // 2 => 1) ⊗ Rep[SU₂](1 // 2 => 1)', Rep[SU₂](0 => 4, 1 => 3)), + ] + Random.seed!(2468) + L = 8 + maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 2:length(ψ)) + + ψ = FiniteMPS(rand, ComplexF64, fill(pcomp, L), Dspace) + @test numind(ψ.AC[L ÷ 2]) == 4 # two physical legs + two virtual legs + + # RandExpand grows the bond while preserving the state (norm-preserving expansion) + ψ_re = changebonds(ψ, RandExpand(; trscheme = truncrank(dim(Dspace) * 2))) + @test numind(ψ_re.AC[L ÷ 2]) == 4 + @test abs(dot(ψ, ψ_re)) ≈ 1 atol = 1.0e-8 + @test maxbond(ψ_re) > maxbond(ψ) + + # SvdCut truncates the enlarged bond back down, leaving a normalized state + ψ_tr = changebonds(ψ_re, SvdCut(; trscheme = truncrank(dim(Dspace)))) + @test maxbond(ψ_tr) < maxbond(ψ_re) + @test abs(dot(ψ_tr, ψ_tr)) ≈ 1 atol = 1.0e-8 +end + @testset "MultilineMPS $(spacetype(pspace))" for (pspace, Dspace) in spacelist o = rand(ComplexF64, pspace * pspace, pspace * pspace) mpo = MultilineMPO(o) diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index 224507a44..a75cb1cb7 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -115,13 +115,34 @@ end plain, = timestep(plain, H, 0.0, dt, TDVP()) # stuck at Dstart end - @test norm(cbe) ≈ 1 atol = 1.0e-8 + @test norm(cbe) ≈ 1 atol = 1.0e-6 @test real(expectation_value(cbe, H)) ≈ E₀ atol = 1.0e-2 @test dim(left_virtualspace(cbe, L ÷ 2)) > Dstart @test dim(left_virtualspace(plain, L ÷ 2)) == Dstart @test abs(dot(ref, cbe)) > abs(dot(ref, plain)) end + # the bond truncation must preserve the norm for real-time evolution (the norm reflects the + # discarded weight) and only renormalize for imaginary-time evolution + @testset "norm handling" begin + Random.seed!(6) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) + # a deliberately lossy cap so the truncation discards weight every step + lossy = TDVP(; alg_expand = OptimalExpand(; trscheme = truncrank(2)), trscheme = truncrank(2)) + + ψrt = ψ₀ + for _ in 1:12 + ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy) # real time + end + @test norm(ψrt) < 1 - 1.0e-3 # truncation loss is not renormalized away + + ψit = ψ₀ + for _ in 1:12 + ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true) # no external normalize! + end + @test norm(ψit) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + @testset "imaginary-time lowers energy" begin Random.seed!(5) ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) @@ -129,8 +150,7 @@ end E₀ = real(expectation_value(ψ₀, H)) ψ = ψ₀ for _ in 1:8 - ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true) - normalize!(ψ) + ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true) # gauge renormalizes end @test real(expectation_value(ψ, H)) < E₀ @test dim(left_virtualspace(ψ, L ÷ 2)) > Dstart From 65c921b4b6297e37a7bf9b5ff31b53786c8e39ab Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 29 Jun 2026 22:19:25 -0400 Subject: [PATCH 14/20] more improvements --- src/algorithms/changebonds/optimalexpand.jl | 4 +- src/algorithms/changebonds/randexpand.jl | 2 +- src/algorithms/changebonds/sketchedexpand.jl | 135 +++++++------------ 3 files changed, 55 insertions(+), 86 deletions(-) diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 43f713e2c..5a44d6326 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -120,11 +120,11 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op if alg.warmstart # seed the new left directions with the (physically scaled) gradient instead of a zero # block, warm-starting the subsequent optimization (alters the state) - nal, nc = qr_compact!(catdomain(left, gradnorm * (NL * U * S))) + nal, nc = left_gauge(catdomain(left, gradnorm * (NL * U * S))) else # embed `left` into the enlarged domain (zero weight in the new directions) nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) - nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 6eacafbe8..f7cf68699 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -102,7 +102,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra ar_re = Vᴴ * NR # embed `left` into the enlarged domain (zero weight in the new directions) nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) - nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) normalize && normalize!(nc) diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index 284296c20..42ec28642 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -1,47 +1,33 @@ """ $(TYPEDEF) -An algorithm that expands the bond dimension by selecting the dominant directions of the -projected two-site update orthogonal to the current state, just like [`OptimalExpand`](@ref), -but at single-site cost using a randomized range-finder (the "shrewd selection" of Controlled -Bond Expansion, Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)). +An algorithm that expands the bond dimension like [`OptimalExpand`](@ref) — selecting the +dominant directions of the projected two-site update orthogonal to the current state — but at +single-site cost using the randomized "shrewd selection" of Controlled Bond Expansion +(Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)). A random sketch of the orthogonal complement +is folded into the effective environment, collapsing the large bond before the two-site update is +ever formed, and the dominant directions are read off a small singular value decomposition. -Rather than forming the full projected two-site update (which requires applying the two-site -effective Hamiltonian), a random sketch of width `rank + oversampling` of the orthogonal -complement is folded into the effective *environment* first, collapsing the large bond before -the two-site contraction is ever materialized. The dominant directions are then recovered from -a small singular value decomposition of the sketched gradient. +The `warmstart` and state-preserving behaviour match [`OptimalExpand`](@ref). !!! note - This strategy is only defined for `FiniteMPS` (through [`changebond!`](@ref)), so that it - can be used both standalone and as the `alg_expand` strategy of [`DMRG`](@ref). - -Like [`OptimalExpand`](@ref), the expansion is by default state-preserving (the added -directions are connected through a zero block, as required for e.g. TDVP). When -`warmstart = true`, the new directions are instead seeded with the (physically scaled) sketched -two-site gradient, warm-starting the subsequent single-site optimization in ground-state search -(e.g. as the `alg_expand` strategy of [`DMRG`](@ref)); this alters the state, and the injected -amplitude scales with the gradient norm so that it vanishes automatically at convergence. - -!!! note - `ϵ_2site` is a randomized estimate of the two-site complement-gradient norm (a diagnostic - only), not the exact ratio `‖g2‖ / ‖AC2‖` reported by [`OptimalExpand`](@ref). The folded - application routes the `MPOHamiltonian` tensors through the generic transfer/one-site - machinery and therefore does not exploit `JordanMPO` sparsity in the MPO bond dimension. + Only defined for `FiniteMPS` (through [`changebond!`](@ref)), so it can be used standalone or + as the `alg_expand` strategy of [`DMRG`](@ref). The reported `ϵ_2site` is a randomized + estimate, and the folded application does not exploit `JordanMPO` sparsity. ## Fields $(TYPEDFIELDS) """ @kwdef struct SketchedExpand{S} <: Algorithm - "algorithm used for the (small) singular value decomposition" + "algorithm used for the singular value decomposition" alg_svd::S = Defaults.alg_svd() "algorithm used for truncating the expanded space" trscheme::TruncationStrategy "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)" - oversampling::Int = 5 + oversampling::Int = 0 "whether to seed the new directions with the sketched two-site gradient (warm start, alters the state) instead of a zero block" warmstart::Bool = false @@ -51,26 +37,24 @@ $(TYPEDFIELDS) end """ - sketch_space(V, alg::SketchedExpand) + sketch_space(V, alg::SketchedExpand) -> Vℓ, Vk -Determine the space of the random sketch drawn from the complement space `V`: the target -subspace selected by `alg.trscheme`, enlarged by `alg.oversampling` extra directions (capped -by `V`). Reuses [`sample_space`](@ref) so the result is sector-correct for graded spaces. +The random-sketch space `Vℓ` drawn from the complement `V`, together with its oversampling-free +target `Vk` (selected by `alg.trscheme`). `Vℓ` enlarges `Vk` by `alg.oversampling` extra +directions (capped by `V`); the selection is truncated back to `Vk`. """ function sketch_space(V, alg::SketchedExpand) Vk = sample_space(V, alg.trscheme) - alg.oversampling == 0 && return Vk + alg.oversampling == 0 && return Vk, Vk Vp = typeof(V)(c => min(dim(V, c), alg.oversampling) for c in sectors(V)) - return infimum(V, Vk ⊕ Vp) + return infimum(V, Vk ⊕ Vp), Vk end """ complement_space(Vfull::ElementarySpace, image::ElementarySpace) -The orthogonal complement `Vfull ⊖ image` of the range `image` inside the fused space `Vfull`, -computed sector-wise from the spaces alone (no factorization). For an isometry `A : Vfull ← image` -this is the space carried by `left_null(A)`; here `Vfull` is the fused codomain/domain of the -isometric MPS tensor and `image` its small virtual leg. Used to size and cap the random sketch. +The orthogonal complement `Vfull ⊖ image`, computed sector-wise from the spaces alone (no +factorization). Used to size and cap the random sketch. """ function complement_space(Vfull::ElementarySpace, image::ElementarySpace) @assert !isdual(Vfull) @@ -81,58 +65,46 @@ end # Finite system # ------------- -# The "shrewd selection" needs the orthogonal complement of the current state on either side of -# the bond. Rather than materializing an explicit complement basis with `left_null`/`right_null`, -# we use the complement projector of the *isometric* MPS tensor `A`: `I - A A'` (left-iso) and -# `I - A' A` (right-iso). The isometric tensor is obtained from the center tensor by a local -# `left_orth`/`right_orth!`, which leaves the MPS gauge (and any environments a caller maintains -# incrementally, e.g. CBE `DMRG`) untouched. On the sketch side this turns the range-finder into a -# narrow random draw projected into the complement; on the selection side the gradient is projected -# and its dominant directions read off directly from the SVD (its singular vectors land in the -# complement automatically). The kept rank is capped at the complement dimension: the projected -# object has only that many genuine singular values, and at edge bonds the complement may be empty, -# where an uncapped `svd_trunc` would otherwise pad the selection with non-orthogonal noise. +# Rather than forming the two-site update, a random sketch of the orthogonal complement is folded +# into the effective environment, and the dominant directions are read off a small SVD of the +# sketched gradient. The complement projectors act with the isometric MPS tensors directly, which +# leaves the MPS gauge (and any incrementally-maintained environments) untouched. function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true) left = ψ.AC[site] right = ψ.AR[site + 1] - AL, _ = left_orth(left) # local left-isometric form (range(AL) = range(left)) - ARtt = _transpose_tail(right) # right-isometric: ARtt * ARtt' = I + AL, _ = left_gauge(left) # local left-isometric form + ARtt = _transpose_tail(right) # AR is already right-isometric # nothing to add when either complement is empty (e.g. edge bonds) compL = complement_space(fuse(codomain(AL)), only(domain(AL))) compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) (dim(compL) == 0 || dim(compR) == 0) && return ψ - # sketch the left complement: Q = (I - AL AL') Ω lies in the complement of `left`, re-orthonormalized - # by a narrow QR into a proper randomized range-finder basis (no full complement basis is formed) - Vℓ = sketch_space(compL, alg) + # random sketch of the left complement, folded into the left environment + Vℓ, Vk = sketch_space(compL, alg) Ω = randisometry(scalartype(left), codomain(AL) ← Vℓ) Q, _ = qr_compact!(project_complement!(Ω, AL)) - - # fold the sketch into the left environment (bra-leg becomes the sketch leg), then apply the - # existing single-site effective Hamiltonian on site+1: this reconstructs `Q† AC2_projection` - # without ever forming the full two-site update GL = leftenv(envs, site, ψ) * TransferMatrix(left, H[site], Q) Hac = MPO_AC_Hamiltonian(GL, H[site + 1], rightenv(envs, site + 1, ψ)) Y = Hac * right - # project onto the right complement with (I - AR' AR): the SVD right-vectors lie in the complement - Ytt = _transpose_tail(Y) - B = project_complement_right!(Ytt, ARtt) - nrm = norm(B) - trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compR)) + # select the dominant directions in the right complement, dropping the oversampling padding + B = project_complement_right!(_transpose_tail(Y), ARtt) + gradnorm = norm(B) + trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(Vk)) U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) - @infov 4 "bond expansion (sketched)" site dir = :right ϵ_select ϵ_2site = nrm + @infov 4 "bond expansion (sketched)" site dir = :right ϵ_select ϵ_2site = gradnorm - # optimal vectors at site+1 (identical bookkeeping to OptimalExpand) + # optimal vectors at site+1 ar_re = Vᴴ if alg.warmstart - # seed the new left directions with the (physically scaled) sketched gradient instead of - # a zero block; `Q * U` maps the sketch-basis left vectors back to the full space - nal, nc = qr_compact!(catdomain(left, nrm * (Q * U * S))) + # seed the new left directions with the (physically scaled) gradient instead of a zero + # block, warm-starting the subsequent optimization (alters the state) + nal, nc = left_gauge(catdomain(left, gradnorm * (Q * U * S))) else + # embed `left` into the enlarged domain (zero weight in the new directions) nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) - nal, nc = qr_compact!(absorb!(zerovector!(similar(left, nal_space)), left)) + nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) @@ -142,43 +114,40 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk return ψ end function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true) - left = ψ.AL[site - 1] # left-isometric, left of center (valid) + left = ψ.AL[site - 1] right = ψ.AC[site] - _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc()) # local right-iso + _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc()) # local right-isometric form # nothing to add when either complement is empty (e.g. edge bonds) compL = complement_space(fuse(codomain(left)), only(domain(left))) compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) (dim(compL) == 0 || dim(compR) == 0) && return ψ - # sketch the right complement: Qr = (I - AR' AR) Ω lies in the complement of `right`, with - # orthonormal rows recovered by a narrow LQ (the randomized range-finder basis) - Vℓ = sketch_space(compR, alg) + # random sketch of the right complement, folded into the right environment + Vℓ, Vk = sketch_space(compR, alg) Ω = adjoint(randisometry(scalartype(right), domain(ARtt) ← Vℓ)) _, Qr_o = lq_compact!(project_complement_right!(Ω, ARtt)) Qr = _transpose_front(Qr_o) - - # fold the sketch into the right environment, then apply the single-site effective - # Hamiltonian on site-1: reconstructs `AC2_projection NR†` (sketched) at single-site cost GR = TransferMatrix(right, H[site], Qr) * rightenv(envs, site, ψ) Hac = MPO_AC_Hamiltonian(leftenv(envs, site - 1, ψ), H[site - 1], GR) Y = Hac * left - # project onto the left complement with (I - AL AL'): the SVD left-vectors lie in the complement + # select the dominant directions in the left complement, dropping the oversampling padding B = project_complement!(Y, left) - nrm = norm(B) - trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(compL)) + gradnorm = norm(B) + trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(Vk)) U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) - @infov 4 "bond expansion (sketched)" site dir = :left ϵ_select ϵ_2site = nrm + @infov 4 "bond expansion (sketched)" site dir = :left ϵ_select ϵ_2site = gradnorm - # optimal vectors at site-1 (identical bookkeeping to OptimalExpand) + # optimal vectors at site-1 Q = U right_tail = _transpose_tail(right) if alg.warmstart - # seed the new right directions with the (physically scaled) sketched gradient instead of - # a zero block; `Vᴴ * Qr_o` maps the sketch-basis right vectors back to the full space - nc, Qr2 = lq_compact!(catcodomain(right_tail, nrm * (S * Vᴴ * Qr_o))) + # seed the new right directions with the (physically scaled) gradient instead of a zero + # block, warm-starting the subsequent optimization (alters the state) + nc, Qr2 = lq_compact!(catcodomain(right_tail, gradnorm * (S * Vᴴ * Qr_o))) else + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) end From 102c93337c5829eb23cd0a8090185f43ba8554d7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 08:18:55 -0400 Subject: [PATCH 15/20] some more clean up --- src/algorithms/changebonds/sketchedexpand.jl | 30 ++++++++++---------- src/algorithms/groundstate/dmrg.jl | 2 +- test/algorithms/timestep.jl | 2 +- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index 42ec28642..a4b870ff0 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -88,22 +88,22 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk Hac = MPO_AC_Hamiltonian(GL, H[site + 1], rightenv(envs, site + 1, ψ)) Y = Hac * right - # select the dominant directions in the right complement, dropping the oversampling padding + # orthonormalize the sketched right complement, truncating away the oversampling padding B = project_complement_right!(_transpose_tail(Y), ARtt) - gradnorm = norm(B) - trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(Vk)) - U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) - @infov 4 "bond expansion (sketched)" site dir = :right ϵ_select ϵ_2site = gradnorm + if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices + L, ar_re = right_orth!(B) + else + L, ar_re = right_orth!(B; trunc = truncspace(Vk)) + end # optimal vectors at site+1 - ar_re = Vᴴ if alg.warmstart # seed the new left directions with the (physically scaled) gradient instead of a zero # block, warm-starting the subsequent optimization (alters the state) - nal, nc = left_gauge(catdomain(left, gradnorm * (Q * U * S))) + nal, nc = left_gauge(catdomain(left, Q * L)) else # embed `left` into the enlarged domain (zero weight in the new directions) - nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(ar_re, 1)) nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) end nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) @@ -132,20 +132,20 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske Hac = MPO_AC_Hamiltonian(leftenv(envs, site - 1, ψ), H[site - 1], GR) Y = Hac * left - # select the dominant directions in the left complement, dropping the oversampling padding + # orthonormalize the sketched left complement, truncating away the oversampling padding B = project_complement!(Y, left) - gradnorm = norm(B) - trunc = alg.trscheme & MatrixAlgebraKit.truncrank(dim(Vk)) - U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(B); trunc, alg = alg.alg_svd) - @infov 4 "bond expansion (sketched)" site dir = :left ϵ_select ϵ_2site = gradnorm + if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices + Q, R = left_orth!(B) + else + Q, R = left_orth!(B; trunc = truncspace(Vk)) + end # optimal vectors at site-1 - Q = U right_tail = _transpose_tail(right) if alg.warmstart # seed the new right directions with the (physically scaled) gradient instead of a zero # block, warm-starting the subsequent optimization (alters the state) - nc, Qr2 = lq_compact!(catcodomain(right_tail, gradnorm * (S * Vᴴ * Qr_o))) + nc, Qr2 = lq_compact!(catcodomain(right_tail, R * Qr_o)) else # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index bf8c8451d..61f3097a0 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -6,7 +6,7 @@ _truncates(::MatrixAlgebraKit.TruncatedAlgorithm) = true """ $(TYPEDEF) -Single-site DMRG algorithm for finding the dominant eigenvector. +Density Matrix Renormalization Group algorithm for finding the dominant eigenvector. Each site update is, in order: (1) an optional bond expansion (`alg_expand`), (2) a single-site eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and a diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index a75cb1cb7..a846b2c5c 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -102,7 +102,7 @@ end # unitary (norm-preserving) and energy-conserving while tracking the bond-adaptive TDVP2 # reference better than fixed-bond single-site TDVP @testset "$(nameof(Exp))" for (Exp, kw) in - ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) + ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) Random.seed!(4) ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) alg = TDVP(; alg_expand = Exp(; trscheme = truncrank(Dstart), kw...), trscheme = truncrank(Dcap)) From cfefd510b074858d9e2e16ea882bb66cf2618671 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 09:37:07 -0400 Subject: [PATCH 16/20] remove unused field --- src/algorithms/changebonds/sketchedexpand.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index a4b870ff0..310b273bd 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -20,8 +20,8 @@ The `warmstart` and state-preserving behaviour match [`OptimalExpand`](@ref). $(TYPEDFIELDS) """ @kwdef struct SketchedExpand{S} <: Algorithm - "algorithm used for the singular value decomposition" - alg_svd::S = Defaults.alg_svd() + "algorithm used to orthonormalize the sketched complement (passed as the `alg` of `left_orth!`/`right_orth!`); `nothing` selects QR without oversampling and an SVD-based decomposition otherwise" + alg_orth::S = nothing "algorithm used for truncating the expanded space" trscheme::TruncationStrategy @@ -91,9 +91,9 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk # orthonormalize the sketched right complement, truncating away the oversampling padding B = project_complement_right!(_transpose_tail(Y), ARtt) if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices - L, ar_re = right_orth!(B) + L, ar_re = right_orth!(B; alg = alg.alg_orth) else - L, ar_re = right_orth!(B; trunc = truncspace(Vk)) + L, ar_re = right_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) end # optimal vectors at site+1 @@ -135,9 +135,9 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske # orthonormalize the sketched left complement, truncating away the oversampling padding B = project_complement!(Y, left) if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices - Q, R = left_orth!(B) + Q, R = left_orth!(B; alg = alg.alg_orth) else - Q, R = left_orth!(B; trunc = truncspace(Vk)) + Q, R = left_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) end # optimal vectors at site-1 From a946d5629a2627e05f6e3fc9e6553886d766b914 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 10:50:02 -0400 Subject: [PATCH 17/20] remove warmstart and fix test errors --- src/algorithms/changebonds/optimalexpand.jl | 43 ++++++-------------- src/algorithms/changebonds/sketchedexpand.jl | 39 ++++++------------ src/algorithms/timestep/tdvp.jl | 5 +-- test/algorithms/changebonds.jl | 10 ----- test/algorithms/groundstate.jl | 21 ++-------- 5 files changed, 31 insertions(+), 87 deletions(-) diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 5a44d6326..c80a0c94c 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -5,14 +5,8 @@ An algorithm that expands the given mps as described in [Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018), by selecting the dominant contributions of a two-site updated MPS tensor, orthogonal to the original ψ. -By default the expansion does not alter the state: the added directions are connected through a -zero block, so that the expanded state is identical to the original one (as required for e.g. -TDVP). When `warmstart = true`, the new directions are instead seeded with the (physically -scaled) two-site gradient itself, so the expanded state already moves toward the optimal -two-site update. This changes the state and is therefore only useful for ground-state search -(e.g. as the `alg_expand` strategy of [`DMRG`](@ref)), where it warm-starts the subsequent -single-site optimization; the injected amplitude scales with the gradient norm and so vanishes -automatically at convergence. +The expansion does not alter the state: the added directions are connected through a zero block, +so that the expanded state is identical to the original one (as required for e.g. TDVP). !!! note [`changebonds!`](@ref) is only defined for `FiniteMPS`, and modifies both the state and its environment. @@ -28,9 +22,6 @@ $(TYPEDFIELDS) "algorithm used for truncating the expanded space" trscheme::TruncationStrategy - "whether to seed the new directions with the two-site gradient (warm start, alters the state) instead of a zero block" - warmstart::Bool = false - "setting for how much information is displayed" verbosity::Int = Defaults.verbosity end @@ -112,20 +103,15 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) gradnorm = norm(g2) - U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) @infov 4 "bond expansion" site dir = :right ϵ_select ϵ_2site = gradnorm / norm(AC2) # optimal vectors at site+1 ar_re = Vᴴ * NR - if alg.warmstart - # seed the new left directions with the (physically scaled) gradient instead of a zero - # block, warm-starting the subsequent optimization (alters the state) - nal, nc = left_gauge(catdomain(left, gradnorm * (NL * U * S))) - else - # embed `left` into the enlarged domain (zero weight in the new directions) - nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) - nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) - end + # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state + # unchanged + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1)) + nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) normalize && normalize!(nc) @@ -146,21 +132,16 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Opt # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) gradnorm = norm(g2) - U, S, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) + U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) @infov 4 "bond expansion" site dir = :left ϵ_select ϵ_2site = gradnorm / norm(AC2) # optimal vectors at site-1 Q = NL * U right_tail = _transpose_tail(right) - if alg.warmstart - # seed the new right directions with the (physically scaled) gradient instead of a zero - # block, warm-starting the subsequent optimization (alters the state) - nc, Qr = lq_compact!(catcodomain(right_tail, gradnorm * (S * Vᴴ * NR))) - else - # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) - nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) - nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) - end + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new + # directions), leaving the state unchanged + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) AL_exp = catdomain(left, Q) normalize && normalize!(nc) diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index 310b273bd..b044dc64e 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -8,7 +8,7 @@ single-site cost using the randomized "shrewd selection" of Controlled Bond Expa is folded into the effective environment, collapsing the large bond before the two-site update is ever formed, and the dominant directions are read off a small singular value decomposition. -The `warmstart` and state-preserving behaviour match [`OptimalExpand`](@ref). +The state-preserving behaviour matches [`OptimalExpand`](@ref). !!! note Only defined for `FiniteMPS` (through [`changebond!`](@ref)), so it can be used standalone or @@ -29,9 +29,6 @@ $(TYPEDFIELDS) "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)" oversampling::Int = 0 - "whether to seed the new directions with the sketched two-site gradient (warm start, alters the state) instead of a zero block" - warmstart::Bool = false - "setting for how much information is displayed" verbosity::Int = Defaults.verbosity end @@ -91,21 +88,16 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk # orthonormalize the sketched right complement, truncating away the oversampling padding B = project_complement_right!(_transpose_tail(Y), ARtt) if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices - L, ar_re = right_orth!(B; alg = alg.alg_orth) + _, ar_re = right_orth!(B; alg = alg.alg_orth) else - L, ar_re = right_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) + _, ar_re = right_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) end # optimal vectors at site+1 - if alg.warmstart - # seed the new left directions with the (physically scaled) gradient instead of a zero - # block, warm-starting the subsequent optimization (alters the state) - nal, nc = left_gauge(catdomain(left, Q * L)) - else - # embed `left` into the enlarged domain (zero weight in the new directions) - nal_space = codomain(left) ← (only(domain(left)) ⊕ space(ar_re, 1)) - nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) - end + # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state + # unchanged + nal_space = codomain(left) ← (only(domain(left)) ⊕ space(ar_re, 1)) + nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left)) nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re)) normalize && normalize!(nc) @@ -135,22 +127,17 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske # orthonormalize the sketched left complement, truncating away the oversampling padding B = project_complement!(Y, left) if dim(Vℓ) == dim(Vk) # no oversampling: a plain QR/LQ suffices - Q, R = left_orth!(B; alg = alg.alg_orth) + Q, _ = left_orth!(B; alg = alg.alg_orth) else - Q, R = left_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) + Q, _ = left_orth!(B; trunc = truncspace(Vk), alg = alg.alg_orth) end # optimal vectors at site-1 right_tail = _transpose_tail(right) - if alg.warmstart - # seed the new right directions with the (physically scaled) gradient instead of a zero - # block, warm-starting the subsequent optimization (alters the state) - nc, Qr2 = lq_compact!(catcodomain(right_tail, R * Qr_o)) - else - # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions) - nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) - nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) - end + # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new + # directions), leaving the state unchanged + nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail) + nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail)) AL_exp = catdomain(left, Q) normalize && normalize!(nc) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 396a68482..e99041a22 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -7,9 +7,8 @@ For finite MPS, setting `alg_expand` to a bond-expansion algorithm (e.g. [`Optim [`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to the current state ahead of each local integration, recovering Controlled Bond Expansion (CBE) TDVP and lifting the fixed-bond limitation of plain single-site TDVP. A truncating `trscheme` is then required to cut -the enlarged bond back down (selecting the truncated-SVD gauge). The expansion must be -state-preserving, so leave the expander's `warmstart = false` (the default); `warmstart = true` -injects a gradient and corrupts the evolution. +the enlarged bond back down (selecting the truncated-SVD gauge). The expansion is +state-preserving, as required for a consistent time evolution. !!! note Real-time evolution preserves the norm: neither the bond expansion nor the truncation diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl index a085f587e..a2aea1c08 100644 --- a/test/algorithms/changebonds.jl +++ b/test/algorithms/changebonds.jl @@ -91,16 +91,6 @@ end ) @test dot(state, state_se) ≈ 1 atol = 1.0e-8 - # `warmstart` seeds the new directions with the two-site gradient: the expansion is no longer - # norm-preserving, but the resulting state is still normalized - for (Exp, kw) in ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) - state_ws, _ = changebonds( - state, H, Exp(; trscheme = truncrank(dim(Dspace) * dim(Dspace)), warmstart = true, kw...) - ) - @test dot(state_ws, state_ws) ≈ 1 atol = 1.0e-8 - @test !isapprox(abs(dot(state, state_ws)), 1; atol = 1.0e-4) - end - state_tr = changebonds(state_oe, SvdCut(; trscheme = truncrank(dim(Dspace)))) @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5)) diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index 1338669ce..605d45814 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -88,11 +88,13 @@ verbosity_conv = 1 end @testset "CBEDMRG (SketchedExpand)" begin - # randomized bond expansion at single-site cost + # randomized bond expansion at single-site cost. The sketch is redrawn every sweep, so an + # aggressive expansion (a large fraction of the bond) keeps the single-site Galerkin error + # noisy; a gentle per-sweep increment lets it converge like the deterministic expanders. Random.seed!(1234) ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) v₀ = variance(ψ₀, H) - expand = SketchedExpand(; trscheme = truncrank(D ÷ 2), oversampling = 4) + expand = SketchedExpand(; trscheme = truncrank(2), oversampling = 4) trscheme = truncrank(D) # test logging @@ -113,21 +115,6 @@ verbosity_conv = 1 @test dim(left_virtualspace(ψ, L ÷ 2)) == D end - @testset "CBEDMRG warmstart $(nameof(Exp))" for (Exp, kw) in - ((OptimalExpand, (;)), (SketchedExpand, (; oversampling = 4))) - # warmstart seeds the expansion with the two-site gradient (alters the state); the - # ground-state search should still converge to the correct low-variance state - Random.seed!(2025) - ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2)) - expand = Exp(; trscheme = truncrank(D ÷ 2), warmstart = true, kw...) - ψ, envs, δ = find_groundstate( - ψ₀, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trscheme = truncrank(D)) - ) - @test sum(δ) ≈ 0 atol = 1.0e-3 - @test variance(ψ, H) < 1.0e-2 - @test dim(left_virtualspace(ψ, L ÷ 2)) == D - end - @testset "GradientGrassmann" begin ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D) v₀ = variance(ψ₀, H) From bb841e61abb9845c1d95c93dc2134f52f3906fa5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 11:28:26 -0400 Subject: [PATCH 18/20] more simplification --- src/algorithms/changebonds/changebonds.jl | 1 - src/algorithms/changebonds/optimalexpand.jl | 17 +++-------- src/algorithms/changebonds/randexpand.jl | 17 +++-------- src/algorithms/changebonds/sketchedexpand.jl | 30 ++++---------------- 4 files changed, 14 insertions(+), 51 deletions(-) diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl index c8b8bd17a..3a5dbeb5e 100644 --- a/src/algorithms/changebonds/changebonds.jl +++ b/src/algorithms/changebonds/changebonds.jl @@ -17,7 +17,6 @@ function changebonds! end changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ Expand a single bond of `ψ` in place by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the enriched bond. -Diagnostic error measures of the expansion (e.g. the discarded selection weight) are reported through the logging system at `@infov 4` rather than returned. The sweep direction `dir` is a `Val(:right)` or `Val(:left)` used for dispatch. For `Val(:right)` the bond `(site, site + 1)` is enriched on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions; for `Val(:left)` the mirror is applied to bond `(site - 1, site)`. diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index c80a0c94c..73da47899 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -21,9 +21,6 @@ $(TYPEDFIELDS) "algorithm used for truncating the expanded space" trscheme::TruncationStrategy - - "setting for how much information is displayed" - verbosity::Int = Defaults.verbosity end # Simple wrapper to convert between diffrent type of InifniteMPS. @@ -102,9 +99,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) - gradnorm = norm(g2) - _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) - @infov 4 "bond expansion" site dir = :right ϵ_select ϵ_2site = gradnorm / norm(AC2) + _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) # optimal vectors at site+1 ar_re = Vᴴ * NR @@ -131,9 +126,7 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Opt # select the dominant directions in the complement of the current state g2 = adjoint(NL) * AC2 * adjoint(NR) - gradnorm = norm(g2) - U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) - @infov 4 "bond expansion" site dir = :left ϵ_select ϵ_2site = gradnorm / norm(AC2) + U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) # optimal vectors at site-1 Q = NL * U @@ -154,10 +147,8 @@ changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ changebonds!(copy(ψ), H, alg, envs) function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ)) - LoggingExtras.withlevel(; alg.verbosity) do - for i in 1:(length(ψ) - 1) - changebond!(i, Val(:right), ψ, H, alg, envs) - end + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, H, alg, envs) end return ψ, envs end diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index f7cf68699..183830e36 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -23,9 +23,6 @@ $(TYPEDFIELDS) "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy] the expanded space" trscheme::TruncationStrategy - - "setting for how much information is displayed" - verbosity::Int = Defaults.verbosity end function changebonds!(ψ::InfiniteMPS, alg::RandExpand) @@ -73,10 +70,8 @@ end function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) # the expansion directions are sampled from a randomized two-site update, so no Hamiltonian # or environments are required - LoggingExtras.withlevel(; alg.verbosity) do - for i in 1:(length(ψ) - 1) - changebond!(i, Val(:right), ψ, nothing, alg, nothing) - end + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, nothing, alg, nothing) end return normalize!(ψ) @@ -94,9 +89,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra # select the dominant directions in the complement of the current state g2 = adjoint(NL) * ac2 * adjoint(NR) - gradnorm = norm(g2) - _, _, Vᴴ, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) - @infov 4 "bond expansion" site dir = :right ϵ_select ϵ_2site = gradnorm / norm(ac2) + _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) # optimal vectors at site+1, zero weight at site ar_re = Vᴴ * NR @@ -122,9 +115,7 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ran # select the dominant directions in the complement of the current state g2 = adjoint(NL) * ac2 * adjoint(NR) - gradnorm = norm(g2) - U, _, _, ϵ_select = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) - @infov 4 "bond expansion" site dir = :left ϵ_select ϵ_2site = gradnorm / norm(ac2) + U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd) # optimal vectors at site-1, zero weight at site Q = NL * U diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl index b044dc64e..9e89b3b81 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -28,9 +28,6 @@ $(TYPEDFIELDS) "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)" oversampling::Int = 0 - - "setting for how much information is displayed" - verbosity::Int = Defaults.verbosity end """ @@ -47,19 +44,6 @@ function sketch_space(V, alg::SketchedExpand) return infimum(V, Vk ⊕ Vp), Vk end -""" - complement_space(Vfull::ElementarySpace, image::ElementarySpace) - -The orthogonal complement `Vfull ⊖ image`, computed sector-wise from the spaces alone (no -factorization). Used to size and cap the random sketch. -""" -function complement_space(Vfull::ElementarySpace, image::ElementarySpace) - @assert !isdual(Vfull) - pairs = [c => (dim(Vfull, c) - dim(image, c)) for c in sectors(Vfull)] - filter!(p -> last(p) > 0, pairs) - return typeof(Vfull)(pairs) -end - # Finite system # ------------- # Rather than forming the two-site update, a random sketch of the orthogonal complement is folded @@ -73,8 +57,8 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk ARtt = _transpose_tail(right) # AR is already right-isometric # nothing to add when either complement is empty (e.g. edge bonds) - compL = complement_space(fuse(codomain(AL)), only(domain(AL))) - compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) + compL = fuse(codomain(AL)) ⊖ only(domain(AL)) + compR = fuse(domain(ARtt)) ⊖ only(codomain(ARtt)) (dim(compL) == 0 || dim(compR) == 0) && return ψ # random sketch of the left complement, folded into the left environment @@ -111,8 +95,8 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ske _, ARtt = right_orth!(_transpose_tail(right; copy = true); trunc = notrunc()) # local right-isometric form # nothing to add when either complement is empty (e.g. edge bonds) - compL = complement_space(fuse(codomain(left)), only(domain(left))) - compR = complement_space(fuse(domain(ARtt)), only(codomain(ARtt))) + compL = fuse(codomain(left)) ⊖ only(domain(left)) + compR = fuse(domain(ARtt)) ⊖ only(codomain(ARtt)) (dim(compL) == 0 || dim(compR) == 0) && return ψ # random sketch of the right complement, folded into the right environment @@ -150,10 +134,8 @@ changebonds(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments( changebonds!(copy(ψ), H, alg, envs) function changebonds!(ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs = environments(ψ, H, ψ)) - LoggingExtras.withlevel(; alg.verbosity) do - for i in 1:(length(ψ) - 1) - changebond!(i, Val(:right), ψ, H, alg, envs) - end + for i in 1:(length(ψ) - 1) + changebond!(i, Val(:right), ψ, H, alg, envs) end return ψ, envs end From 9e64a190ecb9192eb83fc7895cbe17b3187d5423 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 14:06:52 -0400 Subject: [PATCH 19/20] remove remnant code --- src/algorithms/timestep/tdvp.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index e99041a22..9ec1cb499 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -126,10 +126,7 @@ function timestep!( # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) - isnothing(alg.alg_expand) || - LoggingExtras.withlevel(; alg.alg_expand.verbosity) do - changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) - end + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) @@ -155,10 +152,7 @@ function timestep!( # sweep right to left for i in length(ψ):-1:2 # 1. optionally expand the bond ahead of the local update (CBE) - isnothing(alg.alg_expand) || - LoggingExtras.withlevel(; alg.alg_expand.verbosity) do - changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) - end + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) From 5395a896bb6ea40c2e2c6c1a316b4c769e58dcf1 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 30 Jun 2026 16:16:43 -0400 Subject: [PATCH 20/20] Guard CBE bond expansion in finite TDVP timestep! MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Restore the `isnothing(alg.alg_expand)` guard around `changebond!` in the finite TDVP sweep. It was dropped together with the verbosity `withlevel` wrapper during cleanup, so `TDVP()` (default `alg_expand = nothing`) called `changebond!(..., nothing, ...)` — which has no method — erroring for any `AbstractFiniteMPS` (e.g. WindowMPS evolution with an InfiniteMPOHamiltonian). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/timestep/tdvp.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 9ec1cb499..993782e54 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -126,7 +126,8 @@ function timestep!( # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) - changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + isnothing(alg.alg_expand) || + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) @@ -152,7 +153,8 @@ function timestep!( # sweep right to left for i in length(ψ):-1:2 # 1. optionally expand the bond ahead of the local update (CBE) - changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + isnothing(alg.alg_expand) || + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs)