Strategy API: pass all candidates, Λ, and kadd

Project.toml

@@ -1,7 +1,7 @@
 name = "GsvdInitialization"
 uuid = "2ac24108-be9c-42b8-8d78-6a4f62a87e7d"
-authors = ["youdongguo <1010705897@qq.com> and contributors"]
-version = "1.1.2"
+authors = ["youdongguo <1010705897@qq.com>, Tim Holy <tim.holy@gmail.com>, and contributors"]
+version = "2.0.0"
 
 [deps]
 Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"

src/GsvdInitialization.jl

@@ -22,7 +22,7 @@ export gsvdnmf,
     gsvdrecover
 
 @static if VERSION >= v"1.11"
-    eval(Meta.parse("public truncating, joint_nnls"))
+    eval(Meta.parse("public appending, truncating, joint_nnls"))
 end
 
 """
@@ -42,10 +42,17 @@ See [`gsvdrecover`](@ref) for the augmentation step alone, without polishing.
 
 # Arguments
 
-- `strategy`: a callable `(X, W0, H0, Hadd) -> (W_augmented, H_augmented)`
-  that assembles fully non-negative augmented factors from the ranked candidate
-  directions `Hadd`. Defaults to [`GsvdInitialization.truncating`](@ref);
-  [`GsvdInitialization.joint_nnls`](@ref) is the alternative bundled strategy.
+- `strategy`: a callable
+  `(X, W0, H0, Hadd, Λ, kadd) -> (W_augmented, H_augmented)` that assembles
+  fully non-negative augmented factors from the candidate directions (see
+  [`gsvdrecover`](@ref) for the full contract). Defaults to
+  [`GsvdInitialization.truncating`](@ref);
+  [`GsvdInitialization.joint_nnls`](@ref) and
+  [`GsvdInitialization.appending`](@ref) are the alternative bundled
+  strategies. The polishing run is sized to the number of components the
+  strategy actually kept, so a strategy that adds fewer than `n2 - n1` (e.g.
+  [`appending`](@ref)) yields a final factorization with fewer than `n2`
+  components.
 - `X`: non-negative data matrix.
 - `W`: left factor of the existing factorization, of size `(m, n1)`.
 - `H`: right factor of the existing factorization, of size `(n1, p)`.
@@ -101,7 +108,9 @@ function gsvdnmf(
     if alg == :multmse
         W_recover, H_recover = max.(W_recover, truncmult), max.(H_recover, truncmult)
     end
-    result_recover = nnmf(X, n2; kwargs..., alg, init = :custom, tol = tol_final, W0 = copy(W_recover), H0 = copy(H_recover))
+    # An "at most kadd" strategy (e.g. `appending`) may keep fewer than
+    # `n2 - n1` components; size the polishing run to what it kept.
+    result_recover = nnmf(X, size(W_recover, 2); kwargs..., alg, init = :custom, tol = tol_final, W0 = copy(W_recover), H0 = copy(H_recover))
     return result_recover, Λ
 end
 gsvdnmf(X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; kwargs...) =
@@ -185,24 +194,35 @@ generalized SVD between `f` and the current factorization and ranked by
 generalized singular value; `strategy` then assembles the non-negative
 augmented factors.
 
-Return `W_augmented` (`W0` with `kadd` extra columns), `H_augmented` (`H0`
-with `kadd` extra rows), and `Λ`, the generalized singular values that ranked
-the candidate directions.
+Return `W_augmented` (`W0` with the new columns appended), `H_augmented`
+(`H0` with the new rows appended), and `Λ`, the generalized singular values
+of all `n1` candidate directions, in descending order.
 
 # Arguments
 
-- `strategy`: a callable `(X, W0, H0, Hadd) -> (W_augmented, H_augmented)`
-  that assembles fully non-negative augmented factors from the ranked candidate
-  directions `Hadd` (a matrix of `kadd` rows). Defaults to
-  [`GsvdInitialization.truncating`](@ref);
-  [`GsvdInitialization.joint_nnls`](@ref) is the alternative bundled strategy.
+- `strategy`: a callable
+  `(X, W0, H0, Hadd, Λ, kadd) -> (W_augmented, H_augmented)` that assembles
+  fully non-negative augmented factors from the candidate directions. `Hadd`
+  contains *all* `n1` candidates as rows, aligned with `Λ`: `Λ[i]` is the
+  generalized singular value of `Hadd[i, :]`, in descending order. `λ ≈ 1`
+  marks a direction on which the data and the factorization agree; the
+  interesting candidates are those with `λ` far from 1 (`λ ≫ 1`: present in
+  the data but missing from the factorization, `Inf` meaning absent entirely;
+  `λ ≪ 1`: present in the factorization but absent from the data). `kadd` is
+  forwarded from the caller and each strategy defines its interpretation: the
+  bundled [`GsvdInitialization.truncating`](@ref) (the default) and
+  [`GsvdInitialization.joint_nnls`](@ref) add *exactly* `kadd` components,
+  while [`GsvdInitialization.appending`](@ref) adds *at most* `kadd`. The
+  number of components a strategy kept can be read off the width of the
+  returned factors.
 - `X`: non-negative data matrix. `X` need not be an `AbstractMatrix`: any
   object supporting the operations required by the chosen strategy (see
   [`truncating`](@ref) and [`joint_nnls`](@ref)) can be used, e.g. a lazy
   low-rank representation.
 - `W0`: left factor of the existing factorization, of size `(m, n1)`.
 - `H0`: right factor of the existing factorization, of size `(n1, p)`.
-- `kadd`: the number of components to add; must satisfy `1 ≤ kadd ≤ n1`.
+- `kadd`: the requested number of new components, forwarded to `strategy`
+  (see above for its interpretation); must satisfy `1 ≤ kadd ≤ n1`.
 - `f`: a singular value decomposition of `X` with at least `n1` components,
   e.g. from `LinearAlgebra.svd` or `TSVD.tsvd`. Any object whose factors `U`,
   `S`, `V` are indexable as `f[1]`, `f[2]`, `f[3]` works.
@@ -241,20 +261,33 @@ function gsvdrecover(strategy, X, W0::AbstractMatrix, H0::AbstractMatrix, kadd::
     # succeed on the wrong columns; reject it before slicing.
     Base.require_one_based_indexing(U0, S0, V0)
     U0, S0, V0 = U0[:, 1:n], S0[1:n], V0[:, 1:n]
-    Hadd, Λ = init_H(U0, S0, V0, W0, H0, kadd)
-    W, H = strategy(X, W0, H0, Hadd)
+    # The GSVD below treats `Diagonal(S0)` as nonsingular; a numerically
+    # rank-deficient slice means `X` cannot support `n` components and the
+    # candidate directions degenerate.
+    S0[n] > sqrt(eps(float(real(eltype(S0))))) * S0[1] || throw(
+        ArgumentError(
+            "the supplied SVD is numerically rank deficient over its first $n components " *
+                "(s[$n]/s[1] = $(S0[n] / S0[1])): the factorization is overcomplete relative to X, " *
+                "and GSVD augmentation is not defined"
+        )
+    )
+    Hadd, Λ = init_H(U0, S0, V0, W0, H0)
+    ndeficient = count(isinf, Λ)
+    kadd >= ndeficient || @warn "kadd ($kadd) is less than the rank deficiency of W0*H0 ($ndeficient)."
+    W, H = strategy(X, W0, H0, Hadd, Λ, kadd)
     return W, H, Λ
 end
 gsvdrecover(X, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer, f) =
     gsvdrecover(truncating, X, W0, H0, kadd, f)
 
 """
