From c55f6d1c789518b48e7eb31fe1b3a7440ea6f92e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 26 Jul 2026 13:05:33 -0400 Subject: [PATCH 1/7] fix: operator-free two-site MPO derivative contraction The `MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor}` action (the operator-free / overlap two-site derivative) was missing a semicolon in the `rightenv` index, declaring a `(2, 0)` index pattern for a `(1, 1)` bond tensor. This threw a `TensorOperations.IndexError` for the planar bond-tensor case (e.g. `approximate(x, b, DMRG2())` on planar states); the ComplexSpace path happened to avoid the broken method. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/derivatives/mpo_derivatives.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 5bc4bc943..55d1be42c 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -90,7 +90,7 @@ end function (h::MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor})( x::MPOTensor ) - @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2 -3] + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2; -3] return y isa AbstractBlockTensorMap ? only(y) : y end function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})( From 6d7be143bfebbe6569079915ec2be2751198aceb Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 26 Jul 2026 13:05:46 -0400 Subject: [PATCH 2/7] refactor: qualify internal linsolve calls as KrylovKit.linsolve Use the fully-qualified `KrylovKit.linsolve` at the internal transfer-matrix / quasiparticle / fidelity call sites, so that the top-level `linsolve` name can be owned by MPSKit for the upcoming MPS linear-solver category without shadowing these tensor-level solves. Behavior-preserving. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/excitation/exci_transfer_system.jl | 4 ++-- src/algorithms/fidelity_susceptibility.jl | 2 +- src/algorithms/propagator/corvector.jl | 4 ++-- src/environments/infinite_envs.jl | 8 ++++---- src/environments/qp_envs.jl | 4 ++-- src/states/quasiparticle_state.jl | 4 ++-- 6 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/algorithms/excitation/exci_transfer_system.jl b/src/algorithms/excitation/exci_transfer_system.jl index c2b4b217d..c3210be2a 100644 --- a/src/algorithms/excitation/exci_transfer_system.jl +++ b/src/algorithms/excitation/exci_transfer_system.jl @@ -37,7 +37,7 @@ function left_excitation_transfer_system( ) end - found[i], convhist = linsolve( + found[i], convhist = KrylovKit.linsolve( flip(T), found[i], found[i], solver, 1, -cis(-mom * len) ) convhist.converged == 0 && @@ -87,7 +87,7 @@ function right_excitation_transfer_system( ) end - found[i], convhist = linsolve( + found[i], convhist = KrylovKit.linsolve( tm, found[i], found[i], solver, 1, -cis(mom * len) ) convhist.converged < 1 && diff --git a/src/algorithms/fidelity_susceptibility.jl b/src/algorithms/fidelity_susceptibility.jl index 55dbac066..bbb89e4d0 100644 --- a/src/algorithms/fidelity_susceptibility.jl +++ b/src/algorithms/fidelity_susceptibility.jl @@ -29,7 +29,7 @@ function fidelity_susceptibility( @plansor Tos[i][-1 -2; -3 -4] := temp[-1 -2; -4] * help[-3] end - vec, convhist = linsolve(Tos, Tos, GMRES(; maxiter = maxiter, tol = tol)) do x + vec, convhist = KrylovKit.linsolve(Tos, Tos, GMRES(; maxiter = maxiter, tol = tol)) do x return effective_excitation_hamiltonian(H₀, x, environments(x, H₀, x; lenvs = henvs)) end convhist.converged == 0 && @warn "failed to converge: normres = $(convhist.normres)" diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 89466e7e0..bf7057b6d 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -87,7 +87,7 @@ function propagator( H_AC = AC_hamiltonian(i, init, H, init, h_envs) AC = init.AC[i] - AC′, convhist = linsolve(H_AC, -tos, AC, alg.solver, -z, one(z)) + AC′, convhist = KrylovKit.linsolve(H_AC, -tos, AC, alg.solver, -z, one(z)) ϵ = max(ϵ, norm(AC′ - AC)) init.AC[i] = AC′ @@ -163,7 +163,7 @@ function propagator( H1_AC = AC_hamiltonian(i, init, H, init, envs1) H2_AC = AC_hamiltonian(i, init, H2, init, envs2) H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1)) - AC′, convhist = linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1) + AC′, convhist = KrylovKit.linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1) ϵ = max(ϵ, norm(AC′ - init.AC[i])) init.AC[i] = AC′ diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl index 25285068e..33c7f2bc0 100644 --- a/src/environments/infinite_envs.jl +++ b/src/environments/infinite_envs.jl @@ -233,7 +233,7 @@ function compute_leftenvs!( if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization T = regularize(TransferMatrix(above.AL, below.AL), ρ_left, ρ_right) - GLs[1][i], convhist = linsolve(flip(T), GLs[1][i], prev, alg, 1, -1) + GLs[1][i], convhist = KrylovKit.linsolve(flip(T), GLs[1][i], prev, alg, 1, -1) convhist.converged == 0 && @warn "GL$i failed to converge: normres = $(convhist.normres)" @@ -249,7 +249,7 @@ function compute_leftenvs!( if !isemptylevel(operator, i) diag = map(h -> h[i, 1, 1, i], operator[:]) T = TransferMatrix(above.AL, diag, below.AL) - GLs[1][i], convhist = linsolve(flip(T), GLs[1][i], prev, alg, 1, -1) + GLs[1][i], convhist = KrylovKit.linsolve(flip(T), GLs[1][i], prev, alg, 1, -1) convhist.converged == 0 && @warn "GL$i failed to converge: normres = $(convhist.normres)" end @@ -304,7 +304,7 @@ function compute_rightenvs!( if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization # subtract fixpoints T = regularize(TransferMatrix(above.AR, below.AR), ρ_left, ρ_right) - GRs[end][i], convhist = linsolve(T, GRs[end][i], prev, alg, 1, -1) + GRs[end][i], convhist = KrylovKit.linsolve(T, GRs[end][i], prev, alg, 1, -1) convhist.converged == 0 && @warn "GR$i failed to converge: normres = $(convhist.normres)" @@ -319,7 +319,7 @@ function compute_rightenvs!( if !isemptylevel(operator, i) diag = map(b -> b[i, 1, 1, i], operator[:]) T = TransferMatrix(above.AR, diag, below.AR) - GRs[end][i], convhist = linsolve(T, GRs[end][i], prev, alg, 1, -1) + GRs[end][i], convhist = KrylovKit.linsolve(T, GRs[end][i], prev, alg, 1, -1) convhist.converged == 0 && @warn "GR$i failed to converge: normres = $(convhist.normres)" end diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index 9429da3a6..da99e3c13 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -226,7 +226,7 @@ function environments(exci::InfiniteQP, O::InfiniteMPO, above, alg; lenvs, renvs T_LR = regularize(T_LR, lvec, rvec) end - GBL[1], convhist = linsolve( + GBL[1], convhist = KrylovKit.linsolve( flip(T_RL), gbl, gbl, solver, 1, -cis(-length(exci) * exci.momentum) * prod(left_regularization) ) @@ -234,7 +234,7 @@ function environments(exci::InfiniteQP, O::InfiniteMPO, above, alg; lenvs, renvs convhist.converged == 0 && @warn "GBL failed to converge: normres = $(convhist.normres)" - GBR[end], convhist = linsolve( + GBR[end], convhist = KrylovKit.linsolve( T_LR, gbr, gbr, GMRES(), 1, -cis(length(exci) * exci.momentum) * prod(right_regularization) ) diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl index 6934e2b22..e8e2d2d33 100644 --- a/src/states/quasiparticle_state.jl +++ b/src/states/quasiparticle_state.jl @@ -192,7 +192,7 @@ function Base.convert( tm = regularize(tm, l_LR(input.right_gs), r_LR(input.right_gs)) end - rBE, convhist = linsolve( + rBE, convhist = KrylovKit.linsolve( tm, rBs[1], rBs[1], GMRES(), 1, -exp(1im * input.momentum * len) ) convhist.converged == 0 && @warn "failed to converge: normres = $(convhist.normres)" @@ -239,7 +239,7 @@ function Base.convert( tm = regularize(tm, l_RL(input.right_gs), r_RL(input.right_gs)) end - lBE, convhist = linsolve( + lBE, convhist = KrylovKit.linsolve( flip(tm), lBs[end], lBs[end], GMRES(), 1, -1 / exp(1im * input.momentum * len) ) convhist.converged == 0 && @warn "failed to converge: normres = $(convhist.normres)" From b4b0268f1233480ab9281ffe5d15b8a6dabdd6cb Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 26 Jul 2026 13:05:56 -0400 Subject: [PATCH 3/7] feat: optional truncation-error floor for DynamicTol Add a `truncation_factor` field to `DynamicTol` and let `adapt_solver` floor the retuned tolerance at `truncation_factor * eps_trunc`, i.e. `tol = clamp(max(truncation_factor*eps_trunc, tol_factor*g_global/sqrt(iter)), tol_min, tol_max)`. The default `truncation_factor = 0` keeps existing behavior (DMRG's gauge / eigensolver uses are unchanged); a positive value is used by truncating two-site sweeps to avoid solving a bond far below what the SVD truncation will discard. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/utility/dynamictols.jl | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index 5ef4cc9f1..8c39e44cd 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -65,26 +65,40 @@ struct DynamicTol{A} <: Algorithm "tolerance factor for updating relative to the current (global) gradient norm" tol_factor::Float64 - function DynamicTol(alg::A, tol_min::Real, tol_max::Real, tol_factor::Real) where {A} + "factor on the local truncation error, setting a tolerance floor (0 disables it)" + truncation_factor::Float64 + + function DynamicTol( + 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")) + 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 """ - adapt_solver(alg::DynamicTol; iter, g_global, ...) + adapt_solver(alg::DynamicTol; iter, g_global, eps_trunc, ...) 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: +global gradient / convergence-error scalar `g_global`, damped by the iteration count, but never +below the truncation-error floor `truncation_factor · eps_trunc` (relevant for truncating two-site +sweeps; with the default `truncation_factor = 0` the floor is inactive): - tol = clamp(tol_factor · g_global / √iter, tol_min, tol_max) + tol = clamp(max(truncation_factor · eps_trunc, tol_factor · g_global / √iter), tol_min, 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) +function adapt_solver( + alg::DynamicTol; iter::Integer = 1, g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs... + ) + tol = clamp( + max(alg.truncation_factor * eps_trunc, alg.tol_factor * g_global / sqrt(max(iter, 1))), + alg.tol_min, alg.tol_max + ) return _updatetol(alg.alg, tol) end From b4f0fe8d2e3e742cdc71a291564a8d22d96ec927 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 26 Jul 2026 13:06:12 -0400 Subject: [PATCH 4/7] feat: MPS linear solver (linsolve) category Add a first-class `linsolve`/`linsolve!` for finite MPS, solving `(a0 + a1*A)*x = b` with `A` an MPO/MPOHamiltonian and `x`, `b` MPS. The signature is state-first (`linsolve(x0, A, b, [alg], [envs]; a0=0, a1=1, ...)`) and returns `(x, envs, eps)`, matching find_groundstate/approximate; the shift follows KrylovKit's convention. Algorithms are DMRG-style sweeps `DMRGSolve` (single-site) and `DMRGSolve2` (two-site, rank-adaptive), each parameterized by a local `LinsolveFormulation`: `Galerkin` (solve the effective system directly; CG/GMRES/BiCGStab) and `LeastSquares` (normal equations via squared environments). Convergence uses the local Galerkin residual `||(a0 + a1*A_eff)*AC - b_eff|| / ||b||`, formed as a local tensor (machine-precision floor, no catastrophic cancellation). The two-site stop is truncation-aware. Local solves default to adaptive tolerances via `Defaults.alg_linsolve()` (a structure-selected CG/MINRES/GMRES wrapped in `DynamicTol`, with a truncation floor for two-site); pass a plain solver to disable adaptation. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/MPSKit.jl | 5 + src/algorithms/linsolve/flinsolve.jl | 197 +++++++++++++++++++++++++++ src/algorithms/linsolve/linsolve.jl | 188 +++++++++++++++++++++++++ src/utility/defaults.jl | 18 ++- 4 files changed, 407 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/linsolve/flinsolve.jl create mode 100644 src/algorithms/linsolve/linsolve.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 920a8fefa..899a2bbb7 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -42,6 +42,8 @@ export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S export Zipup export propagator export DynamicalDMRG, NaiveInvert, Jeckelmann +export linsolve, linsolve! +export DMRGSolve, DMRGSolve2, Galerkin, LeastSquares export exact_diagonalization, fidelity_susceptibility # toolbox: @@ -203,6 +205,9 @@ include("algorithms/ED.jl") include("algorithms/unionalg.jl") +include("algorithms/linsolve/linsolve.jl") +include("algorithms/linsolve/flinsolve.jl") + include("utility/show.jl") function __init__() diff --git a/src/algorithms/linsolve/flinsolve.jl b/src/algorithms/linsolve/flinsolve.jl new file mode 100644 index 000000000..c8d1d6b5b --- /dev/null +++ b/src/algorithms/linsolve/flinsolve.jl @@ -0,0 +1,197 @@ +# Finite-MPS implementations of the linear solver `linsolve`. +# +# The sweep drivers mirror the variational structure of `find_groundstate!`/`corvector.jl`: at each +# site (single-site) or bond (two-site) we build the local effective operator from the operator +# environments and the local right-hand side from the overlap environments, solve a small local +# linear problem with `KrylovKit.linsolve`, and install the result. +# +# Convergence is measured, as in `find_groundstate!` (via `calc_galerkin`), by the local residual of +# the *original* system `(a₀ + a₁·A)·x = b`: at each site the norm of `(a₀ + a₁·A_eff)·AC − b_eff` +# (relative to `‖b‖`), evaluated on the current tensor before the local solve, maximized over the +# sweep. It is formed as a local tensor, so it resolves to machine precision even for +# shifted/resolvent systems (no `√eps` cancellation) and needs no MPS-level `+`/`-`. +# +# The local solves use adaptive *tolerances* by default: `adapt_solver` retunes the inner solver's +# tolerance per bond from the previous-sweep residual (`g_global`) — a `DynamicTol` wrapper (which +# works for any KrylovKit solver, `GMRES`/`CG`/`BiCGStab`/…, since it only sets `tol`). The Krylov +# budget is left fixed. For a plain (unwrapped) solver `adapt_solver` is the identity. + +# Right-hand-side context +# ----------------------- +# Bundles the environments a formulation needs beyond the operator sandwich `openvs`. `Galerkin` +# only needs the overlap `⟨x|b⟩`; `LeastSquares` additionally needs the squared operator (for +# `A²`) and the mixed sandwich `⟨x|A|b⟩` (for `A·b`). +_rhs_context(::Galerkin, x, A, b, openvs) = (; openvs, rhsenvs = environments(x, b)) +function _rhs_context(::LeastSquares, x, A, b, openvs) + Asq, sqenvs = squaredenvs(x, A, openvs) + # the JordanMPO effective operator requires `below === above`, so the `A·b` term of the + # normal-equation right-hand side is precomputed as an MPS and projected via the overlap form + Ab = A * b + return (; + openvs, rhsenvs = environments(x, b), Asq, sqenvs, Ab, abenvs = environments(x, Ab), + ) +end + +function _warn_unconverged(info, pos) + return info.converged == 0 && + @warn "linsolve: local solve at $pos did not converge (normres = $(info.normres))" +end + +# denominator for the relative residual; guards against a zero right-hand side +_rhs_norm(b) = (n = norm(b); iszero(n) ? one(n) : n) + +# Local single-site solves. Each returns the updated center tensor and the local residual of the +# *original* system (relative-residual convergence uses `res/‖b‖`; adaptation uses the previous +# sweep's residual via `g_global`). +function _local_linsolve(::Galerkin, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0) + A_eff = AC_hamiltonian(pos, x, A, x, ctx.openvs) + b_eff = AC_projection(pos, x, b, ctx.rhsenvs) + AC = x.AC[pos] + res = norm(a₀ * AC + a₁ * (A_eff * AC) - b_eff) + AC′, info = KrylovKit.linsolve(A_eff, b_eff, AC, adapt_solver(solver; iter, g_global), a₀, a₁) + _warn_unconverged(info, "site $pos") + return AC′, res +end +function _local_linsolve(::LeastSquares, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0) + A_eff = AC_hamiltonian(pos, x, A, x, ctx.openvs) + Asq_eff = AC_hamiltonian(pos, x, ctx.Asq, x, ctx.sqenvs) + N_eff = LinearCombination((A_eff, Asq_eff), (2 * real(conj(a₀) * a₁), abs2(a₁))) + b_eff = AC_projection(pos, x, b, ctx.rhsenvs) + Ab_eff = AC_projection(pos, x, ctx.Ab, ctx.abenvs) + AC = x.AC[pos] + # convergence uses the ORIGINAL-system residual, not the normal-equation residual + res = norm(a₀ * AC + a₁ * (A_eff * AC) - b_eff) + rhs = conj(a₀) * b_eff + conj(a₁) * Ab_eff + AC′, info = KrylovKit.linsolve(N_eff, rhs, AC, adapt_solver(solver; iter, g_global), abs2(a₀), one(a₁)) + _warn_unconverged(info, "site $pos") + return AC′, res +end + +# Local two-site solves. Return the updated two-site tensor and the original-system residual. +# `eps_trunc` (the previous discarded weight at this bond) floors the adaptive tolerance so the +# local solve is not driven far below what the SVD truncation will discard. +function _local_linsolve(::Galerkin, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0) + A_eff = AC2_hamiltonian(pos, x, A, x, ctx.openvs) + b_eff = AC2_projection(pos, x, b, ctx.rhsenvs) + ac2 = AC2(x, pos; kind) + res = norm(a₀ * ac2 + a₁ * (A_eff * ac2) - b_eff) + AC2′, info = KrylovKit.linsolve(A_eff, b_eff, ac2, adapt_solver(solver; iter, g_global, eps_trunc), a₀, a₁) + _warn_unconverged(info, "bond $pos") + return AC2′, res +end +function _local_linsolve(::LeastSquares, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0) + A_eff = AC2_hamiltonian(pos, x, A, x, ctx.openvs) + Asq_eff = AC2_hamiltonian(pos, x, ctx.Asq, x, ctx.sqenvs) + N_eff = LinearCombination((A_eff, Asq_eff), (2 * real(conj(a₀) * a₁), abs2(a₁))) + b_eff = AC2_projection(pos, x, b, ctx.rhsenvs) + Ab_eff = AC2_projection(pos, x, ctx.Ab, ctx.abenvs) + ac2 = AC2(x, pos; kind) + res = norm(a₀ * ac2 + a₁ * (A_eff * ac2) - b_eff) + rhs = conj(a₀) * b_eff + conj(a₁) * Ab_eff + AC2′, info = KrylovKit.linsolve(N_eff, rhs, ac2, adapt_solver(solver; iter, g_global, eps_trunc), abs2(a₀), one(a₁)) + _warn_unconverged(info, "bond $pos") + return AC2′, res +end + +# Single-site driver (bond-preserving, so no truncation term in the stop test) +# --------------------------------------------------------------------------- +function linsolve!( + x::AbstractFiniteMPS, A, b, alg::DMRGSolve, + envs = environments(x, A, x); a₀ = 0, a₁ = 1 + ) + ctx = _rhs_context(alg.formulation, x, A, b, envs) + normb = _rhs_norm(b) + ϵ::Float64 = 2 * alg.tol # relative residual, for the stop test + ϵ_global = Inf # previous-sweep absolute residual, drives the adaptive tolerance + log = IterLog("linsolve") + + LoggingExtras.withlevel(; alg.verbosity) do + @infov 2 loginit!(log, ϵ) + for iter in 1:(alg.maxiter) + ϵ = 0.0 + res_max = 0.0 + for pos in [1:(length(x) - 1); length(x):-1:2] + AC′, res = _local_linsolve( + alg.formulation, Val(1), pos, x, A, b, ctx, alg.solver, a₀, a₁; + iter, g_global = ϵ_global + ) + ϵ = max(ϵ, res / normb) + res_max = max(res_max, res) + x.AC[pos] = AC′ + end + + x, envs = alg.finalize(iter, x, A, envs)::Tuple{typeof(x), typeof(envs)} + ϵ_global = res_max + + if ϵ <= alg.tol + @infov 2 logfinish!(log, iter, ϵ) + break + end + if iter == alg.maxiter + @warnv 1 logcancel!(log, iter, ϵ) + else + @infov 3 logiter!(log, iter, ϵ) + end + end + end + + return x, envs, ϵ +end + +# Two-site driver (truncation-aware stop: the residual cannot beat the discarded weight) +# ------------------------------------------------------------------------------------- +function linsolve!( + x::AbstractFiniteMPS, A, b, alg::DMRGSolve2, + envs = environments(x, A, x); a₀ = 0, a₁ = 1 + ) + ctx = _rhs_context(alg.formulation, x, A, b, envs) + normb = _rhs_norm(b) + ϵ_truncs = zeros(length(x) - 1) # per-bond discarded weight + ϵ::Float64 = 2 * alg.tol + ϵ_global = Inf + log = IterLog("linsolve2") + + LoggingExtras.withlevel(; alg.verbosity) do + @infov 2 loginit!(log, ϵ) + for iter in 1:(alg.maxiter) + ϵ = 0.0 + res_max = 0.0 + for pos in 1:(length(x) - 1) + AC2′, res = _local_linsolve( + alg.formulation, Val(2), pos, x, A, b, ctx, alg.solver, a₀, a₁, :ACAR; + iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos] + ) + ϵ = max(ϵ, res / normb) + res_max = max(res_max, res) + x, tr = gauge2!(x, pos, Val(:right), AC2′, alg.alg_gauge; normalize = false) + ϵ_truncs[pos] = tr + end + for pos in (length(x) - 2):-1:1 + AC2′, res = _local_linsolve( + alg.formulation, Val(2), pos, x, A, b, ctx, alg.solver, a₀, a₁, :ALAC; + iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos] + ) + ϵ = max(ϵ, res / normb) + res_max = max(res_max, res) + x, tr = gauge2!(x, pos, Val(:left), AC2′, alg.alg_gauge; normalize = false) + ϵ_truncs[pos] = tr + end + + x, envs = alg.finalize(iter, x, A, envs)::Tuple{typeof(x), typeof(envs)} + ϵ_global = res_max + + # the Galerkin residual cannot drop below the level set by the discarded weight + if ϵ <= max(alg.tol, maximum(ϵ_truncs) / normb) + @infov 2 logfinish!(log, iter, ϵ) + break + end + if iter == alg.maxiter + @warnv 1 logcancel!(log, iter, ϵ) + else + @infov 3 logiter!(log, iter, ϵ) + end + end + end + + return x, envs, ϵ +end diff --git a/src/algorithms/linsolve/linsolve.jl b/src/algorithms/linsolve/linsolve.jl new file mode 100644 index 000000000..cdbf9ecc6 --- /dev/null +++ b/src/algorithms/linsolve/linsolve.jl @@ -0,0 +1,188 @@ +# Flavours +# -------- +""" +$(TYPEDEF) + +Abstract supertype selecting how the local linear subproblem is posed in an MPS [`linsolve`](@ref) +sweep. See [`Galerkin`](@ref) and [`LeastSquares`](@ref). +""" +abstract type LinsolveFormulation end + +""" +$(TYPEDEF) + +Local formulation imposing the Galerkin condition: the residual `(a₀ + a₁·A)·x − b` is projected +orthogonal to the local tangent space, yielding the effective system `(a₀ + a₁·A_eff)·x = b_eff` +which is solved with the configured Krylov `solver`. Use a `CG` solver for hermitian +positive-definite `A`, and `GMRES`/`BiCGStab` for general (non-hermitian or indefinite) `A`. +""" +struct Galerkin <: LinsolveFormulation end + +""" +$(TYPEDEF) + +Local formulation that minimizes `‖(a₀ + a₁·A)·x − b‖²`, i.e. solves the normal equations +`M†M·x = M†b` with `M = a₀ + a₁·A`, built from squared operator environments. The normal operator +is unconditionally positive-definite (solve it with `CG`) at the cost of squaring the condition +number. + +Currently assumes `A` is hermitian (the resolvent / dynamical-DMRG case), for which +`M†M = |a₀|² + 2·Re(ā₀·a₁)·A + |a₁|²·A²` and `M†b = ā₀·b + ā₁·A·b`. This is a generalization of the +[`Jeckelmann`](@ref) dynamical-DMRG functional. +""" +struct LeastSquares <: LinsolveFormulation end + +# Algorithms +# ---------- +""" +$(TYPEDEF) + +Single-site DMRG-style sweeping algorithm for the MPS linear solver [`linsolve`](@ref). It keeps +the bond dimension of the initial guess fixed and, at each site, solves the local linear problem +selected by `formulation` with the local Krylov `solver`. + +# Fields + +$(TYPEDFIELDS) +""" +@kwdef struct DMRGSolve{F <: LinsolveFormulation, S, FIN} <: Algorithm + "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)" + formulation::F = Galerkin() + "local linear solver; a plain KrylovKit solver, or one wrapped in [`DynamicTol`](@ref) for per-bond adaptive tolerances (the default)" + solver::S = Defaults.alg_linsolve() + "tolerance for convergence criterium" + tol::Float64 = Defaults.tol + "maximal amount of iterations" + maxiter::Int = Defaults.maxiter + "setting for how much information is displayed" + verbosity::Int = Defaults.verbosity + "callback function applied after each iteration, of signature `finalize(iter, x, A, envs) -> x, envs`" + finalize::FIN = Defaults._finalize +end + +""" +$(TYPEDEF) + +Two-site DMRG-style sweeping algorithm for the MPS linear solver [`linsolve`](@ref). Each bond +update solves the local two-site linear problem selected by `formulation` and truncates the enlarged +bond back down with `alg_gauge` (a truncated SVD built from `trunc`), making the bond dimension +adaptive. + +# Fields + +$(TYPEDFIELDS) +""" +struct DMRGSolve2{F <: LinsolveFormulation, S, G, FIN} <: Algorithm + "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)" + formulation::F + "local linear solver; a plain KrylovKit solver, or one wrapped in [`DynamicTol`](@ref) for per-bond adaptive tolerances" + solver::S + "tolerance for convergence criterium" + tol::Float64 + "maximal amount of iterations" + maxiter::Int + "setting for how much information is displayed" + verbosity::Int + "truncated SVD used for the post-update gauge" + alg_gauge::G + "callback function applied after each iteration, of signature `finalize(iter, x, A, envs) -> x, envs`" + finalize::FIN +end +function DMRGSolve2(; + formulation = Galerkin(), solver = Defaults.alg_linsolve(), tol = Defaults.tol, + maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, + alg_svd = Defaults.alg_svd(), trunc, finalize = Defaults._finalize + ) + alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc) + return DMRGSolve2(formulation, solver, tol, maxiter, verbosity, alg_gauge, finalize) +end + +# Interface +# --------- +@doc """ + linsolve(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, ϵ) + linsolve!(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, ϵ) + +Solve the linear system `(a₀ + a₁·A)·x = b` for the MPS `x`, where `A` is an MPO / MPOHamiltonian +and `b` is an MPS. The `a₀ + a₁·A` shift follows `KrylovKit.linsolve` and is applied implicitly; +the defaults `a₀ = 0`, `a₁ = 1` give the plain system `A·x = b`. `linsolve!` solves in-place, +overwriting `x₀`. + +The initial guess `x₀` comes first, mirroring [`find_groundstate`](@ref) and [`approximate`](@ref) +(the operator/right-hand side follow), so it sits in the same argument slot in both the copying and +in-place forms. + +Currently only finite MPS are supported. + +# Arguments +- `x₀::AbstractMPS`: initial guess (and, for `linsolve!`, the state overwritten with the solution) +- `A`: operator of the linear system (MPO / MPOHamiltonian) +- `b::AbstractMPS`: right-hand side +- `algorithm`: linear-solve algorithm, see [`DMRGSolve`](@ref) and [`DMRGSolve2`](@ref) +- `[environments]`: MPS environment manager for the operator sandwich `⟨x|A|x⟩` + +# Keyword Arguments +The keyword-based call (no explicit `algorithm`) selects an algorithm from the structure flags: +- `tol::Float64`: tolerance for convergence criterium +- `maxiter::Int`: maximum amount of iterations +- `verbosity::Int`: display progress information +- `a₀`, `a₁`: shift/scale scalars of the system `(a₀ + a₁·A)·x = b` +- `ishermitian::Bool`, `isposdef::Bool`: declare structure of `(a₀ + a₁·A)` so the default local + solver is chosen accordingly (`CG` for positive-definite, `MINRES` for hermitian-indefinite, + `GMRES` otherwise) +- `trunc`: if supplied, a truncated two-site sweep ([`DMRGSolve2`](@ref)) is prepended to + adapt the bond dimension before the single-site algorithm polishes the result + +# Returns +- `x::AbstractMPS`: the (bond-dimension-limited) solution +- `environments`: operator environments corresponding to `x` +- `ϵ::Float64`: final convergence error — the largest local residual `‖(a₀ + a₁·A)·x − b‖` + over the sweep, relative to `‖b‖` (the linear-solve analogue of the Galerkin error used by + [`find_groundstate`](@ref)) +""" linsolve, linsolve! + +# scalar-promote an initial state so it can hold a complex solution when the shift is complex +function _promote_state(x::AbstractMPS, a₀, a₁) + T = promote_type(scalartype(x), typeof(a₀), typeof(a₁)) + return T <: Complex ? complex(x) : x +end + +# default (adaptive) local solver from the declared operator structure, mirroring KrylovKit's rule +function _default_linsolve_algorithm( + x::AbstractMPS; tol, maxiter, verbosity, ishermitian, isposdef, trunc + ) + x isa AbstractFiniteMPS || + throw(ArgumentError("`linsolve` currently only supports finite MPS")) + solver = Defaults.alg_linsolve(; ishermitian, isposdef, tol, maxiter) + alg = DMRGSolve(; solver, tol, maxiter, verbosity) + if !isnothing(trunc) + alg = DMRGSolve2(; + solver, tol = min(1.0e-2, 100tol), maxiter, verbosity, trunc + ) & alg + end + return alg +end + +# keyword form: build a default algorithm from the initial guess and dispatch +function linsolve( + x₀::AbstractMPS, A, b::AbstractMPS; + tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, + a₀ = 0, a₁ = 1, ishermitian = false, isposdef = false, trunc = nothing + ) + alg = _default_linsolve_algorithm( + x₀; tol, maxiter, verbosity, ishermitian, isposdef, trunc + ) + return linsolve(x₀, A, b, alg; a₀, a₁) +end + +# explicit-algorithm form: copy (and scalar-promote) the guess and solve in-place. The `envs` +# splat lets `linsolve!` build the environments from the promoted copy when none are supplied. +function linsolve(x₀, A, b, alg::Union{DMRGSolve, DMRGSolve2}, envs...; a₀ = 0, a₁ = 1) + return linsolve!(_promote_state(copy(x₀), a₀, a₁), A, b, alg, envs...; a₀, a₁) +end + +# sequential chaining of algorithms (e.g. a two-site pass then a single-site polish) +function linsolve(x₀, A, b, alg::UnionAlg, envs...; a₀ = 0, a₁ = 1) + x, newenvs, = linsolve(x₀, A, b, alg.alg1, envs...; a₀, a₁) + return linsolve(x, A, b, alg.alg2, newenvs; a₀, a₁) +end diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index dbc221420..81984de45 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -5,7 +5,7 @@ Some default values and settings for MPSKit. """ module Defaults -import KrylovKit: GMRES, Arnoldi, Lanczos +import KrylovKit: GMRES, CG, MINRES, Arnoldi, Lanczos using OhMyThreads using ..MPSKit: DynamicTol, AdaptiveKrylov using TensorKit: TensorKit @@ -62,6 +62,22 @@ function alg_eigsolve(; return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end +function alg_linsolve(; + ishermitian = false, isposdef = false, tol = tol, maxiter = maxiter, + krylovdim = krylovdim, verbosity = 0, adaptive = true, dynamic_tols = dynamic_tols, + tol_min = tol_min, tol_max = tol_max, tol_factor = eigs_tolfactor, truncation_factor = 1.0e-1 + ) + # pick the local linear solver from the declared operator structure (mirrors KrylovKit's rule) + base = isposdef ? CG(; tol, maxiter, verbosity) : + ishermitian ? MINRES(; tol, maxiter, verbosity) : + GMRES(; tol, maxiter, krylovdim, verbosity) + # linsolve defaults to per-bond adaptive *tolerances* (only the inner tol is retuned, from the + # previous-sweep residual and, for truncating two-site sweeps, floored at the truncation error); + # the Krylov budget is left fixed. Opt out with `adaptive = false`. + return (adaptive && dynamic_tols) ? + DynamicTol(base, tol_min, tol_max, tol_factor, truncation_factor) : base +end + function alg_environments(; tol = tol, maxiter = maxiter, verbosity = 0, krylovdim = krylovdim, eager = true, dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max, From 001ee79035b78561e783c1ef432b2e650a6e8a43 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 26 Jul 2026 13:06:26 -0400 Subject: [PATCH 5/7] test: linsolve on finite MPS Cover the resolvent `(z - H)x = V|gs>` against the analytic `1/(z - E0)` and the existing `propagator` (both Galerkin and LeastSquares, single- and two-site), a positive-definite solve via the keyword interface (CG path), adaptive tolerances including BiCGStab, and a general RHS checked against a dense solve. Residuals are verified densely (MPS-level +/- does not reliably cancel, cf. #473). Co-Authored-By: Claude Opus 4.8 (1M context) --- test/algorithms/linsolve.jl | 102 ++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 test/algorithms/linsolve.jl diff --git a/test/algorithms/linsolve.jl b/test/algorithms/linsolve.jl new file mode 100644 index 000000000..5b9d8b920 --- /dev/null +++ b/test/algorithms/linsolve.jl @@ -0,0 +1,102 @@ +println(" +----------------------------- +| Linear solver tests | +----------------------------- +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using TensorKit +using TensorKit: ℙ +# KrylovKit is not a direct test dependency; reach the solver types through MPSKit's re-export +const GMRES = MPSKit.KrylovKit.GMRES +const CG = MPSKit.KrylovKit.CG +const BiCGStab = MPSKit.KrylovKit.BiCGStab +const DynamicTol = MPSKit.DynamicTol + +verbosity_conv = 1 + +# helper: relative residual of (a₀ + a₁·A)·x = b, computed densely (robust for small systems; +# MPS-level `+`/`-` does not reliably cancel across states of differing bond structure) +function rel_residual(A, x, b, a₀, a₁) + Am = convert(TensorMap, A) + xv = convert(TensorMap, x) + bv = convert(TensorMap, b) + return norm(a₀ * xv + a₁ * (Am * xv) - bv) / norm(bv) +end + +@testset "linsolve FiniteMPS" verbose = true begin + L = 8 + H = force_planar(-transverse_field_ising(; L, g = -4)) + gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^16), H; verbosity = verbosity_conv) + E₀ = real(expectation_value(gs, H)) + + # resolvent of an eigenstate: (z − H)⁻¹|gs⟩ = 1/(z − E₀)|gs⟩, so ⟨gs|x⟩ = 1/(z − E₀). + z = E₀ + 0.5 + 0.3im + predicted = 1 / (z - E₀) + + @testset "resolvent, single-site, formulation $flav" for flav in (Galerkin(), LeastSquares()) + solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12) + alg = DMRGSolve(; formulation = flav, solver, tol = 1.0e-10, verbosity = 0) + # solve (z − H) x = gs ⟺ a₀ = z, a₁ = −1 + x, envs, ϵ = linsolve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1) + @test dot(gs, x) ≈ predicted atol = 1.0e-6 + @test rel_residual(H, x, gs, z, -1) < 1.0e-6 + end + + @testset "resolvent, two-site (adaptive χ)" begin + alg = DMRGSolve2(; + formulation = Galerkin(), solver = GMRES(; tol = 1.0e-12), + trunc = truncrank(16), tol = 1.0e-10, verbosity = 0 + ) + x, = linsolve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1) + @test dot(gs, x) ≈ predicted atol = 1.0e-5 + @test rel_residual(H, x, gs, z, -1) < 1.0e-5 + end + + @testset "matches propagator (correction vector)" begin + alg_ls = DMRGSolve(; solver = GMRES(; tol = 1.0e-12), tol = 1.0e-10, verbosity = 0) + x, = linsolve(complex(copy(gs)), H, gs, alg_ls; a₀ = z, a₁ = -1) + for f in (NaiveInvert(), Jeckelmann()) + g_prop, = propagator(gs, z, H, DynamicalDMRG(; flavour = f, tol = 1.0e-10, verbosity = 0)) + @test dot(gs, x) ≈ g_prop atol = 1.0e-5 + end + end + + @testset "adaptive tolerances (incl. BiCGStab)" begin + # DynamicTol works for any solver, including short-recurrence BiCGStab (retunes only `tol`) + alg_bicg = DMRGSolve(; + solver = DynamicTol(BiCGStab(; tol = 1.0e-12); tol_min = 1.0e-14, tol_max = 1.0e-4), + tol = 1.0e-9, verbosity = 0 + ) + x, = linsolve(complex(copy(gs)), H, gs, alg_bicg; a₀ = z, a₁ = -1) + @test dot(gs, x) ≈ predicted atol = 1.0e-5 + @test rel_residual(H, x, gs, z, -1) < 1.0e-5 + + # adaptive-by-default keyword path (DynamicTol-wrapped GMRES) + x2, envs, ϵ = linsolve(complex(copy(gs)), H, gs; a₀ = z, a₁ = -1, tol = 1.0e-9, verbosity = 0) + @test dot(gs, x2) ≈ predicted atol = 1.0e-5 + @test rel_residual(H, x2, gs, z, -1) < 1.0e-5 + end + + @testset "positive-definite solve via keyword interface" begin + # s ≫ spectrum ⇒ (s − H) is positive-definite; the CG path is selected by `isposdef` + s = E₀ + 50.0 + # eigenstate RHS: exact solution 1/(s − E₀) · gs, real and positive + x, envs, ϵ = linsolve(copy(gs), H, gs; a₀ = s, a₁ = -1, isposdef = true, tol = 1.0e-10, verbosity = 0) + @test dot(gs, x) ≈ 1 / (s - E₀) atol = 1.0e-8 + @test rel_residual(H, x, gs, s, -1) < 1.0e-8 + end + + @testset "general RHS, plain A·x = b against dense" begin + # shift H to a well-conditioned, invertible operator A = s·I − H (s ≫ spectrum) + s = E₀ + 50.0 + b = normalize!(FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^8)) + alg = DMRGSolve2(; + solver = CG(; tol = 1.0e-12), trunc = truncrank(64), tol = 1.0e-9, verbosity = 0 + ) & DMRGSolve(; solver = CG(; tol = 1.0e-12), tol = 1.0e-9, verbosity = 0) + x, = linsolve(complex(copy(b)), H, b, alg; a₀ = s, a₁ = -1) + @test rel_residual(H, x, b, s, -1) < 1.0e-4 + end +end From ea18881621696b4ca59c977b2a456f2efad4b4fc Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 27 Jul 2026 19:03:27 -0400 Subject: [PATCH 6/7] =?UTF-8?q?fix:=20linsolve=20local=20solves,=20and=20f?= =?UTF-8?q?orm=20A=E2=80=A0b=20without=20materializing=20A*b?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three fixes to the local solves: - `Defaults.alg_linsolve(; ishermitian = true)` returned a `MINRES` solver, but KrylovKit exports the type without implementing `linsolve` for it (see its own "TODO: implement MINRES for symmetric but not posdef; for now use GMRES"), so the hermitian-indefinite path threw a `MethodError`. Fall back to GMRES too. - `@doc str a, b` attaches only to the last binding, so exported `linsolve` had no docstring at all. Share it explicitly with `@doc (@doc linsolve) linsolve!`. - A local solve no longer warns about non-convergence when it undershot its own adaptive tolerance but still reached `tol * norm(b)`, the accuracy the outer sweep needs. Adaptive tolerances routinely dip below the round-off floor of the local problem, which made the warning fire on nearly every tightly converged sweep. `LeastSquares` now forms the `A†b` term of the normal equations from the mixed sandwich instead of materializing `A*b` as an MPS, avoiding a state of bond dimension chi_A*chi_b together with its environments. The JordanMPO effective operator assumes `below === above` and branching on that inside `AC_hamiltonian` would make the return type of that hot function a union, so the mixed case goes through the generic sparse MPO derivative in `_mixed_projection` -- one application per bond, so its cost is irrelevant. This also makes the formulation usable on `WindowMPS`, which the `A*b` route never could be: there is no `*(::WindowMPOHamiltonian, ::WindowMPS)`. Also adds the missing in-place `linsolve!(x0, A, b, ::UnionAlg, envs...)`, the counterpart of the existing copying method. Tests cover two-site `LeastSquares` and a `WindowMPS` resolvent for both formulations. Co-Authored-By: Claude Opus 5 (1M context) --- .../derivatives/hamiltonian_derivatives.jl | 6 +- src/algorithms/linsolve/flinsolve.jl | 66 +++++++++++++------ src/algorithms/linsolve/linsolve.jl | 16 +++-- src/utility/defaults.jl | 8 ++- test/algorithms/linsolve.jl | 28 +++++++- 5 files changed, 92 insertions(+), 32 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index f965db9ac..f63e49ce9 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -39,6 +39,10 @@ function AC_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs; prepare::Bool = true ) + # NOTE: the JordanMPO fast path genuinely assumes a single state. A mixed sandwich + # ⟨below|H|above⟩ has to go through the generic sparse contraction, but branching on it here + # would make the return type of this (hot) function a union; see `_mixed_projection` in + # `algorithms/linsolve/flinsolve.jl` for the mixed-sandwich path. @assert below === above "JordanMPO assumptions break" GL = leftenv(envs, site, below) GR = rightenv(envs, site, below) @@ -191,7 +195,7 @@ function AC2_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs; prepare::Bool = true ) - @assert below === above "JordanMPO assumptions break" + @assert below === above "JordanMPO assumptions break" # see `AC_hamiltonian` GL = leftenv(envs, site, below) GR = rightenv(envs, site + 1, below) W1, W2 = operator[site], operator[site + 1] diff --git a/src/algorithms/linsolve/flinsolve.jl b/src/algorithms/linsolve/flinsolve.jl index c8d1d6b5b..895490ca5 100644 --- a/src/algorithms/linsolve/flinsolve.jl +++ b/src/algorithms/linsolve/flinsolve.jl @@ -24,16 +24,40 @@ _rhs_context(::Galerkin, x, A, b, openvs) = (; openvs, rhsenvs = environments(x, b)) function _rhs_context(::LeastSquares, x, A, b, openvs) Asq, sqenvs = squaredenvs(x, A, openvs) - # the JordanMPO effective operator requires `below === above`, so the `A·b` term of the - # normal-equation right-hand side is precomputed as an MPS and projected via the overlap form - Ab = A * b - return (; - openvs, rhsenvs = environments(x, b), Asq, sqenvs, Ab, abenvs = environments(x, Ab), + # the `A·b` term of the normal-equation right-hand side comes from the mixed sandwich + # `⟨x|A|b⟩`, so `A·b` is never materialized as an MPS + return (; openvs, rhsenvs = environments(x, b), Asq, sqenvs, abenvs = environments(x, A, b)) +end + +# `⟨∂x|A|b⟩`, the mixed-sandwich projection supplying the `A†b` term of the normal equations. +# `AC_hamiltonian(pos, x, A, b, envs)` cannot be used for an `MPOHamiltonian`: the JordanMPO +# effective operator assumes `below === above`, and teaching it to branch would make the return type +# of that hot function a union. The generic sparse MPO derivative handles the mixed case, and one +# application per bond makes its cost irrelevant. Window operators store the finite part in the +# environments, so unwrap them the same way the derivative constructors do. +_finite_operator(A) = A +_finite_operator(A::WindowMPOHamiltonian) = A.finite_ham + +function _mixed_projection(::Val{1}, pos, x, A, b, envs) + W = _finite_operator(A) + H = MPO_AC_Hamiltonian(leftenv(envs, pos, x), W[pos], rightenv(envs, pos, x)) + return H * b.AC[pos] +end +function _mixed_projection(::Val{2}, pos, x, A, b, envs) + W = _finite_operator(A) + H = MPO_AC2_Hamiltonian( + leftenv(envs, pos, x), W[pos], W[pos + 1], rightenv(envs, pos + 1, x) ) + return H * AC2(b, pos) end -function _warn_unconverged(info, pos) - return info.converged == 0 && +# KrylovKit's `linsolve` is silent on non-convergence, so check it here. `warn_tol` is the absolute +# residual the outer sweep actually needs (`alg.tol · ‖b‖`): an inner solve that undershoots its own +# adaptive tolerance but still beats that has done its job, and warning about it would fire on nearly +# every tightly-converged sweep (adaptive tolerances routinely dip below the round-off floor of the +# local problem). +function _warn_unconverged(info, pos, warn_tol) + return info.converged == 0 && info.normres > warn_tol && @warn "linsolve: local solve at $pos did not converge (normres = $(info.normres))" end @@ -43,53 +67,53 @@ _rhs_norm(b) = (n = norm(b); iszero(n) ? one(n) : n) # Local single-site solves. Each returns the updated center tensor and the local residual of the # *original* system (relative-residual convergence uses `res/‖b‖`; adaptation uses the previous # sweep's residual via `g_global`). -function _local_linsolve(::Galerkin, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0) +function _local_linsolve(::Galerkin, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0, warn_tol = 0.0) A_eff = AC_hamiltonian(pos, x, A, x, ctx.openvs) b_eff = AC_projection(pos, x, b, ctx.rhsenvs) AC = x.AC[pos] res = norm(a₀ * AC + a₁ * (A_eff * AC) - b_eff) AC′, info = KrylovKit.linsolve(A_eff, b_eff, AC, adapt_solver(solver; iter, g_global), a₀, a₁) - _warn_unconverged(info, "site $pos") + _warn_unconverged(info, "site $pos", warn_tol) return AC′, res end -function _local_linsolve(::LeastSquares, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0) +function _local_linsolve(::LeastSquares, ::Val{1}, pos, x, A, b, ctx, solver, a₀, a₁; iter = 1, g_global = 0.0, warn_tol = 0.0) A_eff = AC_hamiltonian(pos, x, A, x, ctx.openvs) Asq_eff = AC_hamiltonian(pos, x, ctx.Asq, x, ctx.sqenvs) N_eff = LinearCombination((A_eff, Asq_eff), (2 * real(conj(a₀) * a₁), abs2(a₁))) b_eff = AC_projection(pos, x, b, ctx.rhsenvs) - Ab_eff = AC_projection(pos, x, ctx.Ab, ctx.abenvs) + Ab_eff = _mixed_projection(Val(1), pos, x, A, b, ctx.abenvs) AC = x.AC[pos] # convergence uses the ORIGINAL-system residual, not the normal-equation residual res = norm(a₀ * AC + a₁ * (A_eff * AC) - b_eff) rhs = conj(a₀) * b_eff + conj(a₁) * Ab_eff AC′, info = KrylovKit.linsolve(N_eff, rhs, AC, adapt_solver(solver; iter, g_global), abs2(a₀), one(a₁)) - _warn_unconverged(info, "site $pos") + _warn_unconverged(info, "site $pos", warn_tol) return AC′, res end # Local two-site solves. Return the updated two-site tensor and the original-system residual. # `eps_trunc` (the previous discarded weight at this bond) floors the adaptive tolerance so the # local solve is not driven far below what the SVD truncation will discard. -function _local_linsolve(::Galerkin, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0) +function _local_linsolve(::Galerkin, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0, warn_tol = 0.0) A_eff = AC2_hamiltonian(pos, x, A, x, ctx.openvs) b_eff = AC2_projection(pos, x, b, ctx.rhsenvs) ac2 = AC2(x, pos; kind) res = norm(a₀ * ac2 + a₁ * (A_eff * ac2) - b_eff) AC2′, info = KrylovKit.linsolve(A_eff, b_eff, ac2, adapt_solver(solver; iter, g_global, eps_trunc), a₀, a₁) - _warn_unconverged(info, "bond $pos") + _warn_unconverged(info, "bond $pos", warn_tol) return AC2′, res end -function _local_linsolve(::LeastSquares, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0) +function _local_linsolve(::LeastSquares, ::Val{2}, pos, x, A, b, ctx, solver, a₀, a₁, kind; iter = 1, g_global = 0.0, eps_trunc = 0.0, warn_tol = 0.0) A_eff = AC2_hamiltonian(pos, x, A, x, ctx.openvs) Asq_eff = AC2_hamiltonian(pos, x, ctx.Asq, x, ctx.sqenvs) N_eff = LinearCombination((A_eff, Asq_eff), (2 * real(conj(a₀) * a₁), abs2(a₁))) b_eff = AC2_projection(pos, x, b, ctx.rhsenvs) - Ab_eff = AC2_projection(pos, x, ctx.Ab, ctx.abenvs) + Ab_eff = _mixed_projection(Val(2), pos, x, A, b, ctx.abenvs) ac2 = AC2(x, pos; kind) res = norm(a₀ * ac2 + a₁ * (A_eff * ac2) - b_eff) rhs = conj(a₀) * b_eff + conj(a₁) * Ab_eff AC2′, info = KrylovKit.linsolve(N_eff, rhs, ac2, adapt_solver(solver; iter, g_global, eps_trunc), abs2(a₀), one(a₁)) - _warn_unconverged(info, "bond $pos") + _warn_unconverged(info, "bond $pos", warn_tol) return AC2′, res end @@ -101,6 +125,7 @@ function linsolve!( ) ctx = _rhs_context(alg.formulation, x, A, b, envs) normb = _rhs_norm(b) + warn_tol = alg.tol * normb # accuracy the outer sweep actually needs ϵ::Float64 = 2 * alg.tol # relative residual, for the stop test ϵ_global = Inf # previous-sweep absolute residual, drives the adaptive tolerance log = IterLog("linsolve") @@ -113,7 +138,7 @@ function linsolve!( for pos in [1:(length(x) - 1); length(x):-1:2] AC′, res = _local_linsolve( alg.formulation, Val(1), pos, x, A, b, ctx, alg.solver, a₀, a₁; - iter, g_global = ϵ_global + iter, g_global = ϵ_global, warn_tol ) ϵ = max(ϵ, res / normb) res_max = max(res_max, res) @@ -146,6 +171,7 @@ function linsolve!( ) ctx = _rhs_context(alg.formulation, x, A, b, envs) normb = _rhs_norm(b) + warn_tol = alg.tol * normb # accuracy the outer sweep actually needs ϵ_truncs = zeros(length(x) - 1) # per-bond discarded weight ϵ::Float64 = 2 * alg.tol ϵ_global = Inf @@ -159,7 +185,7 @@ function linsolve!( for pos in 1:(length(x) - 1) AC2′, res = _local_linsolve( alg.formulation, Val(2), pos, x, A, b, ctx, alg.solver, a₀, a₁, :ACAR; - iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos] + iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos], warn_tol ) ϵ = max(ϵ, res / normb) res_max = max(res_max, res) @@ -169,7 +195,7 @@ function linsolve!( for pos in (length(x) - 2):-1:1 AC2′, res = _local_linsolve( alg.formulation, Val(2), pos, x, A, b, ctx, alg.solver, a₀, a₁, :ALAC; - iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos] + iter, g_global = ϵ_global, eps_trunc = ϵ_truncs[pos], warn_tol ) ϵ = max(ϵ, res / normb) res_max = max(res_max, res) diff --git a/src/algorithms/linsolve/linsolve.jl b/src/algorithms/linsolve/linsolve.jl index cdbf9ecc6..7e7b70461 100644 --- a/src/algorithms/linsolve/linsolve.jl +++ b/src/algorithms/linsolve/linsolve.jl @@ -48,7 +48,7 @@ $(TYPEDFIELDS) @kwdef struct DMRGSolve{F <: LinsolveFormulation, S, FIN} <: Algorithm "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)" formulation::F = Galerkin() - "local linear solver; a plain KrylovKit solver, or one wrapped in [`DynamicTol`](@ref) for per-bond adaptive tolerances (the default)" + "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances (the default)" solver::S = Defaults.alg_linsolve() "tolerance for convergence criterium" tol::Float64 = Defaults.tol @@ -75,7 +75,7 @@ $(TYPEDFIELDS) struct DMRGSolve2{F <: LinsolveFormulation, S, G, FIN} <: Algorithm "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)" formulation::F - "local linear solver; a plain KrylovKit solver, or one wrapped in [`DynamicTol`](@ref) for per-bond adaptive tolerances" + "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances" solver::S "tolerance for convergence criterium" tol::Float64 @@ -128,8 +128,7 @@ The keyword-based call (no explicit `algorithm`) selects an algorithm from the s - `verbosity::Int`: display progress information - `a₀`, `a₁`: shift/scale scalars of the system `(a₀ + a₁·A)·x = b` - `ishermitian::Bool`, `isposdef::Bool`: declare structure of `(a₀ + a₁·A)` so the default local - solver is chosen accordingly (`CG` for positive-definite, `MINRES` for hermitian-indefinite, - `GMRES` otherwise) + solver is chosen accordingly (`CG` for positive-definite, `GMRES` otherwise) - `trunc`: if supplied, a truncated two-site sweep ([`DMRGSolve2`](@ref)) is prepended to adapt the bond dimension before the single-site algorithm polishes the result @@ -139,7 +138,10 @@ The keyword-based call (no explicit `algorithm`) selects an algorithm from the s - `ϵ::Float64`: final convergence error — the largest local residual `‖(a₀ + a₁·A)·x − b‖` over the sweep, relative to `‖b‖` (the linear-solve analogue of the Galerkin error used by [`find_groundstate`](@ref)) -""" linsolve, linsolve! +""" linsolve + +# NOTE: `@doc str a, b` attaches only to the last binding, so share it explicitly +@doc (@doc linsolve) linsolve! # scalar-promote an initial state so it can hold a complex solution when the shift is complex function _promote_state(x::AbstractMPS, a₀, a₁) @@ -186,3 +188,7 @@ function linsolve(x₀, A, b, alg::UnionAlg, envs...; a₀ = 0, a₁ = 1) x, newenvs, = linsolve(x₀, A, b, alg.alg1, envs...; a₀, a₁) return linsolve(x, A, b, alg.alg2, newenvs; a₀, a₁) end +function linsolve!(x₀, A, b, alg::UnionAlg, envs...; a₀ = 0, a₁ = 1) + x, newenvs, = linsolve!(x₀, A, b, alg.alg1, envs...; a₀, a₁) + return linsolve!(x, A, b, alg.alg2, newenvs; a₀, a₁) +end diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index 81984de45..8119f9105 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -5,7 +5,7 @@ Some default values and settings for MPSKit. """ module Defaults -import KrylovKit: GMRES, CG, MINRES, Arnoldi, Lanczos +import KrylovKit: GMRES, CG, Arnoldi, Lanczos using OhMyThreads using ..MPSKit: DynamicTol, AdaptiveKrylov using TensorKit: TensorKit @@ -67,9 +67,11 @@ function alg_linsolve(; krylovdim = krylovdim, verbosity = 0, adaptive = true, dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max, tol_factor = eigs_tolfactor, truncation_factor = 1.0e-1 ) - # pick the local linear solver from the declared operator structure (mirrors KrylovKit's rule) + # pick the local linear solver from the declared operator structure (mirrors KrylovKit's rule). + # NOTE: KrylovKit exports a `MINRES` type but has no `linsolve` method for it yet (see its own + # "TODO: implement MINRES for symmetric but not posdef; for now use GMRES"), so the + # hermitian-indefinite case falls back to GMRES as well. base = isposdef ? CG(; tol, maxiter, verbosity) : - ishermitian ? MINRES(; tol, maxiter, verbosity) : GMRES(; tol, maxiter, krylovdim, verbosity) # linsolve defaults to per-bond adaptive *tolerances* (only the inner tol is retuned, from the # previous-sweep residual and, for truncating two-site sweeps, floored at the truncation error); diff --git a/test/algorithms/linsolve.jl b/test/algorithms/linsolve.jl index 5b9d8b920..e82fa2d9d 100644 --- a/test/algorithms/linsolve.jl +++ b/test/algorithms/linsolve.jl @@ -45,10 +45,11 @@ end @test rel_residual(H, x, gs, z, -1) < 1.0e-6 end - @testset "resolvent, two-site (adaptive χ)" begin + @testset "resolvent, two-site (adaptive χ), formulation $flav" for flav in + (Galerkin(), LeastSquares()) + solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12) alg = DMRGSolve2(; - formulation = Galerkin(), solver = GMRES(; tol = 1.0e-12), - trunc = truncrank(16), tol = 1.0e-10, verbosity = 0 + formulation = flav, solver, trunc = truncrank(16), tol = 1.0e-10, verbosity = 0 ) x, = linsolve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1) @test dot(gs, x) ≈ predicted atol = 1.0e-5 @@ -100,3 +101,24 @@ end @test rel_residual(H, x, b, s, -1) < 1.0e-4 end end + +@testset "linsolve WindowMPS" verbose = true begin + # `LeastSquares` builds the `A†b` term from the mixed sandwich ⟨x|A|b⟩, which is the only way to + # reach it for a window operator (there is no `*(::WindowMPOHamiltonian, ::WindowMPS)`) + N = 10 + H = transverse_field_ising(; g = -4) + Ω, = find_groundstate(InfiniteMPS(ℂ^2, ℂ^10), H, VUMPS(; verbosity = 0)) + XΩ = WindowMPS(Ω, N) + H_w = WindowMPOHamiltonian(H, 1:N) + + E₀ = expectation_value(XΩ, H_w) + z = E₀ + 0.5 + 0.3im + predicted = 1 / (z - E₀) + + @testset "resolvent, formulation $flav" for flav in (Galerkin(), LeastSquares()) + solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12) + alg = DMRGSolve(; formulation = flav, solver, tol = 1.0e-10, verbosity = 0) + x, = linsolve(XΩ, H_w, XΩ, alg; a₀ = z, a₁ = -1) + @test dot(XΩ, x) ≈ predicted atol = 1.0e-6 + end +end From 38d36f89619e170b4ded0bb9d0aa5403cd75a6f1 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 27 Jul 2026 19:03:43 -0400 Subject: [PATCH 7/7] refactor: propagator/DynamicalDMRG on top of linsolve `corvector.jl` carried its own DMRG sweep, which had drifted from the one in the `linsolve` category on three counts: it converged on the per-site update norm ||AC' - AC|| (which cannot tell converged from stalled) rather than on the residual of the linear system, it used a fixed-tolerance GMRES rather than the structure-selected adaptive solver, and it was single-site only. `propagator` is now a single method for both flavours whose body is a `linsolve!` call. Both flavours keep their variational problem unchanged. `NaiveInvert` is a `Galerkin` solve of `(z - H)x = |psi0>`. `Jeckelmann` still minimizes functional (14) of Jeckelmann2002 and reconstructs G(z) through equation (11): its operator is assembled as `LinearCombination((H^2, H), (1, -2w))` over a hand-built `LazyLincoCache` -- hand-built so that `squaredenvs`' boundary conditions, which matter for `WindowMPS`, survive -- and its stationarity condition `((w - H)^2 + eta^2) psi = -eta |psi0>` is scaled by `-1/eta` so that the right-hand side is `|psi0>` itself. No scaled copy of the right-hand side is needed and `x` is the identical vector; the shift comes out real, so a real initial guess stays real. Consequent changes to the defaults of `DynamicalDMRG`: - `tol` defaults to `Defaults.tol` rather than `10 * Defaults.tol`, and now measures the relative residual of the linear system. - `solver` defaults to `Defaults.alg_linsolve()` (adaptive tolerances), chosen per flavour through an explicit keyword constructor. - a new `trscheme` field prepends a truncated two-site sweep, so the bond dimension no longer has to be fixed by the initial guess. - iteration logs are labelled `linsolve` instead of `DDMRG`. - the `Jeckelmann` flavour throws for real `z` instead of returning `NaN`; the equation-(11) reconstruction divides by `imag(z)`. The docstrings claimed that the second return value is `1/(z - H)|psi0>` for both flavours. That is only true for `NaiveInvert`; for `Jeckelmann` it is the minimizer of the functional, i.e. the imaginary part of the propagator. Say so. Adds a `linsolve` section to the algorithms manual: the category was exported but appeared nowhere in the docs. Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/changelog.md | 33 +++ docs/src/man/algorithms.md | 25 +++ src/algorithms/propagator/corvector.jl | 265 ++++++++++++------------- test/algorithms/dynamical_dmrg.jl | 14 ++ 4 files changed, 199 insertions(+), 138 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 7472c9cfb..38ef42e87 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -32,6 +32,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec Unlike `TDVP` it has no backward-in-time substep (stable for imaginary-time evolution), and passing a truncating `trunc` enables rank-adaptivity (the bond dimension grows and shrinks automatically to track entanglement). +- `linsolve`/`linsolve!` solve the linear system `(a₀ + a₁·A)·x = b` for a finite MPS `x`, with `A` + an MPO/MPOHamiltonian and `b` an MPS. + The shift follows the convention of `KrylovKit.linsolve` and is applied implicitly. + Algorithms are the DMRG-style sweeps `DMRGSolve` (single-site) and `DMRGSolve2` (two-site, + rank-adaptive), each parameterized by a local formulation — `Galerkin` (solve the effective system + directly) or `LeastSquares` (normal equations via squared environments). + Convergence is the local Galerkin residual `‖(a₀ + a₁·A)·x − b‖ / ‖b‖`, and the local solves use + adaptive tolerances by default. +- `DynamicalDMRG` gained a `trunc` field: when supplied, a truncated two-site sweep is prepended, + making `propagator` bond-adaptive instead of fixing the bond dimension of the initial guess. ### Changed @@ -41,6 +51,22 @@ When releasing a new version, move the "Unreleased" changes to a new version sec or the decaying weight in imaginary time). Previously imaginary-time evolution always renormalized every step; **to recover that behavior, pass `normalize = true`** (e.g. for ground-state or thermal-state search via imaginary-time evolution). +- `propagator`/`DynamicalDMRG` are now implemented on top of `linsolve`. + Both flavours keep their variational problem unchanged — `NaiveInvert` is a `Galerkin` solve of + `(z − H)·x = |ψ₀⟩`, and `Jeckelmann` still minimizes functional (14) of Jeckelmann2002 and + reconstructs `G(z)` through equation (11) of that paper. + Three defaults changed as a consequence: convergence is measured by the relative residual of the + linear system rather than by the per-site update norm, `tol` defaults to `Defaults.tol` rather than + `10·Defaults.tol`, and the local solver defaults to `Defaults.alg_linsolve()` (adaptive tolerances) + rather than a fixed-tolerance `GMRES`. + Iteration logs are labelled `linsolve` instead of `DDMRG`. + `propagator` with the `Jeckelmann` flavour now throws for real `z` instead of returning `NaN` + (the equation-(11) reconstruction divides by `imag(z)`). +- The `LeastSquares` formulation of `linsolve` now forms the `A†b` term of the normal equations from + the mixed sandwich `⟨x|A|b⟩` instead of materializing `A·b` as an MPS. + This avoids building a state of bond dimension `χ_A·χ_b` together with its environments, and makes + the formulation applicable to `WindowMPS`, for which no `*(::WindowMPOHamiltonian, ::WindowMPS)` + exists. - `environments` now follows a single positional contract for every state and operator kind: `environments(below, operator, above, alg)`, where `alg` is the environment algorithm (slot 4). The operator form requires an explicit `above`. Auxiliary inputs are keyword-only: @@ -90,6 +116,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec and the right virtual leg of `mpo[end]` and contracted the two — which at length 1 is the *same* tensor, so it returned `O * O` on twice the physical space instead of `O`. ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484)) +- `Defaults.alg_linsolve(; ishermitian = true)` returned a `MINRES` solver, for which KrylovKit has no + `linsolve` method, so the hermitian-indefinite path of `linsolve(…; ishermitian = true)` threw a + `MethodError`. It now falls back to `GMRES`, mirroring KrylovKit's own auto-selection. +- A local `linsolve` solve no longer warns about non-convergence when it undershot its own adaptive + tolerance but still reached the accuracy the outer sweep needs. + Adaptive tolerances routinely dip below the round-off floor of the local problem, which made the + warning fire on nearly every tightly-converged sweep. ### Performance diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 02b078b61..0e3ecc10b 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -339,6 +339,28 @@ state, by a new state. approximate ``` +## `linsolve` + +Solving a linear system ``(a₀ + a₁ A) x = b`` for an MPS ``x``, with ``A`` an MPO and ``b`` an MPS. +The shift ``a₀ + a₁ A`` follows the convention of `KrylovKit.linsolve` and is applied implicitly, so +no shifted operator is ever formed. +The sweeps are DMRG-like: [`DMRGSolve`](@ref) keeps the bond dimension of the initial guess fixed, +while [`DMRGSolve2`](@ref) updates two sites at a time and is therefore rank-adaptive. +Each is parameterized by how the local subproblem is posed — [`Galerkin`](@ref) solves the effective +system directly, [`LeastSquares`](@ref) minimizes ``\|(a₀ + a₁ A) x - b\|²`` through the normal +equations. + +The resolvent is the canonical application: ``\frac{1}{z - H}|ψ₀⟩`` is `linsolve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1)`. + +```@docs; canonical=false +linsolve +linsolve! +DMRGSolve +DMRGSolve2 +Galerkin +LeastSquares +``` + ## Varia What follows is a medley of lesser known (or used) algorithms and don't entirely fit under @@ -350,6 +372,9 @@ Dynamical DMRG has been described in other papers and is a way to find the propa basic idea is that to calculate ``G(z) = ⟨ V | (H-z)^{-1} | V ⟩ `` , one can variationally find ``(H-z) |W ⟩ = | V ⟩ `` and then the propagator simply equals ``G(z) = ⟨ V | W ⟩``. +This is a thin wrapper around [`linsolve`](@ref) — see there for the sweep, the convergence +criterion and the local solver options. + ```@docs; canonical=false propagator DynamicalDMRG diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index bf7057b6d..73f96cf5e 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -8,50 +8,6 @@ abstract type DDMRG_Flavour end """ $(TYPEDEF) -A dynamical DMRG method for calculating dynamical properties and excited states, based on a -variational principle for dynamical correlation functions. - -# Fields - -$(TYPEDFIELDS) - -# See also - -Used as the `algorithm` argument of [`propagator`](@ref). - -# References - -* [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002) -""" -@kwdef struct DynamicalDMRG{F <: DDMRG_Flavour, S} <: Algorithm - "flavour of the algorithm to use, either of type [`NaiveInvert`](@ref) or [`Jeckelmann`](@ref)" - flavour::F = NaiveInvert() - "algorithm used for the linear solvers" - solver::S = Defaults.linearsolver - "tolerance for convergence criterium" - tol::Float64 = Defaults.tol * 10 - "maximal amount of iterations" - maxiter::Int = Defaults.maxiter - "setting for how much information is displayed" - verbosity::Int = Defaults.verbosity -end - -""" - propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ) - -Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG -algorithm. - -# Returns - -- `g`: approximation of the propagator matrix element ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` -- `ψ`: MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩`` -""" -function propagator end - -""" -$(TYPEDEF) - An alternative approach to the dynamical DMRG algorithm, without quadratic terms but with a less controlled approximation. This algorithm minimizes the following cost function @@ -67,50 +23,6 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac """ struct NaiveInvert <: DDMRG_Flavour end -function propagator( - A::AbstractFiniteMPS, z::Number, H, - alg::DynamicalDMRG{NaiveInvert}; init = copy(A) - ) - h_envs = environments(init, H, init) # environments for h - mixedenvs = environments(init, A) # environments for - - ϵ = 2 * alg.tol - log = IterLog("DDMRG") - - LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ) - for iter in 1:(alg.maxiter) - ϵ = 0.0 - - for i in [1:(length(A) - 1); length(A):-1:2] - tos = AC_projection(i, init, A, mixedenvs) - - H_AC = AC_hamiltonian(i, init, H, init, h_envs) - AC = init.AC[i] - AC′, convhist = KrylovKit.linsolve(H_AC, -tos, AC, alg.solver, -z, one(z)) - - ϵ = max(ϵ, norm(AC′ - AC)) - init.AC[i] = AC′ - - convhist.converged == 0 && - @warn "propagator ($i) failed to converge: normres = $(convhist.normres)" - end - - if ϵ <= alg.tol - @infov 2 logfinish!(log, iter, ϵ) - break - end - if iter == alg.maxiter - @warnv 1 logcancel!(log, iter, ϵ) - else - @infov 3 logiter!(log, iter, ϵ) - end - end - end - - return dot(A, init), init -end - """ $(TYPEDEF) @@ -124,10 +36,12 @@ which attains its minimum at ((ω - H)^2 + η^2)|ψ⟩ = -η|ψ₀⟩ ``` -Together with equation (11) from that same paper we can determine the full propagator -``\\frac{1}{z - H}|ψ₀⟩``. +The solution of that equation is the imaginary part of the propagator; together with equation (11) +from that same paper it determines the full ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩``. Because of that +reconstruction step this flavour requires ``η = \\mathrm{Im}(z) ≠ 0``. -Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``. +Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and the minimizer ``|ψ⟩`` of the functional +above (*not* ``\\frac{1}{z - H}|ψ₀⟩`` itself). # See also @@ -139,61 +53,136 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac """ struct Jeckelmann <: DDMRG_Flavour end -function propagator( - A::AbstractFiniteMPS, z::Number, H, - alg::DynamicalDMRG{Jeckelmann}; init = copy(A) +# default local linear solver per flavour, following the structure of the effective operator: +# `z - H` is non-hermitian for complex `z`, while the Jeckelmann operator is hermitian but (after the +# rescaling below) indefinite. Both currently resolve to GMRES, since KrylovKit has no `linsolve` +# method for `MINRES` yet; declaring the structure means Jeckelmann picks it up once it does. +_ddmrg_solver(::NaiveInvert) = Defaults.alg_linsolve() +_ddmrg_solver(::Jeckelmann) = Defaults.alg_linsolve(; ishermitian = true) + +""" +$(TYPEDEF) + +A dynamical DMRG method for calculating dynamical properties and excited states, based on a +variational principle for dynamical correlation functions. + +This is a thin wrapper around [`linsolve`](@ref): the sweep, the convergence criterion and the +adaptive local tolerances are those of [`DMRGSolve`](@ref) / [`DMRGSolve2`](@ref). + +# Fields + +$(TYPEDFIELDS) + +# See also + +Used as the `algorithm` argument of [`propagator`](@ref). + +# References + +* [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002) +""" +struct DynamicalDMRG{F <: DDMRG_Flavour, S, T} <: Algorithm + "flavour of the algorithm to use, either of type [`NaiveInvert`](@ref) or [`Jeckelmann`](@ref)" + flavour::F + "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances (the default)" + solver::S + "tolerance for convergence criterium, measured as the relative residual of the linear system" + tol::Float64 + "maximal amount of iterations" + maxiter::Int + "setting for how much information is displayed" + verbosity::Int + "if supplied, a truncated two-site sweep ([`DMRGSolve2`](@ref)) is prepended to adapt the bond dimension" + trunc::T +end +function DynamicalDMRG(; + flavour = NaiveInvert(), solver = _ddmrg_solver(flavour), tol = Defaults.tol, + maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, trunc = nothing ) - ω = real(z) + return DynamicalDMRG(flavour, solver, tol, maxiter, verbosity, trunc) +end + +# mirrors `_default_linsolve_algorithm`: a loose two-site pass to grow the bond dimension, then a +# single-site polish at the requested tolerance +function _ddmrg_algorithm(alg::DynamicalDMRG) + (; solver, tol, maxiter, verbosity) = alg + alg_1site = DMRGSolve(; solver, tol, maxiter, verbosity) + isnothing(alg.trunc) && return alg_1site + return DMRGSolve2(; + solver, tol = min(1.0e-2, 100tol), maxiter, verbosity, alg.trunc + ) & alg_1site +end + +# The linear system `(a₀ + a₁·A)·x = |ψ₀⟩` solved by each flavour. +# +# `NaiveInvert` is the resolvent itself, `(z - H)·x = |ψ₀⟩`. +# +# `Jeckelmann` is functional (14)'s stationarity condition `((ω - H)² + η²)·ψ = -η|ψ₀⟩`. Scaling that +# equation by `-1/η` puts it in the `(a₀ + a₁·A)·x = b` form with `b = |ψ₀⟩` itself, so no scaled +# copy of the right-hand side is needed and `x` is the very same vector as before. The operator +# `A = H² - 2ω·H` is assembled as a `LinearCombination`, whose `AC_hamiltonian` method reproduces +# exactly the local operator this algorithm used to build by hand. +_ddmrg_shift(::NaiveInvert, z) = (z, -one(z)) +function _ddmrg_shift(::Jeckelmann, z) η = imag(z) + iszero(η) && throw( + ArgumentError( + "`Jeckelmann` requires `imag(z) != 0`; use `NaiveInvert` flavour for real `z`" + ) + ) + return (-abs2(z) / η, -inv(η)) +end - envs1 = environments(init, H, init) # environments for h - H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2 - mixedenvs = environments(init, A) # environments for - - ϵ = 2 * alg.tol - log = IterLog("DDMRG") - - LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ) - for iter in 1:(alg.maxiter) - ϵ = 0.0 - - for i in [1:(length(A) - 1); length(A):-1:2] - tos = AC_projection(i, init, A, mixedenvs) - H1_AC = AC_hamiltonian(i, init, H, init, envs1) - H2_AC = AC_hamiltonian(i, init, H2, init, envs2) - H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1)) - AC′, convhist = KrylovKit.linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1) - - ϵ = max(ϵ, norm(AC′ - init.AC[i])) - init.AC[i] = AC′ - - convhist.converged == 0 && - @warn "propagator ($i) failed to converge: normres $(convhist.normres)" - end - - if ϵ <= alg.tol - @infov 2 logfinish!(log, iter, ϵ) - break - end - if iter == alg.maxiter - @warnv 1 logcancel!(log, iter, ϵ) - else - @infov 3 logiter!(log, iter, ϵ) - end - end - end +_ddmrg_operator(::NaiveInvert, H, x, envs, z) = (H, envs) +function _ddmrg_operator(::Jeckelmann, H, x, envs, z) + H², envs² = squaredenvs(x, H, envs) + A = LinearCombination((H², H), (one(real(z)), -2 * real(z))) + return A, LazyLincoCache(A, (envs², envs)) +end - a = dot(AC_projection(1, init, A, mixedenvs), init.AC[1]) - cb = leftenv(envs1, 1, A) * TransferMatrix(init.AL, H[1:length(A.AL)], A.AL) +# `G(z)` from the solution vector +_ddmrg_value(::NaiveInvert, x, ψ₀, z, H, envs) = dot(ψ₀, x) +function _ddmrg_value(::Jeckelmann, x, ψ₀, z, H, envs) + ω, η = real(z), imag(z) + # equation (11) of Jeckelmann2002: the solve only fixes the imaginary part of the propagator, + # the real part follows from ⟨ψ₀|H|x⟩ + a = dot(ψ₀, x) + cb = leftenv(envs, 1, ψ₀) * TransferMatrix(x.AL, H[1:length(ψ₀.AL)], ψ₀.AL) b = zero(a) for i in 1:length(cb) - b += @plansor cb[i][1 2; 3] * init.C[end][3; 4] * - rightenv(envs1, length(A), A)[i][4 2; 5] * conj(A.C[end][1; 5]) + b += @plansor cb[i][1 2; 3] * x.C[end][3; 4] * + rightenv(envs, length(ψ₀), ψ₀)[i][4 2; 5] * conj(ψ₀.C[end][1; 5]) end + return b / η - ω / η * a + 1im * a +end - v = b / η - ω / η * a + 1im * a - return v, init +""" + propagator(ψ₀::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG; init = ψ₀) -> (g, ψ) + +Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG +algorithm. + +# Returns + +- `g`: approximation of the propagator matrix element ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` +- `ψ`: for [`NaiveInvert`](@ref), the MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``; for + [`Jeckelmann`](@ref), the vector its functional optimizes, + ``-η[(ω - H)^2 + η^2]^{-1}|ψ₀⟩``, i.e. the imaginary part of the propagator, from which `g` is + reconstructed. + +`init` is used as the initial guess and is left untouched. The underlying variational problem is +solved with [`linsolve`](@ref); for full control over the sweep, call that directly, e.g. +`linsolve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1)` for the [`NaiveInvert`](@ref) flavour. +""" +function propagator( + ψ₀::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG; init = ψ₀ + ) + a₀, a₁ = _ddmrg_shift(alg.flavour, z) + x = _promote_state(copy(init), a₀, a₁) + envs = environments(x, H, x) + A, Aenvs = _ddmrg_operator(alg.flavour, H, x, envs, z) + x, = linsolve!(x, A, ψ₀, _ddmrg_algorithm(alg), Aenvs; a₀, a₁) + return _ddmrg_value(alg.flavour, x, ψ₀, z, H, envs), x end function squaredenvs( diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl index 8e07757f4..e16677f76 100644 --- a/test/algorithms/dynamical_dmrg.jl +++ b/test/algorithms/dynamical_dmrg.jl @@ -31,6 +31,20 @@ verbosity_conv = 1 end @test data ≈ predicted atol = 1.0e-8 end + + # `trunc` prepends a two-site sweep, so the bond dimension of the initial guess need not + # already be large enough: start from χ = 2 and let it grow back to the χ = 10 of the solution + @testset "trunc (adaptive χ), flavour $f" for f in (Jeckelmann(), NaiveInvert()) + alg = DynamicalDMRG(; + flavour = f, verbosity = 0, tol = 1.0e-8, trunc = truncrank(16) + ) + init = FiniteMPS(L, ℙ^2, ℙ^2) + data = map(vals) do v + result, = propagator(gs, v + eta, H, alg; init) + return result + end + @test data ≈ predicted atol = 1.0e-6 + end end