feat: standalone linsolve primitives

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
 

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

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__()

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]

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})(

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 &&

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)"

src/algorithms/linsolve/flinsolve.jl

@@ -0,0 +1,223 @@
+# 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 `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
+
+# 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
+
+# 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, 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_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, 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 = _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_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, 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_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, 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 = _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_tol)
+    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)
+    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")
+
+    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, warn_tol
+                )
+                ϵ = 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)
+    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
+    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], warn_tol
+                )
+                ϵ = 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], warn_tol
+                )
+                ϵ = 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