-    truncating(X, W0, H0, Hadd) -> (W_augmented, H_augmented)
+    truncating(X, W0, H0, Hadd, Λ, kadd) -> (W_augmented, H_augmented)
 
-Default [`gsvdrecover`](@ref) strategy. Nonnegative least-squares (NNLS) is
-used only for the rescaling weights `α` of the existing columns; the new
-columns of `W` are computed by ordinary least squares and made non-negative by
-an NNDSVD step, after which all columns are rebalanced.
+Default [`gsvdrecover`](@ref) strategy; adds **exactly** `kadd` components,
+the leading (largest-`Λ`) candidate directions. Nonnegative least-squares
+(NNLS) is used only for the rescaling weights `α` of the existing columns; the
+new columns of `W` are computed by ordinary least squares and made
+non-negative by an NNDSVD step, after which all columns are rebalanced.
 
 This strategy requires only `*` and `sum(abs2, ·)` from `X`, so `X` may be a
 lazy or factored low-rank representation rather than a materialized matrix.
@@ -263,10 +296,10 @@ and the rescaling jointly, at greater cost.
 
 Return non-negative `(W_augmented, H_augmented)`.
 """
-function truncating(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix)
+function truncating(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix, Λ::AbstractVector, kadd::Integer)
     X isa AbstractArray && Base.require_one_based_indexing(X)
     Base.require_one_based_indexing(W0, H0, Hadd)
-    kadd = size(Hadd, 1)
+    Hadd = Hadd[1:kadd, :]
     Wadd, a = init_W(X, W0, H0, Hadd)
     Wadd_nn, Hadd_nn = nndsvd(X, kadd, initdata = (U = Wadd, S = ones(eltype(Wadd), kadd), V = Hadd'))
     W0_1, H0_1 = [a' .* W0 Wadd_nn], [H0; Hadd_nn]
@@ -276,10 +309,87 @@ function truncating(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMat
 end
 
 """
