API reference

NMFMerge

Improve nonnegative matrix factorization (NMF) by over-factorizing and then merging components in pairs. Factorizing X ≈ W * H with more components than ultimately desired and merging the closest pairs tends to reach better and more reproducible optima than factorizing directly into the target number of components (Guo & Holy, IEEE Trans. Signal Process. 2025, doi:10.1109/TSP.2025.3585893).

The high-level entry point is nmfmerge, which runs the whole pipeline. The building blocks are also exported: colnormalize to unit-normalize the columns of W, merge_pq to merge a factorization down to a chosen number of components while recording the merge sequence, and merge_replay to replay a prefix of that sequence.

result = nmfmerge([queuepenalty], X, ncomponents; tol_final=1e-4, tol_intermediate=sqrt(tol_final), W0=nothing, H0=nothing, kwargs...)

Perform "NMF-Merge" on data matrix X: factorize X with more components than ultimately desired and then iteratively merge components in pairs, which tends to reach better and more reproducible optima than a direct factorization (Guo & Holy, IEEE Trans. Signal Process. 2025, doi:10.1109/TSP.2025.3585893).

Return an NMF.Result whose W and H fields hold the final factorization X ≈ W * H with the requested final number of components.

Arguments:

• queuepenalty: a function of the form f(E, h1sq, h2sq) that computes the penalty for merging two components, where E is the merge error (the increase in reconstruction error) and h1sq, h2sq are the squared norms of the corresponding rows of H. Default: ssdpenalty (f(E, h1sq, h2sq) = E), which merges purely in order of increasing reconstruction error.

• X::AbstractMatrix: the data matrix to be factorized.

• ncomponents::Pair{<:Integer,<:Integer}: in the form n1 => n2, merging from n1 components to n2, where n1 is the number of components for the overcomplete NMF and n2 is the number of components for the final NMF. We require n1 >= n2. Alternatively, ncomponents can be an integer denoting the final number of components; in that case nmfmerge defaults to an approximate 20% component excess before merging.

Keyword arguments:

• tol_final: the tolerance of the final NMF.

• tol_intermediate: the tolerance of the initial and overcomplete NMF.

• W0, H0: initialization for the initial NMF. If at least one of W0 and H0 is nothing, NNDSVD is used for initialization.

Other keyword arguments are passed to NMF.nnmf.

Examples

julia> X = abs.(randn(100, 50));        # nonnegative data

julia> result = nmfmerge(X, 5);         # merge down to 5 components

julia> size(result.W), size(result.H)
((100, 5), (5, 50))

julia> result = nmfmerge(X, 8 => 5);    # 8 overcomplete components, merge to 5

Wnormalized, Hnormalized = colnormalize(W, H, p=2)

Normalize the factorization so that each retained column satisfies norm(W[:, i], p) ≈ 1, rescaling the matching row of H by the same factor so that W * H is unchanged. p is the norm order passed to LinearAlgebra.norm, so any real order is accepted (e.g. 1, 2, Inf).

Columns of W with zero norm carry no information and are dropped together with their H rows, so the returned factors may have fewer components than the input.

merge_pq and nmfmerge require unit 2-norm columns, so use the default p=2 when preparing input for the merge. Other orders are available for unrelated normalization needs.

The component axis (columns of W, rows of H) must be one-based; the feature axis (rows of W) and sample axis (columns of H) may have any axes.

Examples

julia> W = [3.0 0.0; 4.0 0.0];   # second column has zero norm

julia> H = [1.0 2.0; 3.0 4.0];

julia> Wn, Hn = colnormalize(W, H);

julia> Wn
2×1 Matrix{Float64}:
 0.6
 0.8

julia> Hn
1×2 Matrix{Float64}:
 5.0  10.0

Wmerge, Hmerge, mergeseq = merge_pq([queuepenalty], W::AbstractArray, H::AbstractArray; nstop=1, errstop=typemax(...))

Merge components in W and H (columns in W and rows in H). Merging stops at whichever of the two criteria nstop and errstop is reached first.

Arguments:

-queuepenalty: The same as in nmfmerge. Default: ssdpenalty (f(E, h1sq, h2sq) = E).

• W::AbstractArray: The basis matrix with unit 2-norm columns (e.g. from colnormalize with the default p=2). Columns that are not 2-normalized throw an ArgumentError.

• H::AbstractArray: The coefficient matrix.

Keyword arguments:

• nstop::Integer (default: 1): a floor on the number of components. Merging never produces fewer than nstop components. The default nstop=1 merges to a single component.

• errstop (default: typemax(...)): a ceiling on the per-merge cost. Merging never performs a merge costing more than errstop. The cost is the queuepenalty value; under the default penalty it is the merge error.

The defaults on these keywords allow merging all the way down to a single component.

The component axis (columns of W, rows of H) is a plain enumeration and must be one-based. The feature axis (rows of W) and sample axis (columns of H) may have any axes and are carried through to the output; the merged component axis is one-based.

Outputs:

Wmerge and Hmerge are the merged results; the number of surviving components is nstop unless errstop stopped merging earlier.

mergeseq is the sequence of merges as (id1, id2, err) tuples, in the order they were performed. id1 and id2 are the ids of the merged components and err is the reconstruction error (the merge penalty) incurred by that merge. Ids larger than the number of columns in W refer to components produced by earlier merges. The accumulating err values let a caller merge all the way down (nstop == 1) and locate a "knee" at which to stop, then replay the corresponding prefix of mergeseq with merge_replay.

Examples

julia> W = [1.0 0.0 0.0; 0.0 1.0 1.0];   # columns 2 and 3 are identical

julia> H = [2.0 5.0; 3.0 7.0; 4.0 11.0];

julia> Wmerge, Hmerge, mergeseq = merge_pq(W, H; nstop=2);

julia> Hmerge
2×2 Matrix{Float64}:
 2.0   5.0
 7.0  18.0

julia> mergeseq
1-element Vector{Tuple{Int64, Int64, Float64}}:
 (2, 3, 0.0)

Wmerge, Hmerge = merge_replay(W, H, mergeseq)

Merge components in W and H (columns in W and rows in H) by replaying the sequence of merge-pair ids in mergeseq, returning the merged factors.

Each entry of mergeseq supplies the pair of ids (id1, id2) to merge; any further fields are ignored, so the (id1, id2, err) triples returned by merge_pq can be replayed directly. Replaying a prefix of a merge_pq schedule reproduces the factors merge_pq would have returned had it stopped at the corresponding number of components.

As in merge_pq, the component axis must be one-based while the feature and sample axes are carried through to the output.

Examples

julia> W = [1.0 0.0 0.0; 0.0 1.0 1.0];

julia> H = [2.0 5.0; 3.0 7.0; 4.0 11.0];

julia> merge_replay(W, H, [(2, 3)])
([1.0 0.0; 0.0 1.0], [2.0 5.0; 7.0 18.0])

ssdpenalty(E, h1sq, h2sq)

The default merge penalty: the merge error E itself, ignoring the squared norms h1sq, h2sq of the two H rows. With this penalty merge_pq and nmfmerge merge purely in order of increasing reconstruction error.

Pass a custom f(E, h1sq, h2sq) as the leading argument to those functions to weight merges differently.

Examples

julia> NMFMerge.ssdpenalty(1.5, 2.0, 3.0)
1.5