From ad1ea8120abe9c482e33068fa6ca5dd08fb914a5 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 12 Jun 2026 13:45:23 -0500 Subject: [PATCH 1/4] Declare 1-based indexing at public entry points MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The pipeline runs through svd, nndsvd, nnmf, and sparse, all of which assume 1-based indexing. Without a declaration, offset-axes inputs fail deep inside LinearAlgebra — or worse: an offset SVD wider than size(W0, 2) makes the 1:n factor slices in gsvdrecover succeed on the wrong columns, silently returning wrong factors. Add Base.require_one_based_indexing at gsvdnmf, gsvdrecover, and the two bundled strategies (X checked only when it is an AbstractArray, preserving the duck-typed factored-matrix contract). Tests assert that offset inputs error cleanly on every public path and that view-wrapped inputs reproduce plain-input results. Co-Authored-By: Claude Fable 5 --- Project.toml | 4 ++- src/GsvdInitialization.jl | 17 +++++++++++++ test/runtests.jl | 53 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 73 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a8d4032..b892a98 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ Kronecker = "0.5" LinearAlgebra = "1" NMF = "1" NonNegLeastSquares = "0.4" +OffsetArrays = "1" SparseArrays = "1" StableRNGs = "1" TSVD = "0.4" @@ -32,8 +33,9 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" NMF = "6ef6ca0d-6ad7-5ff6-b225-e928bfa0a386" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "ExplicitImports", "NMF", "Test", "FileIO", "JLD2", "StableRNGs"] +test = ["Aqua", "ExplicitImports", "NMF", "Test", "FileIO", "JLD2", "OffsetArrays", "StableRNGs"] diff --git a/src/GsvdInitialization.jl b/src/GsvdInitialization.jl index 8ab33fe..0980d85 100644 --- a/src/GsvdInitialization.jl +++ b/src/GsvdInitialization.jl @@ -38,6 +38,8 @@ Arguments: - `f`: SVD (or Truncated SVD) of `X` +All array arguments, including the factors of `f`, must use 1-based indexing. + Keyword arguments: - `tol_final`: the tolerance of the NMF polishing step, default: 1e-4 @@ -62,6 +64,7 @@ function gsvdnmf(strategy, X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatr alg = :cd, truncmult = 1e-5, kwargs...) + Base.require_one_based_indexing(X, W, H) n1 = size(W, 2) kadd = n2 - n1 kadd > 0 || throw(ArgumentError("The number of components to add must be positive; got n2 = $n2, size(W, 2) = $n1")) @@ -109,6 +112,7 @@ Keyword arguments: Other keyword arguments are passed to `NMF.nnmf`. """ function gsvdnmf(strategy, X::AbstractMatrix, ncomponents::Pair{<:Integer,<:Integer}; tol_final=1e-4, tol_intermediate=tol_final, kwargs...) + Base.require_one_based_indexing(X) n1, n2 = ncomponents f = tsvd(X, n2) W0, H0 = nndsvd(X, n1; initdata = (U = f[1], S = f[2], V = f[3])) @@ -152,13 +156,22 @@ Arguments: `kadd`: number of new components `f`: SVD (or Truncated SVD) of `X` + +All array arguments, including the factors of `f`, must use 1-based indexing. """ function gsvdrecover(strategy, X, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer, f) + # `X` may be a non-array factored representation (see the docstring); only + # arrays carry axes to validate. + X isa AbstractArray && Base.require_one_based_indexing(X) + Base.require_one_based_indexing(W0, H0) _, n = size(W0) kadd > 0 || throw(ArgumentError("kadd must be positive; got $kadd")) kadd <= n || throw(ArgumentError("the number of extra columns must be at most size(W0, 2); got kadd = $kadd, size(W0, 2) = $n")) size(first(f), 2) >= n || throw(ArgumentError("the supplied SVD has $(size(first(f), 2)) components but size(W0, 2) = $n are required")) U0, S0, V0 = f + # An offset-axes SVD wider than `n` would make the `1:n` slices below + # 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) @@ -177,6 +190,8 @@ to solve for the new columns of `W` (i.e., `Wadd`). Returns non-negative `(W_augmented, H_augmented)`. """ function truncating(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix) + X isa AbstractArray && Base.require_one_based_indexing(X) + Base.require_one_based_indexing(W0, H0, Hadd) kadd = size(Hadd, 1) 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')) @@ -200,6 +215,8 @@ projecting `Hadd`). Returns non-negative `(W_augmented, H_augmented)`. """ function joint_nnls(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix) + X isa AbstractArray && Base.require_one_based_indexing(X) + Base.require_one_based_indexing(W0, H0, Hadd) Hadd_nn = truncatepos(Hadd', X, W0, H0)' Wadd, a = init_W_joint_nnls(X, W0, H0, Hadd_nn) W0_1, H0_1 = [a' .* W0 Wadd], [H0; Hadd_nn] diff --git a/test/runtests.jl b/test/runtests.jl index 13a2450..8a233bf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using Aqua using ExplicitImports using LinearAlgebra, NMF, FileIO +using OffsetArrays using StableRNGs @testset "Aqua" begin @@ -236,6 +237,58 @@ end @test isapprox(r_doblock.H, r_default.H; rtol = 1e-6) end +@testset "generic axes" begin + # The pipeline runs through `svd`, `nndsvd`, `nnmf`, and `sparse`, all of + # which assume 1-based indexing, so the public entry points declare that + # assumption with `require_one_based_indexing`. Offset-axes inputs must + # fail at entry with a clear error — without the declaration they fail deep + # inside LinearAlgebra, or worse: an offset SVD wider than `size(W0, 2)` + # makes the `1:n` factor slices succeed on the wrong columns, silently + # returning wrong factors. + rng = StableRNG(8) + W, H = rand(rng, 10, 4), rand(rng, 4, 8) + X = W * H + 0.01 * rand(rng, 10, 8) + fs = svd(X) + f = (fs.U, fs.S, fs.V) + msg = "offset arrays are not supported" + + Xo = OffsetArray(X, -2, -3) + Wo = OffsetArray(W, -2, 0) + Ho = OffsetArray(H, 0, -3) + # Full SVD of X has 8 components > size(W, 2) = 4: the silent-wrong-columns + # shape. + fo = (OffsetArray(fs.U, 0, -1), OffsetArray(fs.S, -1), OffsetArray(fs.V, 0, -1)) + + @test_throws msg gsvdnmf(Xo, 3 => 4; alg = :cd) + @test_throws msg gsvdnmf(X, Wo, Ho, f; n2 = 5) + @test_throws msg gsvdrecover(X, Wo, Ho, 1, f) + @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) + + # Lazy wrappers carry no axis shift; they must reproduce plain-input + # results. + vX, vW, vH = view(X, :, :), view(W, :, :), view(H, :, :) + Wr, Hr, Λr = gsvdrecover(X, copy(W), copy(H), 2, f) + Wv, Hv, Λv = gsvdrecover(vX, vW, vH, 2, f) + @test Wv ≈ Wr + @test Hv ≈ Hr + @test Λv ≈ Λr + Wjr, Hjr, _ = gsvdrecover(GsvdInitialization.joint_nnls, X, copy(W), copy(H), 2, f) + Wjv, Hjv, _ = gsvdrecover(GsvdInitialization.joint_nnls, vX, vW, vH, 2, f) + @test Wjv ≈ Wjr + @test Hjv ≈ Hjr + # `gsvdnmf` runs `nnmf`, so two runs of the same pipeline agree only to the + # ~1e-9 nondeterminism of multithreaded BLAS; `rtol = 1e-6` as in the + # cross-call checks above. + r_plain, _ = gsvdnmf(X, copy(W), copy(H), f; n2 = 5, alg = :cd) + r_view, _ = gsvdnmf(vX, vW, vH, f; n2 = 5, alg = :cd) + @test isapprox(r_plain.W, r_view.W; rtol = 1e-6) + @test isapprox(r_plain.H, r_view.H; rtol = 1e-6) +end + @testset "joint optimize W and alpha" begin rng = StableRNG(7) W = W_GT From c5387e4ac79b4b513d992984ab7db7de42561b28 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 12 Jun 2026 13:55:26 -0500 Subject: [PATCH 2/4] Cover the integer-n2 gsvdnmf methods and :multmse branch These were the only uncovered lines in src/: the convenience methods that compute tsvd(X, n2) internally, and the flooring of the augmented factors to truncmult that lets multiplicative updates (which cannot move entries off zero) polish them. src/ line coverage is now 100%. Co-Authored-By: Claude Fable 5 --- test/runtests.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 8a233bf..359cecb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -237,6 +237,32 @@ end @test isapprox(r_doblock.H, r_default.H; rtol = 1e-6) end +@testset "integer-n2 convenience methods and :multmse" begin + rng = StableRNG(9) + X = rand(rng, 30, 20) + W, H = rand(rng, 30, 4), rand(rng, 4, 20) + # `gsvdnmf(X, W, H, n2)` computes `tsvd(X, n2)` itself and forwards to the + # explicit-`f` method; the explicit-strategy form runs the same pipeline, + # so they agree to the usual cross-call rtol. + r_default, _ = gsvdnmf(X, W, H, 5; alg = :cd) + r_strategy, _ = gsvdnmf(GsvdInitialization.truncating, X, W, H, 5; alg = :cd) + @test size(r_default.W, 2) == 5 + @test isapprox(r_default.W, r_strategy.W; rtol = 1e-6) + @test isapprox(r_default.H, r_strategy.H; rtol = 1e-6) + + # :multmse floors the augmented factors to `truncmult` so multiplicative + # updates (which cannot move entries off zero) can polish them. + r_mult, _ = gsvdnmf(X, W, H, 5; alg = :multmse, maxiter = 10^4) + @test size(r_mult.W, 2) == 5 + @test r_mult.converged + @test all(>=(0), r_mult.W) && all(>=(0), r_mult.H) + # A full-rank random X puts the best rank-5 residual near 0.13, so an + # absolute bound is meaningless; require :multmse to land close to :cd. + res_cd = sum(abs2, X - r_default.W * r_default.H) + res_mult = sum(abs2, X - r_mult.W * r_mult.H) + @test res_mult <= 1.1 * res_cd +end + @testset "generic axes" begin # The pipeline runs through `svd`, `nndsvd`, `nnmf`, and `sparse`, all of # which assume 1-based indexing, so the public entry points declare that From 6ea057eed402411fb22ac468798cf9e0b4150c3f Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 12 Jun 2026 14:27:34 -0500 Subject: [PATCH 3/4] Rewrite docstrings and add doctest coverage Add a module docstring, restructure all public docstrings to standard Julia conventions (prose "Return ..." sentence, bulleted arguments, consistent keyword style), and fix errors: the return type is NMF.Result (not NMF.NMFResult), and ncomponents accepts any Pair{<:Integer,<:Integer}. Document the n1 < n2 <= 2n1 constraint, the integer-for-f convenience form, and the duck-typed contract for X in gsvdrecover. The strategy docstrings now cross-reference each other with cost/requirement tradeoffs. Add jldoctest examples to gsvdnmf (both forms) and gsvdrecover, and run them in the test suite via Documenter.doctest. Co-Authored-By: Claude Fable 5 --- Project.toml | 4 +- src/GsvdInitialization.jl | 218 +++++++++++++++++++++++++------------- test/runtests.jl | 6 ++ 3 files changed, 156 insertions(+), 72 deletions(-) diff --git a/Project.toml b/Project.toml index b892a98..e6d0a29 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" [compat] Aqua = "0.8" +Documenter = "1" ExplicitImports = "1.15" FileIO = "1.18" JLD2 = "0.6" @@ -29,6 +30,7 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -38,4 +40,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "ExplicitImports", "NMF", "Test", "FileIO", "JLD2", "OffsetArrays", "StableRNGs"] +test = ["Aqua", "Documenter", "ExplicitImports", "NMF", "Test", "FileIO", "JLD2", "OffsetArrays", "StableRNGs"] diff --git a/src/GsvdInitialization.jl b/src/GsvdInitialization.jl index 0980d85..77b2cc7 100644 --- a/src/GsvdInitialization.jl +++ b/src/GsvdInitialization.jl @@ -1,3 +1,14 @@ +""" +`GsvdInitialization` augments an existing low-rank non-negative matrix +factorization (NMF) `X ≈ W*H` with additional components, using the +generalized singular value decomposition (GSVD) between the current +factorization and an SVD of `X` to discover what the factorization is missing +(the GSVD-NMF method, [doi:10.1016/j.isci.2026.114708](https://doi.org/10.1016/j.isci.2026.114708)). + +The main entry points are [`gsvdnmf`](@ref), which augments a factorization +and polishes the result with NMF, and [`gsvdrecover`](@ref), which performs +the augmentation step alone. +""" module GsvdInitialization using LinearAlgebra: Diagonal, I, Symmetric, UpperTriangular, cholesky, diag, isposdef, svd, tr @@ -16,47 +27,61 @@ end """ result, Λ = gsvdnmf([strategy,] X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; - n2 = size(first(f), 2), - tol_final=1e-4, - kwargs...) + n2 = size(first(f), 2), tol_final = 1e-4, + alg = :cd, truncmult = 1e-5, kwargs...) -Augment `W` and `H` to have `n2` components, subsequently polished by NMF. +Augment the existing factorization `X ≈ W*H` to `n2` components, then polish +the result with NMF. An integer `n2` may be passed in place of `f`, in which +case a truncated SVD with `n2` components is computed internally. + +Return the `NMF.Result` of the polishing run (its `W` and `H` fields hold the +augmented factors) and `Λ`, the generalized singular values that ranked the +candidate augmentation directions. -Arguments: +See [`gsvdrecover`](@ref) for the augmentation step alone, without polishing. -- `strategy`: a callable `(X, W0, H0, Hadd) -> (W_augmented, H_augmented)` that - produces fully-assembled non-negative augmented factors from the candidate - directions `Hadd` ranked by [`init_H`](@ref). Defaults to - [`GsvdInitialization.truncating`](@ref); [`GsvdInitialization.joint_nnls`](@ref) is - the alternative (joint-NNLS) strategy. +# Arguments -- `X`: non-negative data matrix +- `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. +- `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)`. +- `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. -- `W` and `H`: initial NMF factorization +All array arguments, including the factors of `f`, must use 1-based indexing. -- `n2`: the number of components in augmented factorization +# Keyword arguments -- `f`: SVD (or Truncated SVD) of `X` +- `n2`: the number of components after augmentation; `n2 - n1` must satisfy + `1 ≤ n2 - n1 ≤ n1` (at most a doubling per call). +- `tol_final`: convergence tolerance of the NMF polishing run (default: `1e-4`). +- `alg`: NMF algorithm for the polishing run, forwarded to `NMF.nnmf` + (default: `:cd`). With `alg == :multmse`, the augmented factors are first + floored at `truncmult`, because multiplicative updates require strictly + positive factors. +- `truncmult`: flooring level applied when `alg == :multmse` (default: `1e-5`). -All array arguments, including the factors of `f`, must use 1-based indexing. +Remaining keyword arguments are forwarded to `NMF.nnmf`. -Keyword arguments: +# Examples -- `tol_final`: the tolerance of the NMF polishing step, default: 1e-4 +```jldoctest +julia> using LinearAlgebra: svd -- `alg`: the NMF algorithm for the polishing step, forwarded to `NMF.nnmf`, - default: `:cd`. When `alg == :multmse` the augmented factors are floored to - `truncmult` first, because multiplicative updates require strictly positive - factors. +julia> X = Float64[1 0 0 1 0; 0 1 0 1 1; 0 0 1 0 1; 1 1 1 2 2]; # rank 3 -- `truncmult`: the flooring level applied to the augmented factors when - `alg == :multmse`, default: 1e-5 +julia> W0 = Float64[1 0; 0 1; 0 0; 1 1]; H0 = Float64[1 0 0 1 0; 0 1 0 1 1]; # rank-2 factorization of X -Other keyword arguments are passed to `NMF.nnmf`. +julia> result, Λ = gsvdnmf(X, W0, H0, svd(X); n2 = 3); -Returns the `NMF.NMFResult` from the polishing step (its `W` and `H` fields hold -the augmented factors) and `Λ`, the generalized singular values that ranked the -candidate augmentation directions. +julia> size(result.W), size(result.H) +((4, 3), (3, 5)) +``` """ function gsvdnmf(strategy, X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; n2 = size(first(f), 2), @@ -85,31 +110,51 @@ gsvdnmf(X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, n2::Integer; kw gsvdnmf(truncating, X, W, H, n2; kwargs...) """ - result, Λ = gsvdnmf([strategy,] X::AbstractMatrix, ncomponents::Pair{Int,Int}; tol_final=1e-4, tol_intermediate=1e-4, kwargs...) + result, Λ = gsvdnmf([strategy,] X::AbstractMatrix, ncomponents; + tol_final = 1e-4, tol_intermediate = tol_final, kwargs...) + +Perform GSVD-NMF on the non-negative data matrix `X`: compute an NMF with `n1` +components, augment it to `n2` components, and polish with a final NMF run. -Perform "GSVD-NMF" on the data matrix `X`. +The initial factorization is computed by `NMF.nnmf` with NNDSVD initialization +seeded from a truncated SVD of `X`; the same SVD supplies the augmentation +directions. To augment an existing factorization instead, use the four-argument +[`gsvdnmf`](@ref) method. -Arguments: +Return the `NMF.Result` of the final NMF run (its `W` and `H` fields hold the +factors) and `Λ`, the generalized singular values that ranked the candidate +augmentation directions. + +# Arguments - `strategy`: see the four-argument [`gsvdnmf`](@ref) method. Defaults to [`GsvdInitialization.truncating`](@ref). +- `X`: non-negative data matrix; must use 1-based indexing. +- `ncomponents`: a pair `n1 => n2` of integers requesting augmentation from + `n1` to `n2` components, where `n1 < n2 ≤ 2n1`. An integer `n` is shorthand + for `n-1 => n` (add a single component). + +# Keyword arguments -- `X`: non-negative data matrix +- `tol_final`: convergence tolerance of the final NMF run (default: `1e-4`). +- `tol_intermediate`: convergence tolerance of the initial rank-`n1` NMF run + (default: same as `tol_final`). -- `ncomponents`: in the form of `n1 => n2`, augments from `n1` components to `n2`components, - where `n1` is the number of components for initial NMF (under-complete NMF), and `n2` is the number of - components for final NMF. +Remaining keyword arguments are forwarded to `NMF.nnmf`. -Alternatively, `ncomponents` can be an integer denoting the number of components for final NMF. -In this case, `gsvdnmf` defaults to augment components on initial NMF solution by 1. +# Examples -Keyword arguments: +```jldoctest +julia> X = Float64[1 0 0 1 0; 0 1 0 1 1; 0 0 1 0 1; 1 1 1 2 2]; # rank 3 -- `tol_final`: The tolerance of final NMF, default:`10^{-4}` +julia> result, Λ = gsvdnmf(X, 2 => 3); -- `tol_intermediate`: The tolerance of initial NMF (under-complete NMF), default: tol_final +julia> size(result.W), size(result.H) +((4, 3), (3, 5)) -Other keyword arguments are passed to `NMF.nnmf`. +julia> sum(abs2, X - result.W*result.H) < 1e-6 * sum(abs2, X) # near-exact rank-3 fit +true +``` """ function gsvdnmf(strategy, X::AbstractMatrix, ncomponents::Pair{<:Integer,<:Integer}; tol_final=1e-4, tol_intermediate=tol_final, kwargs...) Base.require_one_based_indexing(X) @@ -130,34 +175,55 @@ gsvdnmf(X::AbstractMatrix, ncomponents_final::Integer; kwargs...) = """ W_augmented, H_augmented, Λ = gsvdrecover([strategy,] X, W0, H0, kadd, f) -Augment components for `W0` and `H0` without polishing by NMF. - -`strategy` is a callable `(X, W0, H0, Hadd) -> (W_augmented, H_augmented)` that -produces fully-assembled non-negative augmented factors from the candidate -directions `Hadd` ranked by [`init_H`](@ref). Defaults to -[`GsvdInitialization.truncating`](@ref); [`GsvdInitialization.joint_nnls`](@ref) is the -alternative bundled strategy. - -Outputs: +Augment the factorization `X ≈ W0*H0` with `kadd` additional components. This +is the augmentation step of [`gsvdnmf`](@ref), without the final NMF polish. + +Candidate directions for the new rows of `H` are extracted from the +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. + +# 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. +- `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`. +- `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. -`W_augmented`, `H_augmented`: the full augmented NMF factors (with `kadd` extra -components appended to `W0`/`H0`) - -`Λ`: generalized singular values used to rank the candidate augmentation directions +All array arguments, including the factors of `f`, must use 1-based indexing. -Arguments: +# Examples -`X`: non-negative 2D data matrix +```jldoctest +julia> using LinearAlgebra: svd -`W0`: NMF solution +julia> X = Float64[1 0 0 1 0; 0 1 0 1 1; 0 0 1 0 1; 1 1 1 2 2]; # rank 3 -`H0`: NMF solution +julia> W0 = Float64[1 0; 0 1; 0 0; 1 1]; H0 = Float64[1 0 0 1 0; 0 1 0 1 1]; # rank-2 factorization of X -`kadd`: number of new components +julia> W, H, Λ = gsvdrecover(X, W0, H0, 1, svd(X)); -`f`: SVD (or Truncated SVD) of `X` +julia> size(W), size(H) +((4, 3), (3, 5)) -All array arguments, including the factors of `f`, must use 1-based indexing. +julia> sum(abs2, X - W*H) < sum(abs2, X - W0*H0) # the new component improves the fit +true +``` """ function gsvdrecover(strategy, X, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer, f) # `X` may be a non-array factored representation (see the docstring); only @@ -183,11 +249,17 @@ gsvdrecover(X, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer, f) = """ truncating(X, W0, H0, Hadd) -> (W_augmented, H_augmented) -Default `gsvdrecover` strategy. Restricts the use of nonnegative least-squares (NNLS) -to the component weights `α`, and uses least-squares followed by an NNDSVD step -to solve for the new columns of `W` (i.e., `Wadd`). +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. + +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. +See [`joint_nnls`](@ref) for an alternative that solves for the new columns +and the rescaling jointly, at greater cost. -Returns non-negative `(W_augmented, H_augmented)`. +Return non-negative `(W_augmented, H_augmented)`. """ function truncating(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix) X isa AbstractArray && Base.require_one_based_indexing(X) @@ -204,15 +276,19 @@ end """ joint_nnls(X, W0, H0, Hadd) -> (W_augmented, H_augmented) -Alternative `gsvdrecover` strategy that jointly solves for the new columns of -`W` and the rescaling `α` of existing columns using nonnegative least-squares -(NNLS). `Hadd` is first projected onto the non-negative orthant. +Alternative [`gsvdrecover`](@ref) strategy that solves 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`. -Beyond the `*` and `sum(abs2, ·)` that the default [`truncating`](@ref) strategy -needs from `X`, this path also requires `X - W*H` and `eltype(X)` (used while -projecting `Hadd`). +The joint NNLS problem has one unknown for every entry of the new columns of +`W` plus one rescaling weight per existing column, so this strategy is more +expensive than the default [`truncating`](@ref), especially when `X` has many +rows. Beyond the `*` and `sum(abs2, ·)` that `truncating` needs from `X`, it +also requires `X - W0*H0` and `eltype(X)`. -Returns non-negative `(W_augmented, H_augmented)`. +Return non-negative `(W_augmented, H_augmented)`. """ function joint_nnls(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMatrix) X isa AbstractArray && Base.require_one_based_indexing(X) diff --git a/test/runtests.jl b/test/runtests.jl index 359cecb..1db5900 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using GsvdInitialization using Test using Aqua +using Documenter using ExplicitImports using LinearAlgebra, NMF, FileIO @@ -19,6 +20,11 @@ end all_qualified_accesses_are_public = VERSION >= v"1.11") end +DocMeta.setdocmeta!(GsvdInitialization, :DocTestSetup, :(using GsvdInitialization); recursive=true) +@testset "Doctests" begin + doctest(GsvdInitialization; manual = false) +end + # Minimal `Factorization` subtype that implements just the matrix products and # `sum(abs2, ·)` that `gsvdrecover` calls on `X`. Used to verify that # `gsvdrecover` and its helpers accept any `X` for which those operations From 30049c7cca68dd4da7336145f8a5caf296eacf98 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 12 Jun 2026 14:46:04 -0500 Subject: [PATCH 4/4] Update citation to published article; modernize README The paper is published as Guo & Holy, iScience 29(3):114708 (2026), doi:10.1016/j.isci.2026.114708. CITATION.cff now titles the software itself and carries the article in a preferred-citation block (arXiv preprint kept as an identifier); README cites the article directly. README: add JuliaHub version badge and an installation section, make the demo runnable from an installed package via pkgdir, and replace the hand-copied "Functions" section with a brief API overview that defers to the doctested docstrings. The demo now uses a finite maxiter and checks `converged`, demonstrating that NMF stops at the tolerance rather than the iteration cap. Co-Authored-By: Claude Fable 5 --- CITATION.cff | 34 ++++++++--- README.md | 167 ++++++++++++++++----------------------------------- 2 files changed, 77 insertions(+), 124 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 8cb0505..e8212de 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,11 +1,8 @@ -# This CITATION.cff file was generated with cffinit. -# Visit https://bit.ly/cffinit to generate yours today! - cff-version: 1.2.0 -title: 'GSVD-NMF: Recovering Missing Features in Non-negative Matrix Factorization' +title: GsvdInitialization.jl message: >- - If you use this software, please cite it using the - metadata from this file. + If you use this software, please cite the article listed under + preferred-citation. type: software authors: - given-names: Youdong @@ -15,10 +12,12 @@ authors: family-names: Holy orcid: 'https://orcid.org/0000-0002-2429-1071' identifiers: + - type: doi + value: 10.1016/j.isci.2026.114708 + description: The journal article describing the method. - type: url value: 'https://arxiv.org/abs/2408.08260' - description: The ArXiv deposit of the encompassing paper. -doi: 'https://doi.org/10.48550/arXiv.2408.08260' + description: The arXiv preprint of the article. repository-code: 'https://github.com/HolyLab/GsvdInitialization.jl' abstract: >- Non-negative matrix factorization (NMF) is an important @@ -35,3 +34,22 @@ abstract: >- experimental results demonstrate that GSVD-NMF often recovers missing features from under-complete NMF and helps NMF achieve better local optima. +preferred-citation: + type: article + title: >- + Recovering missing features in nonnegative matrix factorization via + generalized singular value decomposition + authors: + - given-names: Youdong + family-names: Guo + orcid: 'https://orcid.org/0009-0007-7787-3722' + - given-names: Timothy E. + family-names: Holy + orcid: 'https://orcid.org/0000-0002-2429-1071' + journal: iScience + volume: 29 + issue: 3 + start: 114708 + year: 2026 + month: 3 + doi: 10.1016/j.isci.2026.114708 diff --git a/README.md b/README.md index a8703e6..dc006a4 100644 --- a/README.md +++ b/README.md @@ -3,158 +3,93 @@ [![CI](https://github.com/HolyLab/GsvdInitialization.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/GsvdInitialization.jl/actions/workflows/CI.yml) [![codecov](https://codecov.io/gh/HolyLab/GsvdInitialization.jl/graph/badge.svg?token=LxqRCsZIvn)](https://codecov.io/gh/HolyLab/GsvdInitialization.jl) [![Aqua QA](https://juliatesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) +[![version](https://juliahub.com/docs/General/GsvdInitialization/stable/version.svg)](https://juliahub.com/ui/Packages/General/GsvdInitialization) -This package implements the technique in the paper [GSVD-NMF: Recovering Missing Features in -Non-negative Matrix Factorization](https://arxiv.org/abs/2408.08260). -It is used to recover Non-negative matrix factorization (NMF) components from an initial lower-rank factorization by exploiting the generalized singular value decomposition (GSVD) between existing NMF results and the SVD of X. -This method allows the incremental expansion of the number of components, which can be convenient and effective for interactive analysis of large-scale data. +This package implements GSVD-NMF ([Guo & Holy, *iScience* 2026](https://doi.org/10.1016/j.isci.2026.114708)), a method for recovering missing components in non-negative matrix factorization (NMF). +Starting from a lower-rank factorization `X ≈ W*H`, it proposes new components from the generalized singular value decomposition (GSVD) between the existing factorization and the SVD of `X`, then polishes the augmented factorization with further NMF iterations. +Because components can be added incrementally, GSVD-NMF is convenient and effective for interactive analysis of large-scale data. -See also [NMFMerge](https://github.com/HolyLab/NMFMerge.jl) for the converse operation. Together, the two result in a substantial improvement in the quality and consistency of NMF factorization. +See also [NMFMerge](https://github.com/HolyLab/NMFMerge.jl) for the converse operation (merging redundant components). Together, the two substantially improve the quality and consistency of NMF factorizations. ---------------------------- +## Installation -Demo: +GsvdInitialization is a registered package; type `]` at the `julia>` prompt to enter `pkg>` mode and install it with -To run this demo, NMF.jl and LinearAlgebra.jl are also required. - -Install and load packages (type `]` at the `julia>` prompt to enter `pkg>` mode): -```julia -pkg> add GsvdInitialization; -julia> using GsvdInitialization, NMF, LinearAlgebra; ``` - -Generating ground truth with 10 features. - -```julia -julia> include("demo/generate_ground_truth.jl") -julia> W_GT, H_GT = generate_ground_truth(); -julia> X = W_GT*H_GT; +pkg> add GsvdInitialization ``` -Sample Figure +## Demo -Running standard NMF(HALS) using NNDSVD as initialization on X. Here, we're taking a couple of precautions to try to ensure the best possible result from NMF: -- we disable premature convergence by setting `maxiter` to something that is practically infinite -- we use the full `svd`, rather than `rsvd`, for initializing NNDSVD, as `svd` gives higher-quality results than `rsvd` -Despite these precautions, we'll see that the NMF result leaves much to be desired: +The demo below also uses [NMF.jl](https://github.com/JuliaStats/NMF.jl) and the LinearAlgebra standard library: ```julia -julia> result_hals = nnmf(X, 10; init=:nndsvd, alg = :cd, tol = 1e-4, maxiter=10^12, initdata = svd(X)); -julia> sum(abs2, X-result_hals.W*result_hals.H)/sum(abs2, X) -0.0999994991270576 +julia> using GsvdInitialization, NMF, LinearAlgebra ``` -The result is given by - -Sample Figure -This factorization is not perfect as two components are the same and two features share one component. -Then, running GSVD-NMF on X (also using NNSVD as initialization) and computing the new reconstruction error: +Generate a ground truth with 10 features (the script ships with the package): ```julia -julia> result_gsvd, Λ = gsvdnmf(X, 9=>10; alg = :cd, tol_final = 1e-4, tol_intermediate = 1e-2, maxiter = 10^12); -julia> Wgsvd, Hgsvd = result_gsvd.W, result_gsvd.H; -julia> sum(abs2, X-Wgsvd*Hgsvd)/sum(abs2, X) -1.2322603074132593e-10 -``` - -`Λ` is the vector of generalized singular values that ranked the candidate augmentation directions, useful as a diagnostic for understanding which directions the algorithm chose. -An imperfect factorization from `nnmf` alone was augmented by `gsvdnmf` to a perfect factorization. -Here are the new components: - -Sample Figure - - ---------------------------- - -## Functions - -result, Λ = **gsvdnmf**([strategy,] X::AbstractMatrix, ncomponents::Pair{Int,Int}; - tol_final=1e-4, - tol_intermediate=1e-4, - kwargs...) - -Perform "GSVD-NMF" on the data matrix `X`. - -Arguments: - -- `strategy`: optional augmentation strategy `(X, W0, H0, Hadd) -> (W_aug, H_aug)`. - Defaults to `GsvdInitialization.truncating`; pass `GsvdInitialization.joint_nnls` - for the alternative bundled strategy, or supply your own. - -- `X`: non-negative data matrix +julia> include(joinpath(pkgdir(GsvdInitialization), "demo", "generate_ground_truth.jl")); -- `ncomponents`: in the form of `n1 => n2`, augments from `n1` components to `n2`components, - where `n1` is the number of components for initial NMF (under-complete NMF), and `n2` is the number of - components for final NMF. - -Alternatively, `ncomponents` can be an integer denoting the number of components for final NMF. -In this case, `gsvdnmf` defaults to augment components on initial NMF solution by 1. - -Keyword arguments: - -- `tol_final`: The tolerance of final NMF, default:`10^{-4}` - -- `tol_intermediate`: The tolerance of initial NMF (under-complete NMF), default: tol_final - -Other keyword arguments are passed to `NMF.nnmf`. - ------ - -result, Λ = **gsvdnmf**([strategy,] X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; - n2 = size(first(f), 2), - tol_final=1e-4, - kwargs...) - -Augment `W` and `H` to have `n2` components, subsequently polished by NMF. - -Arguments: +julia> W_GT, H_GT = generate_ground_truth(); -- `strategy`: see above. Defaults to `GsvdInitialization.truncating`. +julia> X = W_GT * H_GT; +``` -- `X`: non-negative data matrix +the 10 ground-truth components -- `W` and `H`: initial NMF factorization +First, run standard NMF on `X`, initialized with NNDSVD. Two precautions aim for the best possible result from NMF alone: -- `n2`: the number of components in augmented factorization +- `maxiter` is set generously (and we verify convergence below), so the run stops at the convergence tolerance, not at the iteration limit; +- NNDSVD is seeded with the full `svd`, which gives higher-quality results than a randomized SVD. -- `f`: SVD (or Truncated SVD) of `X` +Despite these precautions, the result leaves much to be desired: -Keyword arguments: +```julia +julia> result_nmf = nnmf(X, 10; init=:nndsvd, alg=:cd, tol=1e-4, maxiter=10^4, initdata=svd(X)); -- `tol_final`: the tolerance of the NMF polishing step, default: 1e-4 +julia> result_nmf.converged # stopped at the tolerance, not the iteration cap +true -Other keyword arguments are passed to `NMF.nnmf`. +julia> sum(abs2, X - result_nmf.W*result_nmf.H) / sum(abs2, X) +0.09999800028665384 +``` ------ +components found by standard NMF -W_augmented, H_augmented, Λ = **gsvdrecover**([strategy,] X, W0, H0, kadd, f) +The factorization is imperfect: two components are identical, and two features share a single component. +Now run GSVD-NMF on `X` (also initialized with NNDSVD) and compute the new reconstruction error: -Augment components for `W0` and `H0` without polishing by NMF. -`strategy` defaults to `GsvdInitialization.truncating`; pass -`GsvdInitialization.joint_nnls` or a user-defined callable for alternative -augmentation paths. +```julia +julia> result_gsvd, Λ = gsvdnmf(X, 9 => 10; alg=:cd, tol_final=1e-4, tol_intermediate=1e-2, maxiter=10^4); -Outputs: +julia> W_gsvd, H_gsvd = result_gsvd.W, result_gsvd.H; -`W_augmented`, `H_augmented`: the full augmented NMF factors (with `kadd` extra -components appended to `W0`/`H0`) +julia> sum(abs2, X - W_gsvd*H_gsvd) / sum(abs2, X) +1.2302340443302435e-10 +``` -`Λ`: generalized singular values used to rank the candidate augmentation directions +An imperfect 9-component factorization was augmented by `gsvdnmf` to an essentially perfect 10-component one. +`Λ` holds the generalized singular values that ranked the candidate augmentation directions, a useful diagnostic of which directions the algorithm chose. +Here are the new components: -Arguments: +components found by GSVD-NMF -`X`: non-negative 2D data matrix +## API overview -`W0`: NMF solution +Complete signatures and doctested examples are in the REPL help (e.g., type `?gsvdnmf`). -`H0`: NMF solution +- `gsvdnmf(X, n1 => n2; ...)` runs the full pipeline: an initial NMF with `n1` components, augmentation to `n2` components (`n1 < n2 ≤ 2n1`), and a final NMF polish. `gsvdnmf(X, n)` is shorthand for `gsvdnmf(X, n-1 => n)`. +- `gsvdnmf(X, W, H, f; n2, ...)` augments an existing factorization `X ≈ W*H` to `n2` components, using a precomputed SVD `f` of `X`, then polishes with NMF. +- `gsvdrecover(X, W0, H0, kadd, f)` performs the augmentation step alone, adding `kadd` components without the NMF polish. -`kadd`: number of new components +Each function optionally takes a leading `strategy` argument controlling how the augmented factors are assembled: `GsvdInitialization.truncating` (the default), `GsvdInitialization.joint_nnls`, or a user-supplied callable `(X, W0, H0, Hadd) -> (W_augmented, H_augmented)`. -`f`: SVD (or Truncated SVD) of `X` +## Citation ------ +If you use this package, please cite the paper: -## Citation +> Youdong Guo and Timothy E. Holy, "Recovering missing features in nonnegative matrix factorization via generalized singular value decomposition," *iScience* 29(3):114708 (2026). https://doi.org/10.1016/j.isci.2026.114708 -Thanks for citing this work! See the "Cite this repository" link in the "About" bar for format options. +GitHub's "Cite this repository" link (in the About sidebar) provides this in BibTeX and APA formats.