-    joint_nnls(X, W0, H0, Hadd) -> (W_augmented, H_augmented)
+    appending(X, W0, H0, Hadd, Λ, kadd) -> (W_augmented, H_augmented)
+    appending(thresh; rtol = nothing) -> strategy
+
+[`gsvdrecover`](@ref) strategy that holds the existing factorization fixed:
+the returned factors contain `W0` and `H0` unmodified, with new components
+appended. It adds **at most** `kadd` components, applying three rejection
+criteria to the candidates (considered in descending-`Λ` order):
+
+- `Λ > thresh`: a direction with `λ ≤ thresh` does not carry enough excess
+  energy in the data, relative to what the factorization already explains, to
+  justify a new component. `appending` itself uses `thresh = 1`;
+  `appending(thresh)` returns a strategy with a stricter (or looser) value.
+- nonzero fitted amplitude: a candidate whose refit amplitude is zero cannot
+  reduce the residual.
+- relative energy above `rtol`: `Λ` is a ratio, so a direction of negligible
+  absolute energy can still have a huge `λ` (e.g. exactly-rank-deficient data
+  under an overcomplete factorization). A kept candidate's fitted contribution
+  must satisfy `‖β·w*h'‖ > rtol·‖X‖` (Frobenius). The default
+  `rtol = √eps(eltype)` rejects numerical dust only.
+
+The new columns of `W` are computed by ordinary least squares on the residual
+`X - W0*H0` (i.e., with the existing-column weights pinned at 1) and made
+non-negative together with the candidate directions by an NNDSVD step, after
+which their amplitudes are refit by NNLS (existing columns again pinned at 1).
+
+Use this strategy when the existing components must not be perturbed — for
+example when they are externally constrained, already validated, or shared
+with a fit over a larger domain. The bundled alternatives [`truncating`](@ref)
+and [`joint_nnls`](@ref) instead rescale the existing columns to rebalance the
+augmented factorization as a whole.
+
+Like [`truncating`](@ref), this strategy requires only `*` and `sum(abs2, ·)`
+from `X`, so `X` may be a lazy or factored low-rank representation.
+
+Return non-negative `(W_augmented, H_augmented)`.
+"""
+appending(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix, Λ::AbstractVector, kadd::Integer) =
+    _appending(1, nothing, X, W0, H0, Hadd, Λ, kadd)
+appending(thresh::Real; rtol::Union{Real, Nothing} = nothing) =
+    (X, W0, H0, Hadd, Λ, kadd) -> _appending(thresh, rtol, X, W0, H0, Hadd, Λ, kadd)
+
+function _appending(thresh::Real, rtol, X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix, Λ::AbstractVector, kadd::Integer)
+    X isa AbstractArray && Base.require_one_based_indexing(X)
+    Base.require_one_based_indexing(W0, H0, Hadd)
+    # `Λ` is in descending order, so the candidates above threshold are the
+    # leading ones.
+    nsel = min(kadd, count(>(thresh), Λ))
+    nsel == 0 && return W0, H0
+    Hadd = Hadd[1:nsel, :]
+    Wadd, _ = init_W(X, W0, H0, Hadd; α = ones(eltype(W0), size(W0, 2)))
+    Wadd_nn, Hadd_nn = nndsvd(X, nsel, initdata = (U = Wadd, S = ones(eltype(Wadd), nsel), V = Hadd'))
+    Wadd_nn, Hadd_nn = abs.(Wadd_nn), abs.(Hadd_nn)
+    # The NNDSVD step retains only one sign quadrant of each rank-1 term, which
+    # perturbs its amplitude.  Refit the amplitudes β of just the new columns
+    # (existing columns pinned at 1): min_{β ≥ 0} ‖X − W0*H0 − Σₖ βₖ wₖ hₖ‖².
+    B = (Wadd_nn' * Wadd_nn) .* (Hadd_nn * Hadd_nn')
+    a = diag(Wadd_nn' * X * Hadd_nn') - diag((Wadd_nn' * W0) * (H0 * Hadd_nn'))
+    if isposdef(Symmetric(B))
+        β = vec(nonneg_lsq(B, a; alg = :fnnls, gram = true))
+        # Reject candidates that cannot reduce the residual (β = 0) or whose
+        # fitted contribution is negligible at the scale of the data:
+        # ‖βₖ wₖ hₖ'‖²_F = βₖ² B[k,k].
+        rt = rtol === nothing ? sqrt(eps(float(real(eltype(B))))) : rtol
+        floor2 = rt^2 * sum(abs2, X)
+        sel = findall(k -> β[k]^2 * B[k, k] > floor2, eachindex(β))
+        Wadd_nn = β[sel]' .* Wadd_nn[:, sel]
+        Hadd_nn = Hadd_nn[sel, :]
+    else
+        # Degenerate Gram matrix (e.g. a zero or duplicated candidate column):
+        # keep the NNDSVD amplitudes rather than risk a singular solve.
+        sum(abs2, B) <= 1.0e-12 || @warn "B is not positive definite; keeping NNDSVD amplitudes for the appended components." maxlog = 1
+    end
+    return [W0 Wadd_nn], [H0; Hadd_nn]
+end
+
+"""
+    joint_nnls(X, W0, H0, Hadd, Λ, kadd) -> (W_augmented, H_augmented)
 
