From 924c5e0215ab013c4c9eda54901f47fb8b84d2ef Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 14:56:42 -0400 Subject: [PATCH 01/18] refactor fixedpoint to avoid warnings --- src/algorithms/fixedpoint.jl | 31 +++++++++++-------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/src/algorithms/fixedpoint.jl b/src/algorithms/fixedpoint.jl index 2248d416d..26b692e11 100644 --- a/src/algorithms/fixedpoint.jl +++ b/src/algorithms/fixedpoint.jl @@ -1,32 +1,23 @@ # wrapper around KrylovKit.jl's eigsolve function +function fixedpoint(A, x₀, which::Symbol; kwargs...) + alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...) + return fixedpoint(A, x₀, which, alg) +end + """ - fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec - fixedpoint(A, x₀, which::Symbol, alg) -> val, vec + fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec, info + fixedpoint(A, x₀, which::Symbol, alg) -> val, vec, info -Compute the fixed point of a linear operator `A` using the specified eigensolver `alg`. The -fixedpoint is assumed to be unique. +Compute the fixed point of a given linear operator `A` with initial guess `x₀`. +The dominant eigenvector is assumed to be unique. """ function fixedpoint(A, x₀, which::Symbol, alg::Lanczos) vals, vecs, info = eigsolve(A, x₀, 1, which, alg) - - info.converged == 0 && - @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])" - - return vals[1], vecs[1] + return vals[1], vecs[1], info end - function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi) TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg) - - info.converged == 0 && - @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])" size(TT, 2) > 1 && TT[2, 1] != 0 && @warnv 1 "non-unique fixed point detected" - - return vals[1], vecs[1] -end - -function fixedpoint(A, x₀, which::Symbol; kwargs...) - alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...) - return fixedpoint(A, x₀, which, alg) + return vals[1], vecs[1], info end From 174d86463e4633117c883761919a8688766eae45 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 9 Jul 2026 12:42:13 -0400 Subject: [PATCH 02/18] refactor DMRG implementation --- src/algorithms/groundstate/dmrg.jl | 127 +++++++++++++++++------------ 1 file changed, 75 insertions(+), 52 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 61f3097a0..fb5cb99cc 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -53,8 +53,10 @@ function DMRG(; 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 + # single-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass + # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default). + alg_eigsolve′ = alg_eigsolve isa NamedTuple ? + Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve # 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) @@ -64,83 +66,100 @@ function DMRG(; return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge) end + +function dmrg_update!( + site, direction, + ψ, O, alg, envs, + ϵ_global, ϵ_trunc, decay_rate, + timeroutput + ) + ϵ_local = calc_galerkin(site, ψ, O, ψ, envs) + + # 1. expand + isnothing(alg.alg_expand) || + @timeit timeroutput "expand" changebond!(site, direction, ψ, O, alg.alg_expand, envs) + + # 2. local update + alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) + λ, AC′, info = @timeit timeroutput "AC_eigsolve" begin + H_effective = AC_hamiltonian(site, ψ, O, ψ, envs) + fixedpoint(H_effective, ψ.AC[site], :SR, alg_eigsolve) + end + + # 3. gauge + ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, AC′, alg.alg_gauge; normalize = true) + + # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) + decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) + + @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local + + return ψ, ϵ_local, ϵ_trunc, decay_rate +end + + function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ)) - ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) - ϵ = maximum(ϵs) log = IterLog("DMRG") timeroutput = TimerOutput("DMRG") alg.verbosity > 3 || disable_timer!(timeroutput) + Tr = real(scalartype(ψ)) + ϵ_locals = ones(Tr, length(ψ)) # local gradient norms + ϵ_global = maximum(ϵ_locals) + ϵ_truncs = zeros(Tr, length(ψ)) # local truncation error + decay_rates = zeros(length(ψ)) # local observed decay rate of eigensolver + LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs)) + @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) for iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) - - zerovector!(ϵs) @timeit timeroutput "sweep" begin # 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)) - - # 1. expand - isnothing(alg.alg_expand) || - @timeit timeroutput "expand" changebond!(pos, Val(:right), ψ, H, alg.alg_expand, envs) - - # 2. local update - @timeit timeroutput "AC_eigsolve" begin - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) - 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) + direction = Val(:right) + for site in 1:(length(ψ) - 1) + ψ, ϵ_locals[site], ϵ_truncs[site], decay_rates[site] = + dmrg_update!( + site, direction, + ψ, H, alg, envs, + ϵ_global, ϵ_truncs[site], decay_rates[site], + timeroutput + ) + ϵ_global = maximum(ϵ_locals) end # 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)) - - # 1. expand - isnothing(alg.alg_expand) || - @timeit timeroutput "expand" changebond!(pos, Val(:left), ψ, H, alg.alg_expand, envs) - - # 2. local update - @timeit timeroutput "AC_eigsolve" begin - Hac = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve) - end - - # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`) - @timeit timeroutput "gauge" right_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true) + direction = Val(:left) + for site in length(ψ):-1:2 + ψ, ϵ_locals[site], ϵ_truncs[site], decay_rates[site] = + dmrg_update!( + site, direction, + ψ, H, alg, envs, + ϵ_global, ϵ_truncs[site], decay_rates[site], + timeroutput + ) + ϵ_global = maximum(ϵ_locals) end end - ϵ = maximum(ϵs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( iter, ψ, H, envs )::Tuple{typeof(ψ), typeof(envs)} - if ϵ <= alg.tol + if ϵ_global <= alg.tol @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) break end if iter == alg.maxiter @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) else - @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) end end end - return ψ, envs, ϵ + return ψ, envs, ϵ_global end - """ $(TYPEDEF) @@ -192,8 +211,10 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm LoggingExtras.withlevel(; alg.verbosity) do for iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) zerovector!(ϵs) + # two-site DMRG uses a plain (non-adaptive) local eigensolver whose tolerance is + # tightened per sweep with the current global error via `updatetol`. + eig_alg = updatetol(alg.alg_eigsolve, iter, ϵ) @timeit timeroutput "sweep" begin # left to right sweep @@ -202,7 +223,8 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm @timeit timeroutput "AC2_eigsolve" begin @plansor ac2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3] Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve) + _, newA2center, info = fixedpoint(Hac2, ac2, :SR, eig_alg) + @debug "DMRG2 local update" pos numops = info.numops normres = first(info.normres) end @timeit timeroutput "svd_trunc" begin al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd) @@ -222,7 +244,8 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environm @timeit timeroutput "AC2_eigsolve" begin @plansor ac2[-1 -2; -3 -4] := ψ.AL[pos][-1 -2; 1] * ψ.AC[pos + 1][1 -4; -3] Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve) + _, newA2center, info = fixedpoint(Hac2, ac2, :SR, eig_alg) + @debug "DMRG2 local update" pos numops = info.numops normres = first(info.normres) end @timeit timeroutput "svd_trunc" begin al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd) From abe50dd75dadb2ae026346a353ce78d75a65e1af Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 9 Jul 2026 12:44:02 -0400 Subject: [PATCH 03/18] refactor gauging --- src/algorithms/changebonds/optimalexpand.jl | 2 +- src/algorithms/changebonds/randexpand.jl | 2 +- src/algorithms/changebonds/sketchedexpand.jl | 4 +- src/states/orthoview.jl | 55 +++++++++++++------- 4 files changed, 41 insertions(+), 22 deletions(-) diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index d88b3aa3b..87d54ee5b 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -114,7 +114,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op # 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)) + 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/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 2ab840ffe..85bbfa81c 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -98,7 +98,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 = left_gauge(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 9e89b3b81..deade493f 100644 --- a/src/algorithms/changebonds/sketchedexpand.jl +++ b/src/algorithms/changebonds/sketchedexpand.jl @@ -53,7 +53,7 @@ end function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true) left = ψ.AC[site] right = ψ.AR[site + 1] - AL, _ = left_gauge(left) # local left-isometric form + 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) @@ -81,7 +81,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk # 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)) + 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/states/orthoview.jl b/src/states/orthoview.jl index bb2338340..c713247b7 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -60,7 +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], v.parent.ARs[j] = right_gauge(v.parent.ACs[j]) + 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 @@ -75,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] = left_gauge(v.parent.ACs[j]) + 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 @@ -88,9 +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], v.parent.ARs[i + 1] = right_gauge(v.parent.AC[i + 1]) + 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] = left_gauge(v.parent.AC[i]) + v.parent.ALs[i], v.parent.Cs[i + 1], _ = left_gauge(v.parent.AC[i]) end end @@ -223,8 +223,8 @@ end # Gauging routines # ---------------- @doc """ - left_gauge(AC, [alg]) -> AL, C - right_gauge(AC, [alg]) -> C, AR + 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 @@ -233,28 +233,38 @@ standard MPS-tensor form. `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. + +Also returns the truncation error `ϵ`: the 2-norm of the discarded singular values from the +truncated SVD, or `0` for a norm-preserving QR/LQ gauge. """ 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)) +function left_gauge(AC, alg) + AL, C = left_orth(AC; alg) + return AL, C, zero(real(scalartype(AC))) +end +function left_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) + U, S, Vᴴ, ϵ = svd_trunc(AC, alg) + C = LinearAlgebra.lmul!(S, Vᴴ) # C = S * Vᴴ, matching `LeftOrthViaSVD` + return U, C, ϵ +end 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) + return C, _transpose_front(AR), zero(real(scalartype(AC))) end function right_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) - C, AR = right_orth(_transpose_tail(AC); alg = MatrixAlgebraKit.RightOrthViaSVD(alg)) - return C, _transpose_front(AR) + U, S, Vᴴ, ϵ = svd_trunc(_transpose_tail(AC), alg) + C = LinearAlgebra.rmul!(U, S) # C = U * S, matching `RightOrthViaSVD` + return C, _transpose_front(Vᴴ), ϵ end @doc """ - left_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ - right_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ + 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, @@ -264,19 +274,28 @@ be a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) to trun 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. + +Also returns the truncation error `ϵ` of the factorization: the 2-norm of the discarded singular +values from a truncated SVD gauge, or `0` for a norm-preserving QR gauge. """ 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) + AL, C, ϵ = left_gauge(AC, alg) normalize && normalize!(C) ψ.AC[pos] = (AL, C) - return ψ + return ψ, ϵ end function right_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false) - C, AR = right_gauge(AC, alg) + C, AR, ϵ = right_gauge(AC, alg) normalize && normalize!(C) ψ.AC[pos] = (C, AR) - return ψ + return ψ, ϵ +end + +function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC, alg = Defaults.alg_orth(); kwargs...) where {Dir} + Dir === :right && return left_gauge!(ψ, pos, AC, alg; kwargs...) + Dir === :left && return right_gauge!(ψ, pos, AC, alg; kwargs...) + return throw(ArgumentError(lazy"invalid direction `$Dir`")) end From 1caab6f8d91f4a7fd12d2f0fee8d9529ff6e08fe Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 9 Jul 2026 12:44:33 -0400 Subject: [PATCH 04/18] add AdaptiveKrylov solver --- src/MPSKit.jl | 2 +- src/utility/defaults.jl | 7 +- src/utility/dynamictols.jl | 147 +++++++++++++++++++++++++++++++++++++ 3 files changed, 153 insertions(+), 3 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index b7387f3da..2ef8cf900 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -52,7 +52,7 @@ export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename # unexported using Compat: @compat -@compat public DynamicTols +@compat public DynamicTols, DynamicTol, AdaptiveKrylov @compat public VERBOSE_NONE, VERBOSE_WARN, VERBOSE_CONV, VERBOSE_ITER, VERBOSE_ALL @compat public IterLog, loginit!, logiter!, logfinish!, logcancel! diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index ddfd02afe..11cf238eb 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -7,7 +7,7 @@ module Defaults import KrylovKit: GMRES, Arnoldi, Lanczos using OhMyThreads -using ..MPSKit: DynamicTol +using ..MPSKit: DynamicTol, AdaptiveKrylov using TensorKit: TensorKit using MatrixAlgebraKit: DefaultAlgorithm, Householder @@ -54,8 +54,11 @@ end function alg_eigsolve(; ishermitian = true, tol = tol, maxiter = maxiter, verbosity = 0, eager = true, krylovdim = krylovdim, dynamic_tols = dynamic_tols, tol_min = tol_min, - tol_max = tol_max, tol_factor = eigs_tolfactor + tol_max = tol_max, tol_factor = eigs_tolfactor, adaptive = false, adaptive_kwargs... ) + # the per-bond adaptive controller (`AdaptiveKrylov`) retunes krylovdim/maxiter/tol per site + # and supersedes the (global, tol-only) `DynamicTol` wrapper; it is the default for `DMRG`. + adaptive && return AdaptiveKrylov(; ishermitian, adaptive_kwargs...) alg = ishermitian ? Lanczos(; tol, maxiter, eager, krylovdim, verbosity) : Arnoldi(; tol, maxiter, eager, krylovdim, verbosity) return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index 7a67a9a01..ade854164 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -4,8 +4,10 @@ import ..MPSKit: Algorithm using Accessors using DocStringExtensions using MatrixAlgebraKit: DefaultAlgorithm +using KrylovKit: KrylovKit, Lanczos, Arnoldi export updatetol, DynamicTol +export AdaptiveKrylov, instantiate_algorithm @doc """ updatetol(alg, iter, ϵ) @@ -77,4 +79,149 @@ function _updatetol(alg::DefaultAlgorithm, tol::Real) return DefaultAlgorithm(; kwargs...) end +# Set several solver parameters at once (tol, krylovdim, maxiter). Used by the adaptive +# controller, which retunes more than just the tolerance. `Lanczos`/`Arnoldi` expose these +# as plain fields, so `Accessors.setproperties` rebuilds the immutable struct in one shot. +function _set_eigsolve_params(alg::Union{Lanczos, Arnoldi}; kwargs...) + return Accessors.setproperties(alg, NamedTuple(kwargs)) +end + +# ============================================================ +# Adaptive Krylov controller (per-bond, stateful) +# ============================================================ + +""" +$(TYPEDEF) + +Adaptive controller for local eigensolvers in DMRG-like algorithms. + +This controller attempts to balance minimizing the number of local function applications with +minimizing the total number of iterations in order to obtain the globally fastest convergence. +This is driven by the local and global gradient norm, the truncation error and the measured +decay rate of previous iterations in an attempt to obtain fast convergence for gapped systems +while avoiding stagnation for gapless ones. + +## Fields + +$(TYPEDFIELDS) + +See also [`instantiate_algorithm`](@ref). +""" +struct AdaptiveKrylov{T, O <: KrylovKit.Orthogonalizer} <: Algorithm + "orthogonalizer passed to the instantiated `Lanczos`/`Arnoldi`" + orth::O + + "minimal Krylov subspace dimension (conservative/cold-start default)" + krylovdim_min::Int + "maximal Krylov subspace dimension" + krylovdim_max::Int + + "minimal number of restart iterations" + iter_min::Int + "maximal number of restart iterations" + iter_max::Int + + "lower bound on the dynamically chosen tolerance" + tol_min::Float64 + "upper bound on the dynamically chosen tolerance" + tol_max::Float64 + "factor on the local truncation error, sets the inner-solve tolerance floor" + truncation_factor::Float64 + "factor on the global gradient norm, sets the convergence-driven tolerance" + tol_factor::Float64 + + "whether to use the eager (early-stopping) solver mode" + eager::Bool + "verbosity of the instantiated eigensolver" + verbosity::Int + + "hermitian flag (`Val(true)` → `Lanczos`, `Val(false)` → `Arnoldi`)" + ishermitian::Val{T} + + function AdaptiveKrylov{T, O}( + orth::O, krylovdim_min, krylovdim_max, iter_min, iter_max, + tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, ishermitian::Val{T} + ) where {T, O <: KrylovKit.Orthogonalizer} + 0 < krylovdim_min <= krylovdim_max || + throw(ArgumentError("need 0 < krylovdim_min ≤ krylovdim_max, got ($krylovdim_min, $krylovdim_max)")) + 1 <= iter_min <= iter_max || + throw(ArgumentError("need 1 ≤ iter_min ≤ iter_max, got ($iter_min, $iter_max)")) + 0 <= tol_min <= tol_max || + throw(ArgumentError("need 0 ≤ tol_min ≤ tol_max, got ($tol_min, $tol_max)")) + truncation_factor >= 0 || + throw(ArgumentError("truncation_factor must be non-negative, got $truncation_factor")) + tol_factor >= 0 || + throw(ArgumentError("tol_factor must be non-negative, got $tol_factor")) + return new{T, O}( + orth, krylovdim_min, krylovdim_max, iter_min, iter_max, + tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, ishermitian + ) + end +end + +function AdaptiveKrylov(; + ishermitian::Bool = true, orth::KrylovKit.Orthogonalizer = KrylovKit.KrylovDefaults.orth, + krylovdim_min::Int = 3, krylovdim_max::Int = 8, + iter_min::Int = 1, iter_max::Int = 2, + tol_min::Real = 1.0e-12, tol_max::Real = 1.0e-2, + truncation_factor::Real = 1.0e-1, tol_factor::Real = 1.0e-1, + eager::Bool = true, verbosity::Int = 0 + ) + return AdaptiveKrylov{ishermitian, typeof(orth)}( + orth, krylovdim_min, krylovdim_max, iter_min, iter_max, + tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, Val(ishermitian) + ) +end + +""" + instantiate_algorithm(alg, decay_rate, g_local, g_global, eps_trunc) + +Turn a (possibly adaptive) local-eigensolver specification `alg` into a concrete KrylovKit +algorithm for the current site, using the measured `decay_rate` of the previous solve, the +local gradient norm `g_local`, the global gradient norm `g_global` and the local truncation +error `eps_trunc`. The generic fallback returns `alg` unchanged (for plain `Lanczos`/`Arnoldi`). +""" +instantiate_algorithm(alg, args...) = alg + +# a `DynamicTol` wrapper is resolved to a concrete solver by tightening its tolerance with the +# global gradient norm (the per-sweep, tol-only "legacy" behaviour, now driven per site). +function instantiate_algorithm( + alg::DynamicTol, decay_rate::Real, g_local::Real, g_global::Real, eps_trunc::Real + ) + tol = clamp(g_global * alg.tol_factor, alg.tol_min, alg.tol_max) + return _updatetol(alg.alg, tol) +end + +function instantiate_algorithm( + alg::AdaptiveKrylov{T}, decay_rate::Real, g_local::Real, g_global::Real, eps_trunc::Real + ) where {T} + # 1. target tolerance: inexact inner solve depending on outer convergence, + # never below the truncation error we already incur. + trunc_tol = alg.truncation_factor * eps_trunc + conv_tol = alg.tol_factor * g_global + tol = clamp(max(trunc_tol, conv_tol), alg.tol_min, alg.tol_max) + + # the measured decay rate must be a genuine contraction in (0, 1); a non-positive or ≥ 1 + # rate signals an uninitialized or hard/stalling bond → fall back to the largest budget. + ρ = clamp(decay_rate, 0.0, 1.0) + if !(0 < ρ < 1) + krylovdim = iszero(ρ) ? alg.krylovdim_min : alg.krylovdim_max + maxiter = iszero(ρ) ? alg.iter_min : alg.iter_max + return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity) + end + + # 2. predicted matvecs from the *measured* decay and local convergence + # estimate through Lanczos bounds: going from g_local to tol in factors of ρ + R = max(g_local, alg.tol_min) / tol + nmatvecs = round(Int, log(R) / log(inv(ρ))) + + # 3. krylovdim and maxiter from matvec budget: + # prioritize krylovdim and compensate with maxiter + # thick restarts keep 3/5 of krylovdim, i.e. ~2/5 new vectors per cycle + krylovdim = clamp(nmatvecs, alg.krylovdim_min, alg.krylovdim_max) + maxiter = clamp(cld(5nmatvecs, 2krylovdim) - 3, alg.iter_min, alg.iter_max) + + return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity) +end + end From 89cf07a5c19bfc03a910f6de13b2bd1cf0363cd8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 9 Jul 2026 20:37:28 -0400 Subject: [PATCH 05/18] refactor DMRG2 to be in line with DMRG --- src/algorithms/groundstate/dmrg.jl | 224 +++++++++++++---------------- src/states/finitemps.jl | 11 +- src/states/orthoview.jl | 31 ++++ 3 files changed, 136 insertions(+), 130 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index fb5cb99cc..40ab74d3d 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -67,9 +67,9 @@ function DMRG(; end -function dmrg_update!( +function local_update!( site, direction, - ψ, O, alg, envs, + ψ, O, alg::DMRG, envs, ϵ_global, ϵ_trunc, decay_rate, timeroutput ) @@ -97,69 +97,6 @@ function dmrg_update!( return ψ, ϵ_local, ϵ_trunc, decay_rate end - -function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ)) - log = IterLog("DMRG") - timeroutput = TimerOutput("DMRG") - alg.verbosity > 3 || disable_timer!(timeroutput) - - Tr = real(scalartype(ψ)) - ϵ_locals = ones(Tr, length(ψ)) # local gradient norms - ϵ_global = maximum(ϵ_locals) - ϵ_truncs = zeros(Tr, length(ψ)) # local truncation error - decay_rates = zeros(length(ψ)) # local observed decay rate of eigensolver - - LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) - for iter in 1:(alg.maxiter) - @timeit timeroutput "sweep" begin - # left-to-right - direction = Val(:right) - for site in 1:(length(ψ) - 1) - ψ, ϵ_locals[site], ϵ_truncs[site], decay_rates[site] = - dmrg_update!( - site, direction, - ψ, H, alg, envs, - ϵ_global, ϵ_truncs[site], decay_rates[site], - timeroutput - ) - ϵ_global = maximum(ϵ_locals) - end - - # right-to-left - direction = Val(:left) - for site in length(ψ):-1:2 - ψ, ϵ_locals[site], ϵ_truncs[site], decay_rates[site] = - dmrg_update!( - site, direction, - ψ, H, alg, envs, - ϵ_global, ϵ_truncs[site], decay_rates[site], - timeroutput - ) - ϵ_global = maximum(ϵ_locals) - end - end - - ψ, envs = @timeit timeroutput "finalize" alg.finalize( - iter, ψ, H, envs - )::Tuple{typeof(ψ), typeof(envs)} - - if ϵ_global <= alg.tol - @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) - break - end - if iter == alg.maxiter - @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) - else - @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) - end - end - end - return ψ, envs, ϵ_global -end - """ $(TYPEDEF) @@ -169,7 +106,7 @@ Two-site DMRG algorithm for finding the dominant eigenvector. $(TYPEDFIELDS) """ -struct DMRG2{A, S, F} <: Algorithm +struct DMRG2{A, G, F} <: Algorithm "tolerance for convergence criterium" tol::Float64 @@ -182,11 +119,8 @@ struct DMRG2{A, S, F} <: Algorithm "algorithm used for the eigenvalue solvers" alg_eigsolve::A - "algorithm used for the singular value decomposition" - alg_svd::S - - "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update" - trscheme::TruncationStrategy + "factorization used for the post-update gauge: a truncated SVD (`alg_svd` with `trscheme`)" + alg_gauge::G "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" finalize::F @@ -197,88 +131,122 @@ function DMRG2(; alg_eigsolve = (;), alg_svd = Defaults.alg_svd(), trscheme, finalize = Defaults._finalize ) - alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : - alg_eigsolve - return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_svd, trscheme, finalize) + # two-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass + # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default). + alg_eigsolve′ = alg_eigsolve isa NamedTuple ? + Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve + # two-site DMRG always truncates the enlarged bond back down, so the gauge is a truncated SVD + alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme) + return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize) +end + +function local_update!( + pos, direction, + ψ, O, alg::DMRG2, envs, + ϵ_global, ϵ_trunc, decay_rate, + timeroutput + ) + Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs) + + kind = direction === Val(:right) ? :ACAR : :ALAC + HAC2 = normalize!(Heff * AC2(ψ, pos; kind)) + AC2′ = copy(HAC2) + project_complement!(AC2′, ψ.AL[pos]) + project_complement_right!(AC2′, _transpose_tail(ψ.AR[pos + 1])) + ϵ_local = norm(AC2′) + + # 1. local two-site update + alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) + newA2center, info = @timeit timeroutput "AC2_eigsolve" begin + _, newA2center, info = fixedpoint(Heff, HAC2, :SR, alg_eigsolve) + (newA2center, info) + end + + # 2. gauge: truncated SVD split back into single-site tensors and install; + # the discarded weight is the truncation error + ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, newA2center, alg.alg_gauge; normalize = true) + + # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) + decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) + + @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc + + return ψ, ϵ_local, ϵ_trunc, decay_rate end -function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environments(ψ, H, ψ)) - ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) - ϵ = maximum(ϵs) - log = IterLog("DMRG2") - timeroutput = TimerOutput("DMRG2") +# Per-algorithm sweep geometry: single-site DMRG updates all `length(ψ)` sites (endpoints once, +# interior twice), whereas two-site DMRG2 updates the `length(ψ) - 1` bonds. `_num_updates` +# gives the number of per-update bookkeeping slots and `_sweep_ranges` the forward/backward +# index ranges; everything else in the sweep is shared. +_num_updates(::DMRG, ψ) = length(ψ) +_num_updates(::DMRG2, ψ) = length(ψ) - 1 + +_sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2) +_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1) + +function find_groundstate!( + ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ) + ) + name = string(nameof(typeof(alg))) + log = IterLog(name) + timeroutput = TimerOutput(name) alg.verbosity > 3 || disable_timer!(timeroutput) + Tr = real(scalartype(ψ)) + n = _num_updates(alg, ψ) + ϵ_locals = ones(Tr, n) # local gradient norms + ϵ_global = maximum(ϵ_locals) + ϵ_truncs = zeros(Tr, n) # local truncation error + decay_rates = zeros(n) # local observed decay rate of eigensolver + fwd, bwd = _sweep_ranges(alg, ψ) + LoggingExtras.withlevel(; alg.verbosity) do + @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) for iter in 1:(alg.maxiter) - zerovector!(ϵs) - # two-site DMRG uses a plain (non-adaptive) local eigensolver whose tolerance is - # tightened per sweep with the current global error via `updatetol`. - eig_alg = updatetol(alg.alg_eigsolve, iter, ϵ) - @timeit timeroutput "sweep" begin - # left to right sweep - for pos in 1:(length(ψ) - 1) - local ac2, newA2center, al, c, ar - @timeit timeroutput "AC2_eigsolve" begin - @plansor ac2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3] - Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, newA2center, info = fixedpoint(Hac2, ac2, :SR, eig_alg) - @debug "DMRG2 local update" pos numops = info.numops normres = first(info.normres) - end - @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) - v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) - ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) - end - @timeit timeroutput "update_AC" begin - ψ.AC[pos] = (al, complex(c)) - ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) - end + # left-to-right + for pos in fwd + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = + local_update!( + pos, Val(:right), + ψ, H, alg, envs, + ϵ_global, ϵ_truncs[pos], decay_rates[pos], + timeroutput + ) + ϵ_global = maximum(ϵ_locals) end - # right to left sweep - for pos in (length(ψ) - 2):-1:1 - local ac2, newA2center, al, c, ar - @timeit timeroutput "AC2_eigsolve" begin - @plansor ac2[-1 -2; -3 -4] := ψ.AL[pos][-1 -2; 1] * ψ.AC[pos + 1][1 -4; -3] - Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, newA2center, info = fixedpoint(Hac2, ac2, :SR, eig_alg) - @debug "DMRG2 local update" pos numops = info.numops normres = first(info.normres) - end - @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) - v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) - ϵs[pos] = max(ϵs[pos], abs(1 - abs(v))) - end - @timeit timeroutput "update_AC" begin - ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) - ψ.AC[pos] = (al, complex(c)) - end + # right-to-left + for pos in bwd + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = + local_update!( + pos, Val(:left), + ψ, H, alg, envs, + ϵ_global, ϵ_truncs[pos], decay_rates[pos], + timeroutput + ) + ϵ_global = maximum(ϵ_locals) end end - ϵ = maximum(ϵs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( iter, ψ, H, envs )::Tuple{typeof(ψ), typeof(envs)} - if ϵ <= alg.tol + if ϵ_global <= alg.tol @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) break end if iter == alg.maxiter @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) else - @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs)) + @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) end end end - return ψ, envs, ϵ + return ψ, envs, ϵ_global end function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index 6f65dbbe5..f38989b9c 100644 --- a/src/states/finitemps.jl +++ b/src/states/finitemps.jl @@ -343,8 +343,15 @@ Base.@propagate_inbounds function Base.getindex(ψ::FiniteMPS, i::Int) end end -# TODO: check where gauge center is to determine efficient kind -AC2(psi::FiniteMPS, site::Int) = psi.AC[site] * _transpose_tail(psi.AR[site + 1]) +function AC2(psi::FiniteMPS, site::Int; kind = :ACAR) + if kind == :ACAR + return psi.AC[site] * _transpose_tail(psi.AR[site + 1]) + elseif kind == :ALAC + return psi.AL[site] * _transpose_tail(psi.AC[site + 1]) + else + throw(ArgumentError("Invalid kind: $kind")) + end +end f_if_not_missing(f, x) = ismissing(x) ? x : f(x) _copy_if_not_missing(x) = f_if_not_missing(copy, x) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index c713247b7..d75fde086 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -299,3 +299,34 @@ function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC, alg = Defaults. Dir === :left && return right_gauge!(ψ, pos, AC, alg; kwargs...) return throw(ArgumentError(lazy"invalid direction `$Dir`")) end + +@doc """ + gauge2!(ψ, pos, direction, AC2, alg; normalize = false) -> ψ, ϵ + +Two-site analogue of [`gauge!`](@ref): factor an updated two-site center tensor `AC2` spanning +sites `pos` and `pos+1` with the truncated SVD `alg` and install the resulting canonical tensors +into `ψ` in one step, shifting the gauge center past the bond. +(To the right for `direction = Val(:right)`, to the left for `Val(:left)`). + +Returns the truncation error `ϵ`, the 2-norm of the discarded singular values. Pass +`normalize = true` to renormalize the bond tensor, so `ψ` stays normalized after a local update +that changed its norm. +""" +function gauge2!( + ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC2, + alg::MatrixAlgebraKit.TruncatedAlgorithm; normalize::Bool = false + ) where {Dir} + al, c, ar, ϵ = svd_trunc!(AC2, alg) + normalize && normalize!(c) + C = complex(c) + if Dir === :right + ψ.AC[pos] = (al, C) + ψ.AC[pos + 1] = (C, _transpose_front(ar)) + elseif Dir === :left + ψ.AC[pos + 1] = (C, _transpose_front(ar)) + ψ.AC[pos] = (al, C) + else + throw(ArgumentError(lazy"invalid direction `$Dir`")) + end + return ψ, ϵ +end From 16797ea83d5be2c7a4116149e3f18278b79d9d09 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 11 Jul 2026 23:21:04 -0400 Subject: [PATCH 06/18] small type stability improvements --- src/operators/jordanmpotensor.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index acea13fe3..d9a763402 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -170,7 +170,11 @@ end TensorKit.space(W::JordanMPOTensor) = space(W.tensors) Base.eltype(::Type{T}) where {T <: JordanMPOTensor} = eltype(fieldtype(T, :tensors)) -function Base.getproperty(W::JordanMPOTensor, sym::Symbol) +Base.size(W::JordanMPOTensor) = size(W.tensors) +Base.size(W::JordanMPOTensor, i::Int) = size(W.tensors, i) +BlockTensorKit.eachspace(W::JordanMPOTensor) = eachspace(W.tensors) + +Base.@constprop :aggressive function Base.getproperty(W::JordanMPOTensor, sym::Symbol) sym === :A && return _jordan_A(W) sym === :B && return _jordan_B(W) sym === :C && return _jordan_C(W) From 69c0d16edf8e7b9f7c6534e64753edf8b9136b5f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 12 Jul 2026 08:18:10 -0400 Subject: [PATCH 07/18] small fixes --- src/algorithms/approximate/fvomps.jl | 2 +- src/states/orthoview.jl | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 5a9ec4d29..d5ffc6cc4 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/src/algorithms/approximate/fvomps.jl @@ -8,7 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment ϵ = 0.0 for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1] AC2′ = AC2_projection(pos, ψ, Oϕ, envs) - al, c, ar = svd_trunc!(AC2′; trunc = alg.trscheme, alg = alg.alg_svd) + al, c, ar, = svd_trunc!(AC2′, alg.alg_gauge) AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1]) ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2)) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index d75fde086..910507b7c 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -294,6 +294,14 @@ function right_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_or return ψ, ϵ end +@doc """ + gauge!(ψ, pos, direction, AC, [alg]; normalize = false) -> ψ, ϵ + +Direction-dispatching wrapper around [`left_gauge!`](@ref) / [`right_gauge!`](@ref): gauge an +updated center tensor `AC` at site `pos` and install it into `ψ`, shifting the gauge center past +`pos` to the right for `direction = Val(:right)` and to the left for `Val(:left)`. `alg` and +`normalize` are forwarded unchanged, and the truncation error `ϵ` is returned as-is. +""" function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC, alg = Defaults.alg_orth(); kwargs...) where {Dir} Dir === :right && return left_gauge!(ψ, pos, AC, alg; kwargs...) Dir === :left && return right_gauge!(ψ, pos, AC, alg; kwargs...) @@ -313,8 +321,8 @@ Returns the truncation error `ϵ`, the 2-norm of the discarded singular values. that changed its norm. """ function gauge2!( - ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC2, - alg::MatrixAlgebraKit.TruncatedAlgorithm; normalize::Bool = false + ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC2, alg; + normalize::Bool = false ) where {Dir} al, c, ar, ϵ = svd_trunc!(AC2, alg) normalize && normalize!(c) From 8d337929f10003ff17c320d158f9ba94c783015c Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 12 Jul 2026 13:28:20 -0400 Subject: [PATCH 08/18] restore old DMRG convergence behavior --- src/algorithms/groundstate/dmrg.jl | 41 +++++++++++++++++++----------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 40ab74d3d..5087ca728 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -81,20 +81,24 @@ function local_update!( # 2. local update alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) + ac_old = ψ.AC[site] λ, AC′, info = @timeit timeroutput "AC_eigsolve" begin H_effective = AC_hamiltonian(site, ψ, O, ψ, envs) - fixedpoint(H_effective, ψ.AC[site], :SR, alg_eigsolve) + fixedpoint(H_effective, ac_old, :SR, alg_eigsolve) end + # convergence error: how much the eigensolve moved the center, `1 - |⟨old|new⟩|`. + ϵ_conv = abs(1 - abs(dot(ac_old, AC′)) / (norm(ac_old) * norm(AC′))) + # 3. gauge ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, AC′, alg.alg_gauge; normalize = true) # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) - @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local + @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_conv - return ψ, ϵ_local, ϵ_trunc, decay_rate + return ψ, ϵ_local, ϵ_trunc, decay_rate, ϵ_conv end """ @@ -149,7 +153,8 @@ function local_update!( Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs) kind = direction === Val(:right) ? :ACAR : :ALAC - HAC2 = normalize!(Heff * AC2(ψ, pos; kind)) + ac2 = normalize!(AC2(ψ, pos; kind)) + HAC2 = normalize!(Heff * ac2) AC2′ = copy(HAC2) project_complement!(AC2′, ψ.AL[pos]) project_complement_right!(AC2′, _transpose_tail(ψ.AR[pos + 1])) @@ -162,6 +167,9 @@ function local_update!( (newA2center, info) end + # convergence error: how much the eigensolve moved the two-site center, `1 - |⟨old|new⟩|`. + ϵ_conv = abs(1 - abs(dot(ac2, newA2center)) / norm(newA2center)) + # 2. gauge: truncated SVD split back into single-site tensors and install; # the discarded weight is the truncation error ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, newA2center, alg.alg_gauge; normalize = true) @@ -169,9 +177,9 @@ function local_update!( # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) - @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc + @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc ϵ_conv - return ψ, ϵ_local, ϵ_trunc, decay_rate + return ψ, ϵ_local, ϵ_trunc, decay_rate, ϵ_conv end # Per-algorithm sweep geometry: single-site DMRG updates all `length(ψ)` sites (endpoints once, @@ -194,19 +202,21 @@ function find_groundstate!( Tr = real(scalartype(ψ)) n = _num_updates(alg, ψ) - ϵ_locals = ones(Tr, n) # local gradient norms + ϵ_locals = ones(Tr, n) # local gradient norms (drive the adaptive eigensolver tolerance) ϵ_global = maximum(ϵ_locals) + ϵ_convs = ones(Tr, n) # local convergence errors (update magnitude; drive the stop test) + ϵ_conv = maximum(ϵ_convs) ϵ_truncs = zeros(Tr, n) # local truncation error decay_rates = zeros(n) # local observed decay rate of eigensolver fwd, bwd = _sweep_ranges(alg, ψ) LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) + @infov 2 loginit!(log, ϵ_conv, expectation_value(ψ, H, envs)) for iter in 1:(alg.maxiter) @timeit timeroutput "sweep" begin # left-to-right for pos in fwd - ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos], ϵ_convs[pos] = local_update!( pos, Val(:right), ψ, H, alg, envs, @@ -218,7 +228,7 @@ function find_groundstate!( # right-to-left for pos in bwd - ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos], ϵ_convs[pos] = local_update!( pos, Val(:left), ψ, H, alg, envs, @@ -228,25 +238,26 @@ function find_groundstate!( ϵ_global = maximum(ϵ_locals) end end + ϵ_conv = maximum(ϵ_convs) ψ, envs = @timeit timeroutput "finalize" alg.finalize( iter, ψ, H, envs )::Tuple{typeof(ψ), typeof(envs)} - if ϵ_global <= alg.tol + if ϵ_conv <= alg.tol @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) + @infov 2 logfinish!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) break end if iter == alg.maxiter @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) + @warnv 1 logcancel!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) else - @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) + @infov 3 logiter!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) end end end - return ψ, envs, ϵ_global + return ψ, envs, ϵ_conv end function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) From 5c20e8237e1ee0fef45d02d3b0cfbe932052748e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 16:48:05 -0400 Subject: [PATCH 09/18] update default values --- src/utility/dynamictols.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index ade854164..cbe9b5812 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -161,8 +161,8 @@ end function AdaptiveKrylov(; ishermitian::Bool = true, orth::KrylovKit.Orthogonalizer = KrylovKit.KrylovDefaults.orth, - krylovdim_min::Int = 3, krylovdim_max::Int = 8, - iter_min::Int = 1, iter_max::Int = 2, + krylovdim_min::Int = 3, krylovdim_max::Int = 16, + iter_min::Int = 1, iter_max::Int = 1, tol_min::Real = 1.0e-12, tol_max::Real = 1.0e-2, truncation_factor::Real = 1.0e-1, tol_factor::Real = 1.0e-1, eager::Bool = true, verbosity::Int = 0 From 96636604db732357fb1bfcc92058bf5e229e8621 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 12:25:35 -0400 Subject: [PATCH 10/18] Use Galerkin error as DMRG/DMRG2 convergence measure MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Switch the DMRG/DMRG2 stop test from the overlap-based convergence error `ϵ_conv = 1 - |⟨old|new⟩|` to the Galerkin error (already computed each local update as `ϵ_local`/`ϵ_global`), matching VUMPS/iDMRG. Drop the now-unused `ϵ_conv` bookkeeping from both `local_update!` methods and `find_groundstate!`. Perturb the default `FiniteExcited` initial guess with a small per-site random component. Under the Galerkin criterion an unperturbed ground-state guess can already be (near) an eigenvector of the shifted operator, so convergence would trigger on the first sweep and return the ground state instead of climbing to the excited state. The perturbation preserves spaces (symmetry sector) and bond dimensions. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/excitation/dmrgexcitation.jl | 20 +++++++++-- src/algorithms/groundstate/dmrg.jl | 38 ++++++++------------- 2 files changed, 32 insertions(+), 26 deletions(-) diff --git a/src/algorithms/excitation/dmrgexcitation.jl b/src/algorithms/excitation/dmrgexcitation.jl index 8bdb934c6..edeeabb07 100644 --- a/src/algorithms/excitation/dmrgexcitation.jl +++ b/src/algorithms/excitation/dmrgexcitation.jl @@ -18,12 +18,26 @@ $(TYPEDFIELDS) weight::Float64 = 10.0 end +# Default initial guess for the next excited state: the previous state's center tensors with a +# small random perturbation added per site. Without the perturbation the guess can already be +# (near) an eigenvector of the shifted operator, so the Galerkin convergence measure is tiny on +# the first sweep and the optimizer returns immediately instead of climbing to the excited state. +# The perturbation preserves the physical/virtual spaces (hence the symmetry sector) and bond +# dimensions; `FiniteMPS` re-gauges and normalizes the result. +function _perturbed_state(ψ::FiniteMPS; ϵ = 1.0e-2) + return FiniteMPS( + map(1:length(ψ)) do i + A = ψ.AC[i] + noise = randomize!(similar(A)) + return A + (ϵ / norm(noise)) * noise + end + ) +end + function excitations( H::FiniteMPOHamiltonian, alg::FiniteExcited, states::Tuple{T, Vararg{T}}; - init = FiniteMPS( - [ copy(first(states).AC[i]) for i in 1:length(first(states)) ] - ), num = 1 + init = _perturbed_state(first(states)), num = 1 ) where {T <: FiniteMPS} num == 0 && return (scalartype(T)[], T[]) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 5087ca728..c88cdb9bf 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -87,18 +87,15 @@ function local_update!( fixedpoint(H_effective, ac_old, :SR, alg_eigsolve) end - # convergence error: how much the eigensolve moved the center, `1 - |⟨old|new⟩|`. - ϵ_conv = abs(1 - abs(dot(ac_old, AC′)) / (norm(ac_old) * norm(AC′))) - # 3. gauge ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, AC′, alg.alg_gauge; normalize = true) # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) - @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_conv + @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local - return ψ, ϵ_local, ϵ_trunc, decay_rate, ϵ_conv + return ψ, ϵ_local, ϵ_trunc, decay_rate end """ @@ -167,9 +164,6 @@ function local_update!( (newA2center, info) end - # convergence error: how much the eigensolve moved the two-site center, `1 - |⟨old|new⟩|`. - ϵ_conv = abs(1 - abs(dot(ac2, newA2center)) / norm(newA2center)) - # 2. gauge: truncated SVD split back into single-site tensors and install; # the discarded weight is the truncation error ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, newA2center, alg.alg_gauge; normalize = true) @@ -177,9 +171,9 @@ function local_update!( # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) - @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc ϵ_conv + @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc - return ψ, ϵ_local, ϵ_trunc, decay_rate, ϵ_conv + return ψ, ϵ_local, ϵ_trunc, decay_rate end # Per-algorithm sweep geometry: single-site DMRG updates all `length(ψ)` sites (endpoints once, @@ -202,21 +196,19 @@ function find_groundstate!( Tr = real(scalartype(ψ)) n = _num_updates(alg, ψ) - ϵ_locals = ones(Tr, n) # local gradient norms (drive the adaptive eigensolver tolerance) - ϵ_global = maximum(ϵ_locals) - ϵ_convs = ones(Tr, n) # local convergence errors (update magnitude; drive the stop test) - ϵ_conv = maximum(ϵ_convs) + ϵ_locals = ones(Tr, n) # local Galerkin errors (drive both the eigensolver tol and the stop test) + ϵ_global = maximum(ϵ_locals) # sweep-wide Galerkin error, the convergence measure ϵ_truncs = zeros(Tr, n) # local truncation error decay_rates = zeros(n) # local observed decay rate of eigensolver fwd, bwd = _sweep_ranges(alg, ψ) LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ_conv, expectation_value(ψ, H, envs)) + @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) for iter in 1:(alg.maxiter) @timeit timeroutput "sweep" begin # left-to-right for pos in fwd - ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos], ϵ_convs[pos] = + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = local_update!( pos, Val(:right), ψ, H, alg, envs, @@ -228,7 +220,7 @@ function find_groundstate!( # right-to-left for pos in bwd - ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos], ϵ_convs[pos] = + ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] = local_update!( pos, Val(:left), ψ, H, alg, envs, @@ -238,26 +230,26 @@ function find_groundstate!( ϵ_global = maximum(ϵ_locals) end end - ϵ_conv = maximum(ϵ_convs) + ϵ_global = maximum(ϵ_locals) ψ, envs = @timeit timeroutput "finalize" alg.finalize( iter, ψ, H, envs )::Tuple{typeof(ψ), typeof(envs)} - if ϵ_conv <= alg.tol + if ϵ_global <= alg.tol @infov 4 timeroutput - @infov 2 logfinish!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) + @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) break end if iter == alg.maxiter @infov 4 timeroutput - @warnv 1 logcancel!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) + @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) else - @infov 3 logiter!(log, iter, ϵ_conv, expectation_value(ψ, H, envs)) + @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) end end end - return ψ, envs, ϵ_conv + return ψ, envs, ϵ_global end function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) From d0f039748c94de9b0db37603215efee105363e74 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 12:14:23 -0400 Subject: [PATCH 11/18] fix: seed DMRG2 eigensolver with AC, not H*AC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DMRG2 reused `HAC2 = Heff * ac2` (needed for the Galerkin `ϵ_local`) as the eigensolver seed. Seeding with H*AC is power-iteration intuition (it amplifies the largest-magnitude eigenvector), which is wrong for `:SR` (smallest algebraic): it de-weights the ground state, and when the effective ground eigenvalue is 0 it annihilates the ground-state component exactly. Since H preserves the ground state's orthogonal complement, the Krylov solve then converges to the smallest nonzero eigenvalue instead — the wrong state. This bites frustration-free/parent Hamiltonians and any H - E0 shift. Seed with the plain two-site center `ac2` instead, consistent with single-site DMRG's `ac_old`. `HAC2` is still computed for the Galerkin gradient only. Verified on a finite TFIM (L=16, g=1) shifted so E0=0: old seeding gives a non-convergent 4.73 (ϵ_conv≈0.88), fixed seeding gives 5e-9≈0 (ϵ_conv≈2e-15); unshifted DMRG2 still matches single-site DMRG to 1e-14. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/groundstate/dmrg.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index c88cdb9bf..fd8533a1f 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -151,6 +151,12 @@ function local_update!( kind = direction === Val(:right) ? :ACAR : :ALAC ac2 = normalize!(AC2(ψ, pos; kind)) + # `HAC2 = Heff * ac2` only feeds the Galerkin gradient `ϵ_local` below; it must NOT seed the + # eigensolver. Seeding with `Heff * ac2` is power-iteration intuition (amplifies the + # largest-magnitude eigenvector), which is wrong for `:SR` (smallest algebraic): it de-weights + # the ground state and, when the effective ground eigenvalue is 0, annihilates it exactly so + # `:SR` converges to the smallest nonzero eigenvalue. Seed with the plain center `ac2` instead + # (consistent with single-site DMRG's `ac_old`). HAC2 = normalize!(Heff * ac2) AC2′ = copy(HAC2) project_complement!(AC2′, ψ.AL[pos]) @@ -160,7 +166,7 @@ function local_update!( # 1. local two-site update alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) newA2center, info = @timeit timeroutput "AC2_eigsolve" begin - _, newA2center, info = fixedpoint(Heff, HAC2, :SR, alg_eigsolve) + _, newA2center, info = fixedpoint(Heff, ac2, :SR, alg_eigsolve) (newA2center, info) end From 4ee1bc014fa92d7f91f4e4524068febdf9ef0ed5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 16:52:42 -0400 Subject: [PATCH 12/18] rework and unify adaptive solver interface - use kwargs to easily support different metrics and ignore ones that aren't present --- src/algorithms/approximate/idmrg.jl | 4 +- src/algorithms/approximate/vomps.jl | 8 +- src/algorithms/groundstate/dmrg.jl | 4 +- src/algorithms/groundstate/idmrg.jl | 6 +- src/algorithms/groundstate/vumps.jl | 12 +-- src/algorithms/statmech/idmrg.jl | 8 +- src/algorithms/statmech/vomps.jl | 10 +-- src/states/ortho.jl | 4 +- src/utility/dynamictols.jl | 118 +++++++++++++++++----------- 9 files changed, 102 insertions(+), 72 deletions(-) diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index f437a55fa..de53ad77a 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -52,7 +52,7 @@ function approximate!( end # TODO: immediately compute in-place - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) + alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ) ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) copy!(ψ, ψ′) # ensure output destination is unchanged @@ -179,7 +179,7 @@ function approximate!( end # TODO: immediately compute in-place - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) + alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ) ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) copy!(ψ, ψ′) # ensure output destination is unchanged diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 455265256..70f3335b0 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -17,7 +17,7 @@ function approximate( log = IterLog("VOMPS") iter = 0 ϵ = calc_galerkin(mps, toapprox..., envs) - alg_environments = updatetol(alg.alg_environments, iter, ϵ) + alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ) recalculate!(envs, mps, toapprox..., alg_environments) state = VOMPSState(mps, toapprox, envs, iter, ϵ) @@ -65,7 +65,7 @@ end function localupdate_step!( it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler ) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) alg_orth = alg_gauge.alg_orth ACs = similar(state.mps.AC) @@ -89,7 +89,7 @@ end function localupdate_step!( it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler ) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) alg_orth = alg_gauge.alg_orth ACs = similar(state.mps.AC) @@ -119,6 +119,6 @@ function localupdate_step!( end function envs_step!(it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, mps) - alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) + alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ) return recalculate!(state.envs, mps, state.operator..., alg_environments) end diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index fd8533a1f..6cb68d40d 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -80,7 +80,7 @@ function local_update!( @timeit timeroutput "expand" changebond!(site, direction, ψ, O, alg.alg_expand, envs) # 2. local update - alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) + alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc) ac_old = ψ.AC[site] λ, AC′, info = @timeit timeroutput "AC_eigsolve" begin H_effective = AC_hamiltonian(site, ψ, O, ψ, envs) @@ -164,7 +164,7 @@ function local_update!( ϵ_local = norm(AC2′) # 1. local two-site update - alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) + alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc) newA2center, info = @timeit timeroutput "AC2_eigsolve" begin _, newA2center, info = fixedpoint(Heff, ac2, :SR, alg_eigsolve) (newA2center, info) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index ed83e664a..66143b3bf 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -110,7 +110,7 @@ function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, @infov 3 logiter!(log, it.iter, ϵ, ΔE) end - alg_gauge = updatetol(alg.alg_gauge, it.state.iter, it.state.ϵ) + alg_gauge = adapt_solver(alg.alg_gauge; iter = it.state.iter, g_global = it.state.ϵ) ψ′ = InfiniteMPS(it.state.mps.AR; alg_gauge.tol, alg_gauge.maxiter) envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′) return ψ′, envs, it.state.ϵ @@ -151,7 +151,7 @@ end function localupdate_step!( it::IterativeSolver{<:IDMRG}, state ) - alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) return _localupdate_sweep_idmrg!( state.mps, state.operator, state.envs, alg_eigsolve, state.timeroutput, ) @@ -160,7 +160,7 @@ end function localupdate_step!( it::IterativeSolver{<:IDMRG2}, state ) - alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) return _localupdate_sweep_idmrg2!( state.mps, state.operator, state.envs, alg_eigsolve, it.trscheme, it.alg_svd, state.timeroutput, diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index aa1c413cb..1f1cf51f7 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -63,7 +63,7 @@ function dominant_eigsolve( mps = copy(mps) ϵ = calc_galerkin(mps, operator, mps, envs) - alg_environments = updatetol(alg.alg_environments, iter, ϵ) + alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ) recalculate!(envs, mps, operator, mps, alg_environments; timeroutput) state = VUMPSState(mps, operator, envs, iter, ϵ, which, timeroutput) @@ -118,8 +118,8 @@ end function localupdate_step!( it::IterativeSolver{<:VUMPS}, state, scheduler = Defaults.scheduler[] ) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) - alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) alg_orth = alg_gauge.alg_orth mps = state.mps @@ -186,7 +186,7 @@ function _localupdate_vumps_step!( end function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractVector) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) mps = gaugefix!( state.mps, ACs, state.mps.C[end]; order = :R, timeroutput = state.timeroutput, alg_gauge..., @@ -195,11 +195,11 @@ function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractVector) return mps end function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractMatrix) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge.tol, alg_gauge.maxiter, alg_gauge.alg_orth) end function envs_step!(it::IterativeSolver{<:VUMPS}, state, mps) - alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) + alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ) return recalculate!(state.envs, mps, state.operator, mps, alg_environments; state.timeroutput) end diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 614e9f6bf..208da3664 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -8,7 +8,7 @@ function leading_boundary( LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ, expectation_value(ψ, operator, envs)) for outer iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) + alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ) C_current = ψ.C[:, 0] # left to right sweep @@ -54,7 +54,7 @@ function leading_boundary( end end - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) + alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ) ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) @@ -72,7 +72,7 @@ function leading_boundary( LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) for outer iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) + alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ) C_current = ψ.C[:, 0] # sweep from left to right @@ -189,7 +189,7 @@ function leading_boundary( end end - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) + alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ) ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 9295f2b28..935d9111f 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -54,7 +54,7 @@ function dominant_eigsolve( log = IterLog("VOMPS") iter = 0 ϵ = calc_galerkin(mps, operator, mps, envs) - alg_environments = updatetol(alg.alg_environments, iter, ϵ) + alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ) recalculate!(envs, mps, operator, mps, alg_environments) state = VOMPSState(mps, operator, envs, iter, ϵ) @@ -100,7 +100,7 @@ end function localupdate_step!( it::IterativeSolver{<:VOMPS}, state, scheduler = Defaults.scheduler[] ) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) alg_orth = alg_gauge.alg_orth mps = state.mps ACs = similar(mps.AC) @@ -135,15 +135,15 @@ function _localupdate_vomps_step!( end function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractVector) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) return InfiniteMPS(ACs, state.mps.C[end]; alg_gauge...) end function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractMatrix) - alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ) return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge...) end function envs_step!(it::IterativeSolver{<:VOMPS}, state, mps) - alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) + alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ) return recalculate!(state.envs, mps, state.operator, mps, alg_environments) end diff --git a/src/states/ortho.jl b/src/states/ortho.jl index 0b4cbd6ca..d7901de08 100644 --- a/src/states/ortho.jl +++ b/src/states/ortho.jl @@ -246,7 +246,7 @@ end function gauge_eigsolve_step!(it::IterativeSolver{LeftCanonical}, state) (; AL, C, A, iter, ϵ) = state if iter ≥ it.eig_miniter - alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2) + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = 1, g_global = ϵ^2) _, vec = fixedpoint(flip(TransferMatrix(A, AL)), C[end], :LM, alg_eigsolve) _, C[end] = left_orth!(vec; alg = it.alg_orth) end @@ -308,7 +308,7 @@ end function gauge_eigsolve_step!(it::IterativeSolver{RightCanonical}, state) (; AR, C, A, iter, ϵ) = state if iter ≥ it.eig_miniter - alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2) + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = 1, g_global = ϵ^2) _, vec = fixedpoint(TransferMatrix(A, AR), C[end], :LM, alg_eigsolve) C[end], _ = right_orth!(vec; alg = it.alg_orth) end diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index cbe9b5812..c16e4f990 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -6,16 +6,33 @@ using DocStringExtensions using MatrixAlgebraKit: DefaultAlgorithm using KrylovKit: KrylovKit, Lanczos, Arnoldi -export updatetol, DynamicTol -export AdaptiveKrylov, instantiate_algorithm +export adapt_solver, DynamicTol, AdaptiveKrylov +# deprecated, kept for backwards compatibility (see the bottom of this file) +export updatetol + +# ============================================================ +# Unified solver-adaptation interface +# ============================================================ @doc """ - updatetol(alg, iter, ϵ) + adapt_solver(alg; iter, decay_rate, g_local, g_global, eps_trunc) + +Resolve a (possibly adaptive) solver specification `alg` into a concrete solver, using +whichever adaptation signals are supplied as keyword arguments. Numeric signals default to +`0.0` (and `iter` to `1`) — i.e. "no contribution" — so each caller passes only what it has in +scope: -Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`. -""" updatetol + - `iter` — outer iteration count (drives the `1/√iter` damping) + - `decay_rate` — measured per-matvec contraction factor of the previous solve + - `g_local` — local (per-bond) gradient / Galerkin norm + - `g_global` — global (sweep-wide) gradient / convergence-error scalar + - `eps_trunc` — local truncation error -updatetol(alg, iter::Integer, ϵ::Real) = alg +The generic fallback returns `alg` unchanged (for plain `Lanczos`/`Arnoldi`/… solvers, which +are not adapted). `DynamicTol` and `AdaptiveKrylov` provide the adaptive implementations. +""" adapt_solver + +adapt_solver(alg; kwargs...) = alg # Wrapper for dynamic tolerance adjustment # ---------------------------------------- @@ -23,13 +40,15 @@ updatetol(alg, iter::Integer, ϵ::Real) = alg """ $(TYPEDEF) -Algorithm wrapper with dynamically adjusted tolerances. +Algorithm wrapper with dynamically adjusted tolerances. Only the wrapped solver's tolerance is +retuned; its Krylov budget (if any) is left fixed — this is the simpler counterpart to +[`AdaptiveKrylov`](@ref). ## Fields $(TYPEDFIELDS) -See also [`updatetol`](@ref). +See also [`adapt_solver`](@ref). """ struct DynamicTol{A} <: Algorithm "parent algorithm" @@ -41,33 +60,48 @@ struct DynamicTol{A} <: Algorithm "maximal value of the dynamic tolerance" tol_max::Float64 - "tolerance factor for updating relative to current algorithm error" + "tolerance factor for updating relative to the current (global) gradient norm" tol_factor::Float64 + "factor on the local truncation error, sets an inner-solve tolerance floor (`0` ⇒ truncation-agnostic)" + truncation_factor::Float64 + function DynamicTol( - alg::A, tol_min::Real, tol_max::Real, tol_factor::Real + alg::A, tol_min::Real, tol_max::Real, tol_factor::Real, truncation_factor::Real = 0.0 ) where {A} 0 <= tol_min <= tol_max || throw(ArgumentError("tol_min must be between 0 and tol_max")) - return new{A}(alg, tol_min, tol_max, tol_factor) + truncation_factor >= 0 || + throw(ArgumentError("truncation_factor must be non-negative, got $truncation_factor")) + return new{A}(alg, tol_min, tol_max, tol_factor, truncation_factor) end end -function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1) - return DynamicTol(alg, tol_min, tol_max, tol_factor) +function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1, truncation_factor = 0.0) + return DynamicTol(alg, tol_min, tol_max, tol_factor, truncation_factor) end """ - updatetol(alg::DynamicTol, iter, ϵ) + adapt_solver(alg::DynamicTol; iter, g_global, eps_trunc, ...) + +Tighten only the wrapped solver's tolerance (its Krylov budget, if any, is left fixed). The +target combines the truncation-error floor with the global-gradient-driven convergence target, +optionally damped by the iteration count: + + tol = clamp(max(truncation_factor·eps_trunc, tol_factor·g_global) / √iter, tol_min, tol_max) -Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`, -where the new tolerance is given by - - new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) +Per-bond callers (finite DMRG) supply `g_global`/`eps_trunc` and leave `iter = 1` (no damping); +per-sweep/global callers (VUMPS/iDMRG/…) supply the global error as `g_global` together with +`iter`, and `eps_trunc` defaults to `0`. """ -function updatetol(alg::DynamicTol, iter::Integer, ϵ::Real) - iter = max(iter, one(iter)) - new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max) - return _updatetol(alg.alg, new_tol) +function adapt_solver( + alg::DynamicTol; + iter::Integer = 1, g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs... + ) + trunc_tol = alg.truncation_factor * eps_trunc + conv_tol = alg.tol_factor * g_global + tol = clamp(max(trunc_tol, conv_tol) / sqrt(max(iter, 1)), alg.tol_min, alg.tol_max) + + return _updatetol(alg.alg, tol) end # default implementation with Accessors.jl, but can be hooked into @@ -105,7 +139,7 @@ while avoiding stagnation for gapless ones. $(TYPEDFIELDS) -See also [`instantiate_algorithm`](@ref). +See also [`adapt_solver`](@ref). """ struct AdaptiveKrylov{T, O <: KrylovKit.Orthogonalizer} <: Algorithm "orthogonalizer passed to the instantiated `Lanczos`/`Arnoldi`" @@ -174,32 +208,22 @@ function AdaptiveKrylov(; end """ - instantiate_algorithm(alg, decay_rate, g_local, g_global, eps_trunc) + adapt_solver(alg::AdaptiveKrylov; decay_rate, g_local, g_global, eps_trunc, ...) -Turn a (possibly adaptive) local-eigensolver specification `alg` into a concrete KrylovKit -algorithm for the current site, using the measured `decay_rate` of the previous solve, the -local gradient norm `g_local`, the global gradient norm `g_global` and the local truncation -error `eps_trunc`. The generic fallback returns `alg` unchanged (for plain `Lanczos`/`Arnoldi`). +Build a concrete `Lanczos`/`Arnoldi` for the current site. The target tolerance is the same as +[`DynamicTol`](@ref)'s per-bond tolerance, but the Krylov `krylovdim`/`maxiter` are additionally +predicted from the measured `decay_rate` and the local gradient norm `g_local`. A `decay_rate` +of `0` (the default) signals a cold start and falls back to the minimal budget. """ -instantiate_algorithm(alg, args...) = alg - -# a `DynamicTol` wrapper is resolved to a concrete solver by tightening its tolerance with the -# global gradient norm (the per-sweep, tol-only "legacy" behaviour, now driven per site). -function instantiate_algorithm( - alg::DynamicTol, decay_rate::Real, g_local::Real, g_global::Real, eps_trunc::Real - ) - tol = clamp(g_global * alg.tol_factor, alg.tol_min, alg.tol_max) - return _updatetol(alg.alg, tol) -end - -function instantiate_algorithm( - alg::AdaptiveKrylov{T}, decay_rate::Real, g_local::Real, g_global::Real, eps_trunc::Real +function adapt_solver( + alg::AdaptiveKrylov{T}; decay_rate::Real = 0.0, g_local::Real = 0.0, + g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs... ) where {T} # 1. target tolerance: inexact inner solve depending on outer convergence, # never below the truncation error we already incur. - trunc_tol = alg.truncation_factor * eps_trunc - conv_tol = alg.tol_factor * g_global - tol = clamp(max(trunc_tol, conv_tol), alg.tol_min, alg.tol_max) + tol = clamp( + max(alg.truncation_factor * eps_trunc, alg.tol_factor * g_global), alg.tol_min, alg.tol_max + ) # the measured decay rate must be a genuine contraction in (0, 1); a non-positive or ≥ 1 # rate signals an uninitialized or hard/stalling bond → fall back to the largest budget. @@ -224,4 +248,10 @@ function instantiate_algorithm( return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity) end +# ============================================================ +# Deprecated entry points (subsumed by `adapt_solver`) +# ============================================================ + +Base.@deprecate updatetol(alg, iter::Integer, ϵ::Real) adapt_solver(alg; iter = iter, g_global = ϵ) + end From fd9ba087ab705eee8fcca2cced41bc5f91db4756 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 16:54:23 -0400 Subject: [PATCH 13/18] update defaults once more --- src/utility/defaults.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index 11cf238eb..8af641bcc 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -54,14 +54,13 @@ end function alg_eigsolve(; ishermitian = true, tol = tol, maxiter = maxiter, verbosity = 0, eager = true, krylovdim = krylovdim, dynamic_tols = dynamic_tols, tol_min = tol_min, - tol_max = tol_max, tol_factor = eigs_tolfactor, adaptive = false, adaptive_kwargs... + tol_max = tol_max, tol_factor = eigs_tolfactor, truncation_factor = 0.0, + adaptive = false, adaptive_kwargs... ) - # the per-bond adaptive controller (`AdaptiveKrylov`) retunes krylovdim/maxiter/tol per site - # and supersedes the (global, tol-only) `DynamicTol` wrapper; it is the default for `DMRG`. adaptive && return AdaptiveKrylov(; ishermitian, adaptive_kwargs...) alg = ishermitian ? Lanczos(; tol, maxiter, eager, krylovdim, verbosity) : Arnoldi(; tol, maxiter, eager, krylovdim, verbosity) - return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg + return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor, truncation_factor) : alg end function alg_environments(; From 6d2b2f8141bf4c0485e18785253d57f4d14921b7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 17:21:06 -0400 Subject: [PATCH 14/18] revert DynamicTol behavior --- src/utility/defaults.jl | 5 ++--- src/utility/dynamictols.jl | 38 ++++++++++---------------------------- 2 files changed, 12 insertions(+), 31 deletions(-) diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index 8af641bcc..dbc221420 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -54,13 +54,12 @@ end function alg_eigsolve(; ishermitian = true, tol = tol, maxiter = maxiter, verbosity = 0, eager = true, krylovdim = krylovdim, dynamic_tols = dynamic_tols, tol_min = tol_min, - tol_max = tol_max, tol_factor = eigs_tolfactor, truncation_factor = 0.0, - adaptive = false, adaptive_kwargs... + tol_max = tol_max, tol_factor = eigs_tolfactor, adaptive = false, adaptive_kwargs... ) adaptive && return AdaptiveKrylov(; ishermitian, adaptive_kwargs...) alg = ishermitian ? Lanczos(; tol, maxiter, eager, krylovdim, verbosity) : Arnoldi(; tol, maxiter, eager, krylovdim, verbosity) - return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor, truncation_factor) : alg + return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end function alg_environments(; diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index c16e4f990..6c28a95ce 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -63,44 +63,26 @@ struct DynamicTol{A} <: Algorithm "tolerance factor for updating relative to the current (global) gradient norm" tol_factor::Float64 - "factor on the local truncation error, sets an inner-solve tolerance floor (`0` ⇒ truncation-agnostic)" - truncation_factor::Float64 - - function DynamicTol( - alg::A, tol_min::Real, tol_max::Real, tol_factor::Real, truncation_factor::Real = 0.0 - ) where {A} + function DynamicTol(alg::A, tol_min::Real, tol_max::Real, tol_factor::Real) where {A} 0 <= tol_min <= tol_max || throw(ArgumentError("tol_min must be between 0 and tol_max")) - truncation_factor >= 0 || - throw(ArgumentError("truncation_factor must be non-negative, got $truncation_factor")) - return new{A}(alg, tol_min, tol_max, tol_factor, truncation_factor) + return new{A}(alg, tol_min, tol_max, tol_factor) end end -function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1, truncation_factor = 0.0) - return DynamicTol(alg, tol_min, tol_max, tol_factor, truncation_factor) +function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1) + return DynamicTol(alg, tol_min, tol_max, tol_factor) end """ - adapt_solver(alg::DynamicTol; iter, g_global, eps_trunc, ...) + adapt_solver(alg::DynamicTol; iter, g_global, ...) -Tighten only the wrapped solver's tolerance (its Krylov budget, if any, is left fixed). The -target combines the truncation-error floor with the global-gradient-driven convergence target, -optionally damped by the iteration count: +Tighten only the wrapped solver's tolerance (its Krylov budget, if any, is left fixed), from the +global gradient / convergence-error scalar `g_global`, damped by the iteration count: - tol = clamp(max(truncation_factor·eps_trunc, tol_factor·g_global) / √iter, tol_min, tol_max) - -Per-bond callers (finite DMRG) supply `g_global`/`eps_trunc` and leave `iter = 1` (no damping); -per-sweep/global callers (VUMPS/iDMRG/…) supply the global error as `g_global` together with -`iter`, and `eps_trunc` defaults to `0`. + tol = clamp(tol_factor·g_global / √iter, tol_min, tol_max) """ -function adapt_solver( - alg::DynamicTol; - iter::Integer = 1, g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs... - ) - trunc_tol = alg.truncation_factor * eps_trunc - conv_tol = alg.tol_factor * g_global - tol = clamp(max(trunc_tol, conv_tol) / sqrt(max(iter, 1)), alg.tol_min, alg.tol_max) - +function adapt_solver(alg::DynamicTol; iter::Integer = 1, g_global::Real = 0.0, kwargs...) + tol = clamp(alg.tol_factor * g_global / sqrt(max(iter, 1)), alg.tol_min, alg.tol_max) return _updatetol(alg.alg, tol) end From f267171394792529a05ace45825f418485cde853 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 16 Jul 2026 17:47:49 -0400 Subject: [PATCH 15/18] guard two-site gauge complexification on state scalartype `gauge2!` unconditionally promoted the SVD center matrix to complex (`complex(c)`), which for a real-valued state forced a complexify then coerce-back-to-real roundtrip when installing the tensors. Only complexify when the state's scalartype is itself complex. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/states/orthoview.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 910507b7c..efa394b8a 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -326,7 +326,9 @@ function gauge2!( ) where {Dir} al, c, ar, ϵ = svd_trunc!(AC2, alg) normalize && normalize!(c) - C = complex(c) + # the SVD center `c` (singular values) is real-valued; promote it to the state's scalartype so + # the installed tensors stay consistent, without needlessly forcing a real-valued state complex. + C = scalartype(ψ) <: Complex ? complex(c) : c if Dir === :right ψ.AC[pos] = (al, C) ψ.AC[pos + 1] = (C, _transpose_front(ar)) From e7c9fdb0a7ad458cf1739284df328a652569c8fc Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 17 Jul 2026 10:26:11 -0400 Subject: [PATCH 16/18] fix DMRG2 convergence measure; add truncation-aware stop MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The two-site Galerkin error was projected onto BOTH the left and right null spaces (the (null,null) bond-expansion corner), which is exactly what the truncated SVD discards. It therefore collapsed to the truncation floor after a single local solve and spuriously signalled convergence after one sweep, leaving DMRG2 under-converged (e.g. the multifusion Ising ground state stopped 2.6e-6 above exact with variance 3e-5 after one sweep). Project onto the LEFT null space of `AL[pos]` only, matching single-site `calc_galerkin`, so the measure faithfully tracks convergence. The faithful measure floors at the truncation error (the gradient of a truncated state toward the untruncated optimum), so `ϵ_global <= tol` is unreachable under truncation. Converge instead at `max(tol, maxϵ_trunc)` — once the gradient reaches the level the truncation allows. With no truncation (single-site / QR gauge, `ϵ_trunc = 0`) this reduces to the plain `<= tol`. Note: penalized runs (`FiniteExcited`) can still floor above the truncation and hit maxiter; a more general stop criterion is left to a follow-up. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/groundstate/dmrg.jl | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 6cb68d40d..c06fd87d5 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -151,16 +151,20 @@ function local_update!( kind = direction === Val(:right) ? :ACAR : :ALAC ac2 = normalize!(AC2(ψ, pos; kind)) - # `HAC2 = Heff * ac2` only feeds the Galerkin gradient `ϵ_local` below; it must NOT seed the - # eigensolver. Seeding with `Heff * ac2` is power-iteration intuition (amplifies the - # largest-magnitude eigenvector), which is wrong for `:SR` (smallest algebraic): it de-weights - # the ground state and, when the effective ground eigenvalue is 0, annihilates it exactly so - # `:SR` converges to the smallest nonzero eigenvalue. Seed with the plain center `ac2` instead - # (consistent with single-site DMRG's `ac_old`). - HAC2 = normalize!(Heff * ac2) - AC2′ = copy(HAC2) + # Two-site Galerkin gradient `ϵ_local`, the direct analogue of single-site `calc_galerkin`: the + # part of `Heff * ac2` that points off the current left-canonical manifold, obtained by + # projecting onto the LEFT null space of `AL[pos]` only. Projecting onto BOTH the left and right + # null spaces would keep only the (null, null) bond-expansion corner — the directions a + # truncated SVD discards — which collapses to the truncation floor after each local solve and + # so spuriously signals convergence after a single sweep. + # + # `Heff * ac2` only feeds this gradient; it must NOT seed the eigensolver. Seeding with + # `Heff * ac2` is power-iteration intuition (amplifies the largest-magnitude eigenvector), which + # is wrong for `:SR` (smallest algebraic): it de-weights the ground state and, when the + # effective ground eigenvalue is 0, annihilates it exactly so `:SR` converges to the smallest + # nonzero eigenvalue. Seed with the plain center `ac2` instead (as single-site DMRG's `ac_old`). + AC2′ = normalize!(Heff * ac2) project_complement!(AC2′, ψ.AL[pos]) - project_complement_right!(AC2′, _transpose_tail(ψ.AR[pos + 1])) ϵ_local = norm(AC2′) # 1. local two-site update @@ -242,7 +246,11 @@ function find_groundstate!( iter, ψ, H, envs )::Tuple{typeof(ψ), typeof(envs)} - if ϵ_global <= alg.tol + # Truncation-aware convergence: the Galerkin gradient cannot drop below the level set by + # the discarded weight, so a truncating scheme converges once `ϵ_global` reaches the + # truncation error rather than the (unreachable) bare `tol`. With no truncation + # (`ϵ_truncs .= 0`, e.g. single-site/QR gauge) this reduces to the plain `ϵ_global ≤ tol`. + if ϵ_global <= max(alg.tol, maximum(ϵ_truncs)) @infov 4 timeroutput @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs)) break From 4b58685c24029cd3552b11cb6e86d1c8eedc0b7d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 17 Jul 2026 10:26:11 -0400 Subject: [PATCH 17/18] restore updatetol docstring `updatetol` kept a docstring before the `adapt_solver` refactor; restore one for the deprecated shim so the exported name stays documented. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/utility/dynamictols.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index 6c28a95ce..e52f1531c 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -236,4 +236,11 @@ end Base.@deprecate updatetol(alg, iter::Integer, ϵ::Real) adapt_solver(alg; iter = iter, g_global = ϵ) +@doc """ + updatetol(alg, iter, ϵ) + +Deprecated: superseded by [`adapt_solver`](@ref). Equivalent to +`adapt_solver(alg; iter, g_global = ϵ)`. +""" updatetol + end From e3e701260a5148678ff9e2c40911e5d84306af94 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 17 Jul 2026 13:55:14 -0400 Subject: [PATCH 18/18] cleanup --- src/algorithms/excitation/dmrgexcitation.jl | 13 ++++--------- src/algorithms/groundstate/dmrg.jl | 14 +------------- src/utility/dynamictols.jl | 2 +- 3 files changed, 6 insertions(+), 23 deletions(-) diff --git a/src/algorithms/excitation/dmrgexcitation.jl b/src/algorithms/excitation/dmrgexcitation.jl index edeeabb07..72c47a401 100644 --- a/src/algorithms/excitation/dmrgexcitation.jl +++ b/src/algorithms/excitation/dmrgexcitation.jl @@ -18,18 +18,13 @@ $(TYPEDFIELDS) weight::Float64 = 10.0 end -# Default initial guess for the next excited state: the previous state's center tensors with a -# small random perturbation added per site. Without the perturbation the guess can already be -# (near) an eigenvector of the shifted operator, so the Galerkin convergence measure is tiny on -# the first sweep and the optimizer returns immediately instead of climbing to the excited state. -# The perturbation preserves the physical/virtual spaces (hence the symmetry sector) and bond -# dimensions; `FiniteMPS` re-gauges and normalizes the result. -function _perturbed_state(ψ::FiniteMPS; ϵ = 1.0e-2) +# Initialize excited state by perturbing the current eigenvector to avoid local minima +function _perturbed_state(ψ::FiniteMPS; atol = 1.0e-2) return FiniteMPS( map(1:length(ψ)) do i - A = ψ.AC[i] + A = i == length(ψ) ? ψ.AC[i] : ψ.AL[i] noise = randomize!(similar(A)) - return A + (ϵ / norm(noise)) * noise + return add!(noise, A, 1, atol / norm(noise)) end ) end diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index c06fd87d5..de8ac4bd8 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -150,19 +150,7 @@ function local_update!( Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs) kind = direction === Val(:right) ? :ACAR : :ALAC - ac2 = normalize!(AC2(ψ, pos; kind)) - # Two-site Galerkin gradient `ϵ_local`, the direct analogue of single-site `calc_galerkin`: the - # part of `Heff * ac2` that points off the current left-canonical manifold, obtained by - # projecting onto the LEFT null space of `AL[pos]` only. Projecting onto BOTH the left and right - # null spaces would keep only the (null, null) bond-expansion corner — the directions a - # truncated SVD discards — which collapses to the truncation floor after each local solve and - # so spuriously signals convergence after a single sweep. - # - # `Heff * ac2` only feeds this gradient; it must NOT seed the eigensolver. Seeding with - # `Heff * ac2` is power-iteration intuition (amplifies the largest-magnitude eigenvector), which - # is wrong for `:SR` (smallest algebraic): it de-weights the ground state and, when the - # effective ground eigenvalue is 0, annihilates it exactly so `:SR` converges to the smallest - # nonzero eigenvalue. Seed with the plain center `ac2` instead (as single-site DMRG's `ac_old`). + ac2 = AC2(ψ, pos; kind) AC2′ = normalize!(Heff * ac2) project_complement!(AC2′, ψ.AL[pos]) ϵ_local = norm(AC2′) diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index e52f1531c..fe52d3188 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -177,7 +177,7 @@ end function AdaptiveKrylov(; ishermitian::Bool = true, orth::KrylovKit.Orthogonalizer = KrylovKit.KrylovDefaults.orth, - krylovdim_min::Int = 3, krylovdim_max::Int = 16, + krylovdim_min::Int = 3, krylovdim_max::Int = 20, iter_min::Int = 1, iter_max::Int = 1, tol_min::Real = 1.0e-12, tol_max::Real = 1.0e-2, truncation_factor::Real = 1.0e-1, tol_factor::Real = 1.0e-1,