From 266e94454907c359bb546356513c7a7edc9a7fa1 Mon Sep 17 00:00:00 2001 From: leburgel Date: Sun, 4 Oct 2026 10:05:51 +0200 Subject: [PATCH 1/4] Solve the truncated svd and eigh pullbacks by conjugate gradients --- src/MatrixAlgebraKit.jl | 1 + src/common/pullbacks.jl | 36 -------- src/common/stein.jl | 176 ++++++++++++++++++++++++++++++++++++++++ src/pullbacks/eigh.jl | 15 ++-- src/pullbacks/svd.jl | 24 ++++-- 5 files changed, 203 insertions(+), 49 deletions(-) create mode 100644 src/common/stein.jl diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 78d36bfe5..9fa6f6644 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -92,6 +92,7 @@ include("common/defaults.jl") include("common/householder.jl") include("common/initialization.jl") include("common/pullbacks.jl") +include("common/stein.jl") include("common/safemethods.jl") include("common/view.jl") include("common/regularinv.jl") diff --git a/src/common/pullbacks.jl b/src/common/pullbacks.jl index c99f03859..43bba4c45 100644 --- a/src/common/pullbacks.jl +++ b/src/common/pullbacks.jl @@ -36,42 +36,6 @@ iterating over `ind`, so that this also works for an `ind` that lives on a devic is_leading_index(ind::AbstractRange, p::Int) = ind == 1:p is_leading_index(ind::AbstractVector, p::Int) = length(ind) == p && all(ind .== 1:p) -""" - accelerative_smith_iteration!(X, Xₙ, G, w, atol, maxiter) - -Solve `X = B + G * X * Diagonal(w)` by summing the Neumann series -`X = Σₖ Gᵏ * B * Diagonal(w)ᵏ` by doubling (Smith's method), i.e. by repeatedly adding -`G^(2ʲ) * X * Diagonal(w)^(2ʲ)` to `X` until the norm of that increment drops below `atol`, -for at most `maxiter` steps. - -On entry, `X` contains `B`, and it is overwritten with the result. `Xₙ` is used as a buffer, -and `G` and `w` are overwritten. `w` is normalized such that `maximum(abs, w) == 1`, so that -squaring it can only shrink it; `G` is scaled by the inverse factor to compensate. - -Reference: https://doi.org/10.1016/j.aml.2009.01.012. -""" -function accelerative_smith_iteration!(X, Xₙ, G, w, atol, maxiter) - Gₙ = similar(G) - wmax = maximum(abs, w) - w ./= wmax - G .*= wmax - for k in 1:maxiter - Xₙ = rmul!(mul!(Xₙ, G, X), Diagonal(w)) - if maximum(abs, Xₙ) < atol - break - end - X .+= Xₙ - if k == maxiter - @warn "Sylvester iteration did not converge after $k iterations, final norm of X: $(maximum(abs, X))" - break - end - w .= w .^ 2 - Gₙ = mul!(Gₙ, G, G) - G, Gₙ = Gₙ, G - end - return X -end - """ antihermitian_columns!(X, ind) diff --git a/src/common/stein.jl b/src/common/stein.jl new file mode 100644 index 000000000..e5734cb77 --- /dev/null +++ b/src/common/stein.jl @@ -0,0 +1,176 @@ +# Solvers for the Stein equation X - G X Diagonal(w) = B, i.e. (1 - wᵢ G) xᵢ = bᵢ per column, of the +# truncated pullbacks. svd: G = P Pᴴ or Pᴴ P (the smaller), wᵢ = 1/σᵢ²; eigh: G = P, wᵢ = 1/λᵢ; with +# P the part of A outside the kept vectors and γᵢ < 1 the spectral radius of wᵢ G. + +""" + accelerative_smith_iteration!(X, Xₙ, G, w, atol, maxiter) + +Solve `X = B + G * X * Diagonal(w)` by summing the Neumann series +`X = Σₖ Gᵏ * B * Diagonal(w)ᵏ` by doubling (Smith's method), i.e. by repeatedly adding +`G^(2ʲ) * X * Diagonal(w)^(2ʲ)` to `X` until the norm of that increment drops below `atol`, +for at most `maxiter` steps. + +On entry, `X` contains `B`, and it is overwritten with the result. `Xₙ` is used as a buffer, +and `G` and `w` are overwritten. `w` is normalized such that `maximum(abs, w) == 1`, so that +squaring it can only shrink it; `G` is scaled by the inverse factor to compensate. + +Reference: https://doi.org/10.1016/j.aml.2009.01.012. +""" +function accelerative_smith_iteration!(X, Xₙ, G, w, atol, maxiter) + Gₙ = similar(G) + wmax = maximum(abs, w) + w ./= wmax + G .*= wmax + for k in 1:maxiter + Xₙ = rmul!(mul!(Xₙ, G, X), Diagonal(w)) + if maximum(abs, Xₙ) < atol + break + end + X .+= Xₙ + if k == maxiter + @warn "Sylvester iteration did not converge after $k iterations, final norm of X: $(maximum(abs, X))" + break + end + w .= w .^ 2 + Gₙ = mul!(Gₙ, G, G) + G, Gₙ = Gₙ, G + end + return X +end + +# smallest Ritz value: smallest eigenvalue of the Lanczos tridiagonal from the CG coefficients `α`, `β` +function _cg_ritz_min(α, β) + l = length(α) + d = similar(α) + d[1] = 1 / α[1] + for j in 2:l + d[j] = 1 / α[j] + β[j - 1] / α[j - 1] + end + e = sqrt.(β[1:(l - 1)]) ./ α[1:(l - 1)] + return LinearAlgebra.eigmin(LinearAlgebra.SymTridiagonal(d, e)) +end + +# (1 - wᵢ G) applied to the columns of Z +_stein_op(G, w, Z) = Z .- (G * Z) .* transpose(w) + +# iterations until the squared residual norm `r` drops below `tol²`, at the faster of its rates over +# the last `l` iterations (from `r₀`) and all `m` (from `rᵢ`) +function _cg_remaining(r, r₀, rᵢ, tol, l, m) + q = min(log(r / r₀) / l, log(r / rᵢ) / m) # logarithm of the reduction per iteration + return q < 0 ? log(tol^2 / r) / q : oftype(float(r), Inf) +end + +# real parts of the column-wise inner products of A and B, on the device of A and B. CPU arrays use +# `dot` per column, which avoids the temporary. This is restricted to `Matrix` rather than +# `StridedMatrix`, which GPU arrays also are: there it would return a CPU vector, one `dot` call +# each, which the GPU broadcasts in `hermitian_stein_cg!` cannot mix with the device arrays. +_coldots(A, B) = vec(real(sum(conj.(A) .* B; dims = 1))) +const _CPUMatrix = Union{Matrix, SubArray{<:Any, 2, <:Matrix}} +_coldots(A::_CPUMatrix, B::_CPUMatrix) = [real(LinearAlgebra.dot(view(A, :, j), view(B, :, j))) for j in axes(A, 2)] + +""" + hermitian_stein_cg!(X, applyG!, formG, w, atol, maxiter; cost_apply, cost_apply_formed, cost_form, cost_square = nothing) + +Solve `X - G * X * Diagonal(w) = B` for Hermitian `G` by conjugate gradients on all columns at +once, in at most `maxiter` iterations: O(n² k) per iteration for k columns, against O(n³) per +doubling step. Only the products wᵢ G enter, so neither needs normalizing. `X` contains `B` on +entry and is overwritten with the result. + +`applyG!(Y, Z)` sets `Y = G * Z` without forming `G`, and `formG()` returns `G`. `G` is formed +once the predicted remaining applications save more than `cost_form`, given the costs per column +`cost_apply` (by `applyG!`) and `cost_apply_formed` (by `G`); `cost_form = 0` forms it at once. +With `cost_square`, for indefinite `G` (eigh), the solver restarts on the residual equation +multiplied by 1 + wᵢ G, (1 - wᵢ² G²) xᵢ = (1 + wᵢ G) bᵢ, once half the predicted remaining cost +exceeds `cost_square` (forming G² and the new right-hand side): its spectrum [1 - γᵢ², 1] needs +about half the iterations. A column stops once its residual is below `atol` times the smallest +Ritz value of the slowest column, so that its error is below about `atol`. +""" +function hermitian_stein_cg!( + X, applyG!, formG, w, atol, maxiter; + cost_apply, cost_apply_formed = cost_apply, cost_form, cost_square = nothing, nprobe::Int = 5 + ) + RT = real(eltype(X)) + G = iszero(cost_form) ? formG() : nothing + tol = RT(atol) # stopping residual 2-norm, refined by the Ritz values + ρ = _coldots(X, X) + cols = findall(ρ .> tol^2) # active columns, kept contiguous at the front + nact = length(cols) + iszero(nact) && return fill!(X, zero(eltype(X))) + R = X[:, cols] # X keeps B until the end + P = zero(R) + Q = similar(R) + Xc = zero(R) + wc = w[cols] + ρ = ρ[cols] + β = zero(ρ) + ρ₀ = copy(ρ) # squared residual norms at the previous check + ρᵢ = copy(ρ) # and at the start + αs = [RT[] for _ in 1:nact] + βs = [RT[] for _ in 1:nact] + for numiter in 1:maxiter + if numiter > 1 && (numiter - 1) % nprobe == 0 + c = argmax(abs.(view(wc, 1:nact))) # the slowest column + tol = atol * _cg_ritz_min(αs[c], βs[c]) + nact = _cg_compact!( + view(ρ, 1:nact), tol, view(R, :, 1:nact), view(P, :, 1:nact), view(Xc, :, 1:nact), view(wc, 1:nact), + view(β, 1:nact), view(ρ₀, 1:nact), view(ρᵢ, 1:nact), view(cols, 1:nact), view(αs, 1:nact), view(βs, 1:nact) + ) + iszero(nact) && break + # predicted column applications + napply = sum(_cg_remaining.(view(ρ, 1:nact), view(ρ₀, 1:nact), view(ρᵢ, 1:nact), tol, nprobe, numiter - 1)) + if isnothing(G) && napply * (cost_apply - cost_apply_formed) > cost_form + G = formG() + end + if !isnothing(cost_square) && napply * cost_apply / 2 > cost_square # restart on the squared equation + isnothing(G) && (G = formG()) + X₁ = zero(X) # the current solution + X₁[:, cols] .= Xc + R = X .- _stein_op(G, w, X₁) + X .= R .+ (G * R) .* transpose(w) + G² = G * G + hermitian_stein_cg!(X, nothing, () -> G², w .^ 2, atol, maxiter; cost_apply, cost_form = 0, nprobe) + return X .+= X₁ + end + ρ₀ .= ρ + end + Pₐ, Rₐ, Qₐ, wₐ = view(P, :, 1:nact), view(R, :, 1:nact), view(Q, :, 1:nact), view(wc, 1:nact) + ρₐ, βₐ = view(ρ, 1:nact), view(β, 1:nact) + Pₐ .= Rₐ .+ Pₐ .* transpose(βₐ) + isnothing(G) ? applyG!(Qₐ, Pₐ) : mul!(Qₐ, G, Pₐ) + Qₐ .= Pₐ .- Qₐ .* transpose(wₐ) # q = (1 - wᵢ G) p + α = ρₐ ./ _coldots(Pₐ, Qₐ) + view(Xc, :, 1:nact) .+= Pₐ .* transpose(α) + Rₐ .-= Qₐ .* transpose(α) + βₐ .= ρₐ # ρold + ρₐ .= _coldots(Rₐ, Rₐ) + βₐ .= ρₐ ./ βₐ + for (j, a, b) in zip(1:nact, Array(α), Array(βₐ)) + push!(αs[j], a) + push!(βs[j], b) + end + nact = _cg_compact!( + ρₐ, tol, Rₐ, Pₐ, view(Xc, :, 1:nact), wₐ, + βₐ, view(ρ₀, 1:nact), view(ρᵢ, 1:nact), view(cols, 1:nact), view(αs, 1:nact), view(βs, 1:nact) + ) + iszero(nact) && break + end + iszero(nact) || @warn "conjugate gradients did not converge in $maxiter iterations, largest residual norm: $(sqrt(maximum(view(ρ, 1:nact))))" + fill!(X, zero(eltype(X))) + X[:, cols] .= Xc + return X +end + +# move the columns with `ρ > tol²` to the front of all arrays (in order) and return their number +function _cg_compact!(ρ, tol, arrays...) + keep = Array(ρ) .> tol^2 + all(keep) && return length(keep) + perm = vcat(findall(keep), findall(.!keep)) + for a in (ρ, arrays...) + if a isa AbstractMatrix + a .= a[:, perm] + else + a .= a[perm] + end + end + return count(keep) +end diff --git a/src/pullbacks/eigh.jl b/src/pullbacks/eigh.jl index 2b3afacd0..95e044e96 100755 --- a/src/pullbacks/eigh.jl +++ b/src/pullbacks/eigh.jl @@ -143,16 +143,16 @@ function eigh_trunc_pullback!( ΔA::AbstractMatrix, A, DV, ΔDV; degeneracy_atol::Real = default_pullback_rank_atol(DV[1]), gauge_atol::Real = default_pullback_gauge_atol(ΔDV[2]), - maxiter::Int = 100 # TODO: better default, depending on expected number of steps using quadratic convergence? + maxiter::Int = 10 * size(ΔA, 1) # conjugate-gradient iterations ) # Basic size checks and determination Dmat, V = DV - (n, p) = size(V) + (n, k) = size(V) D = diagview(Dmat) - p == length(D) || throw(DimensionMismatch()) + k == length(D) || throw(DimensionMismatch()) (n, n) == size(ΔA) || throw(DimensionMismatch()) - iszero(p) && return ΔA + iszero(k) && return ΔA ΔDmat, ΔV = ΔDV VᴴΔAV, ΔV₊ = check_and_prepare_eigh_cotangents( @@ -162,13 +162,16 @@ function eigh_trunc_pullback!( if !iszerotangent(ΔV₊) X₀ = rdiv!(ΔV₊, Diagonal(D)) AP = mul!(copy(A), V * Dmat, V', -1, 1) - X = accelerative_smith_iteration!(X₀, similar(X₀), AP, inv.(D), degeneracy_atol, maxiter) + X = hermitian_stein_cg!( + X₀, nothing, () -> AP, inv.(D), degeneracy_atol, maxiter; + cost_apply = n^2, cost_form = 0, cost_square = n^3 + n^2 * k + ) Z .+= X # we cannot directly multiply Z * V' into ΔA, because we have to # take the Hermitian part, and cannot apply project_hermitian! to # the current contents of ΔA # TODO: add an `add_project_hermitian!` - # recycle AP's storage, but overwrite it: `accelerative_smith_iteration!` may leave a power of AP in it + # recycle AP's storage ΔA′ = project_hermitian!(mul!(AP, Z, V')) ΔA .+= ΔA′ else diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl index 576997b69..8b4b9ea26 100755 --- a/src/pullbacks/svd.jl +++ b/src/pullbacks/svd.jl @@ -223,17 +223,17 @@ function svd_trunc_pullback!( rank_atol::Real = 0, degeneracy_atol::Real = default_pullback_rank_atol(USVᴴ[2]), gauge_atol::Real = default_pullback_gauge_atol(ΔUSVᴴ...), - maxiter::Int = 100 # TODO: better default, depending on expected number of steps using quadratic convergence? + maxiter::Int = 10 * minimum(size(ΔA)) # conjugate-gradient iterations ) # Extract the SVD components U, Smat, Vᴴ = USVᴴ m, n = size(U, 1), size(Vᴴ, 2) (m, n) == size(ΔA) || throw(DimensionMismatch(lazy"size of ΔA ($(size(ΔA))) does not match size of USVᴴ ($m, $n)")) S = diagview(Smat) - p = length(S) - p == size(U, 2) || throw(DimensionMismatch(lazy"U has $p columns but S has $(length(S)) singular values")) - p == size(Vᴴ, 1) || throw(DimensionMismatch(lazy"Vᴴ has $p rows but S has $(length(S)) singular values")) - iszero(p) && return ΔA + k = length(S) + k == size(U, 2) || throw(DimensionMismatch(lazy"U has $k columns but S has $(length(S)) singular values")) + k == size(Vᴴ, 1) || throw(DimensionMismatch(lazy"Vᴴ has $k rows but S has $(length(S)) singular values")) + iszero(k) && return ΔA # Extract and check the cotangents ΔU, ΔSmat, ΔVᴴ = ΔUSVᴴ @@ -257,7 +257,12 @@ function svd_trunc_pullback!( if m ≤ n X = rmul!(AP * Y₀ᴴ', Diagonal(S⁻¹)) X .+= X₀ - X = accelerative_smith_iteration!(X, X₀, AP * AP', S⁻¹ .^ 2, degeneracy_atol, maxiter) # recycle X₀ + APᴴZ = similar(X, n, k) # for applying AP AP' without forming it + X = hermitian_stein_cg!( + X, (GZ, Z) -> mul!(GZ, AP, mul!(view(APᴴZ, :, axes(Z, 2)), AP', Z)), () -> AP * AP', + S⁻¹ .^ 2, degeneracy_atol, maxiter; + cost_apply = 2 * m * n, cost_apply_formed = m^2, cost_form = m^2 * n + ) Yᴴ = lmul!(Diagonal(S⁻¹), X' * AP) Yᴴ .+= Y₀ᴴ ΔA = mul!(ΔA, X, Vᴴ, 1, 1) @@ -265,7 +270,12 @@ function svd_trunc_pullback!( else Y = rmul!(AP' * X₀, Diagonal(S⁻¹)) Y .+= Y₀ᴴ' - Y = accelerative_smith_iteration!(Y, similar(Y), AP' * AP, S⁻¹ .^ 2, degeneracy_atol, maxiter) + APZ = similar(Y, m, k) # for applying AP' AP without forming it + Y = hermitian_stein_cg!( + Y, (GZ, Z) -> mul!(GZ, AP', mul!(view(APZ, :, axes(Z, 2)), AP, Z)), () -> AP' * AP, + S⁻¹ .^ 2, degeneracy_atol, maxiter; + cost_apply = 2 * m * n, cost_apply_formed = n^2, cost_form = n^2 * m + ) X = rmul!(AP * Y, Diagonal(S⁻¹)) X .+= X₀ ΔA = mul!(ΔA, X, Vᴴ, 1, 1) From 9220744924bb8487169cce2687a936802dd5881a Mon Sep 17 00:00:00 2001 From: Lander Burgelman <39218680+leburgel@users.noreply.github.com> Date: Wed, 7 Oct 2026 15:54:29 +0200 Subject: [PATCH 2/4] Update src/common/stein.jl Co-authored-by: Jutho --- src/common/stein.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/common/stein.jl b/src/common/stein.jl index e5734cb77..2f74dc6ef 100644 --- a/src/common/stein.jl +++ b/src/common/stein.jl @@ -1,6 +1,7 @@ -# Solvers for the Stein equation X - G X Diagonal(w) = B, i.e. (1 - wᵢ G) xᵢ = bᵢ per column, of the -# truncated pullbacks. svd: G = P Pᴴ or Pᴴ P (the smaller), wᵢ = 1/σᵢ²; eigh: G = P, wᵢ = 1/λᵢ; with -# P the part of A outside the kept vectors and γᵢ < 1 the spectral radius of wᵢ G. +# Solvers for the Stein equation X - G X Diagonal(w) = B, i.e. (1 - wᵢ G) xᵢ = bᵢ per column. +# Used in the pullback of truncated decompositions, with G = P Pᴴ or Pᴴ P, wᵢ = 1/σᵢ² for the case of SVD, +# and G = P, wᵢ = 1/λᵢ for the case of EIG(H). Here, P the part of A outside the kept vectors. +# Naive iteration requires γᵢ, the spectral radius of wᵢ G, to be smaller than 1. """ accelerative_smith_iteration!(X, Xₙ, G, w, atol, maxiter) From 591657a69670d6fa2737aae348cb3cb1b9d6f594 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 7 Oct 2026 16:28:42 +0200 Subject: [PATCH 3/4] Undo rogue `p` -> `k` rename --- src/pullbacks/eigh.jl | 8 ++++---- src/pullbacks/svd.jl | 12 ++++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/pullbacks/eigh.jl b/src/pullbacks/eigh.jl index 95e044e96..6f8063f58 100755 --- a/src/pullbacks/eigh.jl +++ b/src/pullbacks/eigh.jl @@ -148,11 +148,11 @@ function eigh_trunc_pullback!( # Basic size checks and determination Dmat, V = DV - (n, k) = size(V) + (n, p) = size(V) D = diagview(Dmat) - k == length(D) || throw(DimensionMismatch()) + p == length(D) || throw(DimensionMismatch()) (n, n) == size(ΔA) || throw(DimensionMismatch()) - iszero(k) && return ΔA + iszero(p) && return ΔA ΔDmat, ΔV = ΔDV VᴴΔAV, ΔV₊ = check_and_prepare_eigh_cotangents( @@ -164,7 +164,7 @@ function eigh_trunc_pullback!( AP = mul!(copy(A), V * Dmat, V', -1, 1) X = hermitian_stein_cg!( X₀, nothing, () -> AP, inv.(D), degeneracy_atol, maxiter; - cost_apply = n^2, cost_form = 0, cost_square = n^3 + n^2 * k + cost_apply = n^2, cost_form = 0, cost_square = n^3 + n^2 * p ) Z .+= X # we cannot directly multiply Z * V' into ΔA, because we have to diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl index 8b4b9ea26..8fa8904dc 100755 --- a/src/pullbacks/svd.jl +++ b/src/pullbacks/svd.jl @@ -230,10 +230,10 @@ function svd_trunc_pullback!( m, n = size(U, 1), size(Vᴴ, 2) (m, n) == size(ΔA) || throw(DimensionMismatch(lazy"size of ΔA ($(size(ΔA))) does not match size of USVᴴ ($m, $n)")) S = diagview(Smat) - k = length(S) - k == size(U, 2) || throw(DimensionMismatch(lazy"U has $k columns but S has $(length(S)) singular values")) - k == size(Vᴴ, 1) || throw(DimensionMismatch(lazy"Vᴴ has $k rows but S has $(length(S)) singular values")) - iszero(k) && return ΔA + p = length(S) + p == size(U, 2) || throw(DimensionMismatch(lazy"U has $p columns but S has $(length(S)) singular values")) + p == size(Vᴴ, 1) || throw(DimensionMismatch(lazy"Vᴴ has $p rows but S has $(length(S)) singular values")) + iszero(p) && return ΔA # Extract and check the cotangents ΔU, ΔSmat, ΔVᴴ = ΔUSVᴴ @@ -257,7 +257,7 @@ function svd_trunc_pullback!( if m ≤ n X = rmul!(AP * Y₀ᴴ', Diagonal(S⁻¹)) X .+= X₀ - APᴴZ = similar(X, n, k) # for applying AP AP' without forming it + APᴴZ = similar(X, n, p) # for applying AP AP' without forming it X = hermitian_stein_cg!( X, (GZ, Z) -> mul!(GZ, AP, mul!(view(APᴴZ, :, axes(Z, 2)), AP', Z)), () -> AP * AP', S⁻¹ .^ 2, degeneracy_atol, maxiter; @@ -270,7 +270,7 @@ function svd_trunc_pullback!( else Y = rmul!(AP' * X₀, Diagonal(S⁻¹)) Y .+= Y₀ᴴ' - APZ = similar(Y, m, k) # for applying AP' AP without forming it + APZ = similar(Y, m, p) # for applying AP' AP without forming it Y = hermitian_stein_cg!( Y, (GZ, Z) -> mul!(GZ, AP', mul!(view(APZ, :, axes(Z, 2)), AP, Z)), () -> AP' * AP, S⁻¹ .^ 2, degeneracy_atol, maxiter; From 965274f04b499ea207bba9894deb23a13011ff9a Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 8 Oct 2026 08:23:36 +0200 Subject: [PATCH 4/4] Remove trailing whitespace --- src/common/stein.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common/stein.jl b/src/common/stein.jl index 2f74dc6ef..e0b579a7c 100644 --- a/src/common/stein.jl +++ b/src/common/stein.jl @@ -1,5 +1,5 @@ # Solvers for the Stein equation X - G X Diagonal(w) = B, i.e. (1 - wᵢ G) xᵢ = bᵢ per column. -# Used in the pullback of truncated decompositions, with G = P Pᴴ or Pᴴ P, wᵢ = 1/σᵢ² for the case of SVD, +# Used in the pullback of truncated decompositions, with G = P Pᴴ or Pᴴ P, wᵢ = 1/σᵢ² for the case of SVD, # and G = P, wᵢ = 1/λᵢ for the case of EIG(H). Here, P the part of A outside the kept vectors. # Naive iteration requires γᵢ, the spectral radius of wᵢ G, to be smaller than 1.