-Alternative [`gsvdrecover`](@ref) strategy that solves for the new columns of
-`W` and the rescaling `α` of the existing columns jointly, as a single
+Alternative [`gsvdrecover`](@ref) strategy that adds **exactly** `kadd`
+components (the leading, largest-`Λ` candidates), solving for the new columns
+of `W` and the rescaling `α` of the existing columns jointly, as a single
 nonnegative least-squares (NNLS) problem. `Hadd` is first projected onto the
 non-negative orthant, keeping whichever sign of each candidate direction
 better matches the non-negative part of the residual `X - W0*H0`.
@@ -292,34 +402,35 @@ also requires `X - W0*H0` and `eltype(X)`.
 
 Return non-negative `(W_augmented, H_augmented)`.
 """
-function joint_nnls(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix)
+function joint_nnls(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix, Λ::AbstractVector, kadd::Integer)
     X isa AbstractArray && Base.require_one_based_indexing(X)
     Base.require_one_based_indexing(W0, H0, Hadd)
-    Hadd_nn = truncatepos(Hadd', X, W0, H0)'
+    Hadd_nn = truncatepos(Hadd[1:kadd, :]', X, W0, H0)'
     Wadd, a = init_W_joint_nnls(X, W0, H0, Hadd_nn)
     W0_1, H0_1 = [a' .* W0 Wadd], [H0; Hadd_nn]
     return abs.(W0_1), abs.(H0_1)
 end
 
-function init_H(U0::AbstractMatrix, S0::AbstractVector, V0::AbstractMatrix, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer)
+function init_H(U0::AbstractMatrix, S0::AbstractVector, V0::AbstractMatrix, W0::AbstractMatrix, H0::AbstractMatrix)
     _, _, Q, D1, D2, R = svd(Matrix(Diagonal(S0)), (U0' * W0) * (H0 * V0))
     r0 = size(U0, 2)
     k = findfirst(x -> x != 0, D2[1, :])
     k = (k === nothing) ? r0 : k - 1
-    kadd >= k || @warn "kadd ($kadd) is less than the rank deficiency of W0*H0 ($k)."
     F = (diag(D1[(k + 1):r0, (k + 1):r0]) ./ diag(D2[1:(r0 - k), (k + 1):r0])) .^ 2
     Λ = vcat(fill(Inf, k), F)
-    H_index = sortperm(Λ, rev = true)[1:kadd]
-    # Columns of inv(R*Q') = Q*inv(R) selected by H_index, via a triangular
-    # backsolve on just those right-hand sides: the GSVD's Q is orthogonal and
-    # R is square upper triangular (Diagonal(S0) is nonsingular, so k+l = r0).
-    E = zeros(eltype(R), size(R, 1), length(H_index))
+    H_index = sortperm(Λ, rev = true)
+    # Columns of inv(R*Q') = Q*inv(R) in descending-Λ order, via a triangular
+    # backsolve: the GSVD's Q is orthogonal and R is square upper triangular
+    # (Diagonal(S0) is nonsingular, so k+l = r0).
+    E = zeros(eltype(R), r0, r0)
     for (j, idx) in enumerate(H_index)
         E[idx, j] = 1
     end
     Hadd = Q * (UpperTriangular(R) \ E)
     Hadd_1 = V0 * Hadd
-    return Hadd_1', Λ
+    # The outputs are aligned: both follow descending-Λ order, so `Λ[i]` is
+    # the generalized singular value of `Hadd[i, :]`.
+    return Hadd_1', Λ[H_index]
 end
 
 function init_W_joint_nnls(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}) where {T}