*! version 1.0 16sep2022

Title

nmf -- matrix decomposition using non-negative matrix factorization (NMF).

Syntax

Perform matrix decomposition using NMF:

nmf
[varlist]
, k(#) [options]

Main
* k(#) specifes the rank, k, of the factorisation.

epoch(#) maximum number of iterations over which to minimise the error function.

initial(string) declares the matrix initialisation method.

loss(string) declares the loss function used to calculate ||A - WH||.

stop(#) declares the early stopping delta threshold for convergence.

nograph suppress graph of epoch vs. loss function.

* k(#) is required.

by is not allowed. fweights are not allowed.

Description

The aim of nmf is to decompose matrix A into a lower-rank matrix approximation such that
X ≈ WH, where A, W and H are made up of non-negative, real numbers. Missing values are
permitted. If X is an n x m matrix, the dimensions of W and H will be n x k and k x m,
respectively, whereby k represents the rank of the decomposition. In many cases, a good
approximatation for A may be achieved with k << rank(A).

NMF is a NP-hard problem. As such, multiplicative updates of matrices W and H are iteratively
performed in order to minimise the generalized error function, ||A - WH||. This
implements the methods first reported by Paatero and Tapper[1] and later popularised by Lee and
Seung[2, 3].

Running NMF results in the generation of three new frames, W, H and error that
store the basis and coefficient matrices and a summary of the error over each epoch, respectively.
These can be accessed using: frame change W, frame change H and frame change error.

Options

Main

k(#) is required and specifies the rank for the decomposition. The value of k must be greater
than or equal to 2, and must be less than the minimum of the number of columns and rows of the initial
matrix to be factorised.

Options

epoch(#) sets the maximum number of epochs (iterations) over which the decomposition is optimized.
If this is not specified, NMF runs for 200 epochs by default. If convergence has not been reached by the
number of epochs specified in epoch() (or by 200 epochs, if no value is specified), an error message will
advise that a greater number of epochs should be used.

initial(option)
indicates how matrices W and H are initialized. The available options are:

randomu,
the default, specifies that W and H are initialised by sampling values from a uniform distribution, in the range [min(A), max(A)].

randomn,
specifies that W and H are initialised by sampling values from a normal distribution, in the range [1, 2].

nndsvd,
specifies that W and H are initialised using non-negative double singular value decomposition (NNDSVD)[4]. This is deterministic (non-random) and was designed
to ehance the initialisation stage of NMF.

nndsvda,
specifies that W and H are initialised using NNDSVD, in which zeroes are filled with the average of A. It is suggested that this may be better when sparsity is not desired.

nndsvdar,
specifies that W and H are initialised using NNDSVD, in which zeroes are filled with small random values.

loss(option)
indicates the loss function that will be minimised during the matrix decomposition. This is the error function for ||A - WH||.
Options include: eu, is and kl. Note that the Itakura-Saito divergence (is) requires no missing data to be present in A.

eu,
the default, specifies the Frobenius (Euclidean) distance. This may be preferred when matrix A contains continous data.

kl,
specifies specifies the generalized Kullback-Leibler divergence. This may be preferred when matrix A contains binary or count data.

is,
specifies the Itakura-Saito divergence.

stop(#)
sets the early stopping threshold for convergence; if stop(0) is set, optimisation will continue for the set number of epochs without early stopping.
If ((previous error - current error) / error at initiation) < stop tolerance, convergence has occured and nmf terminates.
The default value is 1.0e-4.

nograph
suppresses plotting of a line graph to depict the loss function decreasing with each successive epoch.

Remarks

A number of demonstration .do files are available, complete with sample data. These are metagenes.do, imputation.do, faces.do and trajectories.do and can be found at: https://github.com/jonathanbatty/stata-nmf/tree/main/examples

References

[1] Paatero, P. and Tapper U (1994). Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values. Environmetrics, 5: pp. 111-126.

[2] Lee, D. and Seung, H (1999). Learning the parts of objects by non-negative matrix factorization. Nature, 401, pp. 788–791.

[3] Lee, D. and Seung, H (2000). Algorithms for Non-negative Matrix Factorization. Advances in Neural Information Processing Systems 13: Proceedings of the 2000 Conference, MIT Press. pp. 556–562.

[4] Boutsidis, C. and Gallopoulos, E (2008). SVD-based initialization: A head start for nonnegative matrix factorization. Pattern Recognition, 41(4): pp. 1350-1362, .

Examples

. nmf p*, k(5) epoch(100)

. nmf k*,	k(15) epoch(100) initial(randomu) stop(1.0e-4) loss(kl) nograph

Acknowledgements

This project was funded by the AI for Science and Government Fund, via the Alan Turing Institute Health Programme (TU/ASG/R-SPEH-114).
Jonathan Batty received funding from the Wellcome Trust 4ward North Clinical Research Training Fellowship (227498/Z/23/Z; R127002).

Suggested citation

Batty, J. A. (2024). Stata package ``nmf'': an implementation of non-negative matrix factorisation in Stata (Version 1.0) [Computer software]. https://github.com/jonathanbatty/stata-nmf