diff --git a/Project.toml b/Project.toml index e6d0a29..8e9d519 100644 --- a/Project.toml +++ b/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 , and contributors"] +version = "2.0.0" [deps] Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e" diff --git a/src/GsvdInitialization.jl b/src/GsvdInitialization.jl index e7954fb..b8e02c3 100644 --- a/src/GsvdInitialization.jl +++ b/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} diff --git a/test/runtests.jl b/test/runtests.jl index e8335cf..5c50339 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -95,6 +95,12 @@ svdX = load_svd_of_gt() @test_throws "must be at most the initial number of components" gsvdnmf(X, Wfit, Hfit, (f.U, f.S, f.V); n2 = 11) @test_throws "must be at most size(W0, 2)" gsvdrecover(X, Wfit, Hfit, 6, (f.U, f.S, f.V)) + # A factorization that is overcomplete relative to X makes the GSVD + # ill-posed; the rank check fires before the factor slices are used. + Xlow = Wfit * Hfit # exact rank 5 + W6, H6 = rand(rng, 30, 6), rand(rng, 6, 20) + @test_throws "numerically rank deficient" gsvdrecover(Xlow, W6, H6, 1, svd(Xlow)) + # `alg` must reach the polishing `nnmf`. An unknown algorithm name is # rejected by `nnmf` only if `alg` is forwarded; the polishing call is the # sole `nnmf` invocation in this four-argument path. @@ -150,8 +156,10 @@ end # be the largest array per call. @testset "non-AbstractArray X" begin rng = StableRNG(3) - U = rand(rng, 10, 3) - V = rand(rng, 3, 8) + # Rank 4 ≥ size(W0, 2): `gsvdrecover` rejects an SVD that is numerically + # rank deficient over the components the factorization claims. + U = rand(rng, 10, 4) + V = rand(rng, 4, 8) Xdense = U * V Xfact = MockFactored(U, V) W0, H0 = rand(rng, 10, 4), rand(rng, 4, 8) @@ -190,6 +198,58 @@ end @test Hd_j ≈ Hf_j end +@testset "appending strategy" begin + rng = StableRNG(7) + # Plant a rank-3 X but hand `gsvdrecover` only 2 of the 3 ground-truth + # components; `appending` must supply the missing one without touching the + # two it was given. + Wgt, Hgt = rand(rng, 12, 3), rand(rng, 3, 9) + X = Wgt * Hgt + W0, H0 = Wgt[:, 1:2], Hgt[1:2, :] + fs = svd(X) + f = (fs.U, fs.S, fs.V) + Wa, Ha, Λ = gsvdrecover(GsvdInitialization.appending, X, copy(W0), copy(H0), 1, f) + @test size(Wa, 2) == 3 && size(Ha, 1) == 3 + # The existing factors are returned exactly, not merely approximately. + @test Wa[:, 1:2] == W0 + @test Ha[1:2, :] == H0 + @test all(>=(0), Wa) && all(>=(0), Ha) + # The appended component improves the fit. + @test sum(abs2, X - Wa * Ha) < sum(abs2, X - W0 * H0) + # `Λ` covers all candidates, aligned with `Hadd` in descending order. + @test length(Λ) == 2 + @test issorted(Λ; rev = true) + + # Factored X agrees with dense X (the strategy needs only `*` and + # `sum(abs2, ·)` from X). + Xfact = MockFactored(Wgt, Hgt) + Wf, Hf, _ = gsvdrecover(GsvdInitialization.appending, Xfact, copy(W0), copy(H0), 1, f) + @test Wf ≈ Wa + @test Hf ≈ Ha + + # "At most kadd": when the factorization is already complete, no candidate + # survives the Λ threshold and amplitude refit, and the factors are + # returned with their original widths. + Wc, Hc, _ = gsvdrecover(GsvdInitialization.appending, X, copy(Wgt), copy(Hgt), 1, svd(X)) + @test size(Wc, 2) == 3 && size(Hc, 1) == 3 + @test Wc == Wgt + @test Hc == Hgt + + # `appending(thresh)` curries the Λ threshold: an impossible threshold + # rejects even genuinely missing components (`Inf > Inf` is false, so this + # also covers exact rank deficiencies). + Wt, Ht, _ = gsvdrecover(GsvdInitialization.appending(Inf), X, copy(W0), copy(H0), 1, f) + @test Wt == W0 && Ht == H0 + + # Eltype is preserved. + W32, H32 = rand(rng, Float32, 12, 4), rand(rng, Float32, 4, 9) + X32 = W32 * H32 + fs32 = svd(X32) + Wa32, Ha32, _ = gsvdrecover(GsvdInitialization.appending, X32, copy(W32), copy(H32), 2, (fs32.U, fs32.S, fs32.V)) + @test eltype(Wa32) === Float32 + @test eltype(Ha32) === Float32 +end + @testset "eltype genericity" begin # The default (`truncating`) path must preserve the input eltype rather than # silently promoting to `Float64`. @@ -238,8 +298,8 @@ end @test isapprox(r_default.H, r_explicit.H; rtol = 1.0e-6) # do-block form: anonymous strategy that simply forwards to `truncating` - r_doblock, _ = gsvdnmf(X, 9 => 10; alg = :cd) do X0, W0, H0, Hadd - GsvdInitialization.truncating(X0, W0, H0, Hadd) + r_doblock, _ = gsvdnmf(X, 9 => 10; alg = :cd) do X0, W0, H0, Hadd, Λ, kadd + GsvdInitialization.truncating(X0, W0, H0, Hadd, Λ, kadd) end @test isapprox(r_doblock.W, r_default.W; rtol = 1.0e-6) @test isapprox(r_doblock.H, r_default.H; rtol = 1.0e-6) @@ -299,8 +359,9 @@ end @test_throws msg gsvdrecover(Xo, W, H, 1, f) @test_throws msg gsvdrecover(X, W, H, 1, fo) Hadd = rand(rng, 1, 8) - @test_throws msg GsvdInitialization.truncating(X, Wo, H, Hadd) - @test_throws msg GsvdInitialization.joint_nnls(X, Wo, H, Hadd) + @test_throws msg GsvdInitialization.truncating(X, Wo, H, Hadd, [Inf], 1) + @test_throws msg GsvdInitialization.joint_nnls(X, Wo, H, Hadd, [Inf], 1) + @test_throws msg GsvdInitialization.appending(X, Wo, H, Hadd, [Inf], 1) # Lazy wrappers carry no axis shift; they must reproduce plain-input # results.