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/Project.toml b/Project.toml
index a8d4032..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"
@@ -20,6 +21,7 @@ Kronecker = "0.5"
LinearAlgebra = "1"
NMF = "1"
NonNegLeastSquares = "0.4"
+OffsetArrays = "1"
SparseArrays = "1"
StableRNGs = "1"
TSVD = "0.4"
@@ -28,12 +30,14 @@ 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"
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", "Documenter", "ExplicitImports", "NMF", "Test", "FileIO", "JLD2", "OffsetArrays", "StableRNGs"]
diff --git a/README.md b/README.md
index a8703e6..dc006a4 100644
--- a/README.md
+++ b/README.md
@@ -3,158 +3,93 @@
[](https://github.com/HolyLab/GsvdInitialization.jl/actions/workflows/CI.yml)
[](https://codecov.io/gh/HolyLab/GsvdInitialization.jl)
[](https://github.com/JuliaTesting/Aqua.jl)
+[](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
```
-
+## 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
-
-
-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:
-
-
-
-
----------------------------
-
-## 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
+
-- `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
+```
------
+
-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:
+
-`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.
diff --git a/src/GsvdInitialization.jl b/src/GsvdInitialization.jl
index 8ab33fe..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,45 +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 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.
-Augment `W` and `H` to have `n2` components, subsequently polished by NMF.
+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`).
-Keyword arguments:
+Remaining keyword arguments are forwarded to `NMF.nnmf`.
-- `tol_final`: the tolerance of the NMF polishing step, default: 1e-4
+# Examples
-- `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.
+```jldoctest
+julia> using LinearAlgebra: svd
-- `truncmult`: the flooring level applied to the augmented factors when
- `alg == :multmse`, default: 1e-5
+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
-Other keyword arguments are passed to `NMF.nnmf`.
+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
-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> result, Λ = gsvdnmf(X, W0, H0, svd(X); n2 = 3);
+
+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),
@@ -62,6 +89,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"))
@@ -82,33 +110,54 @@ 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 data matrix `X`.
+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.
-Arguments:
+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.
+
+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)
n1, n2 = ncomponents
f = tsvd(X, n2)
W0, H0 = nndsvd(X, n1; initdata = (U = f[1], S = f[2], V = f[3]))
@@ -126,39 +175,69 @@ 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.
+Augment the factorization `X ≈ W0*H0` with `kadd` additional components. This
+is the augmentation step of [`gsvdnmf`](@ref), without the final NMF polish.
-`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.
+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.
-Outputs:
+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.
-`W_augmented`, `H_augmented`: the full augmented NMF factors (with `kadd` extra
-components appended to `W0`/`H0`)
+# Arguments
-`Λ`: generalized singular values used to rank the candidate augmentation directions
+- `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.
-Arguments:
+All array arguments, including the factors of `f`, must use 1-based indexing.
-`X`: non-negative 2D data matrix
+# Examples
-`W0`: NMF solution
+```jldoctest
+julia> using LinearAlgebra: svd
-`H0`: 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
-`kadd`: number of new components
+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
-`f`: SVD (or Truncated SVD) of `X`
+julia> W, H, Λ = gsvdrecover(X, W0, H0, 1, svd(X));
+
+julia> size(W), size(H)
+((4, 3), (3, 5))
+
+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
+ # 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)
@@ -170,13 +249,21 @@ 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)
+ 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'))
@@ -189,17 +276,23 @@ 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)
+ 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..1db5900 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,11 @@
using GsvdInitialization
using Test
using Aqua
+using Documenter
using ExplicitImports
using LinearAlgebra, NMF, FileIO
+using OffsetArrays
using StableRNGs
@testset "Aqua" begin
@@ -18,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
@@ -236,6 +243,84 @@ 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
+ # 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