L0-Penalised Generalised Best-Rank NMF

Factorises a non-negative data matrix \(D\) into \(W\,H\) under any exponential-family noise model while **simultaneously selecting the best rank** via hard (L0 or grouped L0) best-subset penalties. This NMF can for example be used to separate (deconvolute) overlapping peaks in GC/MS (gas chromatography - mass spectrometry) datasets, allowing one to simultaneously estimate compound elution profiles (W), electron impact mass spectra (H) and the true number of compounds that are present in the chromatogram (the best rank).

Usage

L0BestRankGNMF(
  D,
  k,
  Wstart = NULL,
  grouped_Hupdate = FALSE,
  IC = "gic",
  max_iter = 5,
  tol = 1e-04,
  method = 1L,
  parallelize = FALSE,
  ncores = if (parallelize) {
     round(parallel::detectCores()/2)
 } else {
     1
 },
  family = gaussian(link = "identity")
)

Arguments

D

Numeric (or Matrix) \(n\times p\) input matrix with non-negative entries.

k

Initial rank of the non-negative double SVD (NNDSVD) that is used for initialisation if Wstart=NULL. Allowed to over-estimate the true rank.

Wstart

Optional \(n\times k_0\) non-negative matrix with starting values for \(W\). If supplied, k is ignored and set to ncol(Wstart).

grouped_Hupdate

Logical; if TRUE a **grouped** best-subset penalty is used in the \(H\) update (one latent component on/off for all outcome channels, i.e. best subset of group selection).

IC

Information criterion for the penalty preset; see information_criteria.

max_iter

Maximum number of outer iterations.

tol

Relative change in penalised objective to determine convergence.

method

Eigen solver to use, as in `L0glm.fit()` (e.g. 1L=Cholesky LLT, 4L=least squares conjugate gradient).

parallelize

Use parallelisation? May only be faster for larger problems.

ncores

Number of cores to use when "parallelize=TRUE".

family

A family object defining the exponential family and link; identity-link Poisson provides Lee-Seung KL-NMF but with additional L0 or grouped L0 regularisation.

Value

A list with components

W

Estimated \(n\times \hat k\) W matrix.

H

Estimated \(\hat k\times p\) H matrix.

dev

Final penalised objective value, i.e. scaled Bregman (deviance) loss (=-2*log-likelihood), with added elementwise \(L_0\) or grouped row-\(L_0\)/\(L_{2,0}\) penalties).

k

Discovered rank \(\hat k\).

iter

Number of outer iterations run.

dev_trace

Vector of penalised objective values over iterations.

lambdaW

Lambda penalty used for W.

lambdaH

Lambda penalty used for H.

Details

The routine alternates two calls to `L0glm.fit()` and works via cyclical block coordinate descent:

H-update

multivariate GLM with L0 or optional **grouped** L0 penalty (selects or discards an entire latent component across all outcome channels).

W-update

plain L0 GLM row-wise, warm-started from the previous iterate.

After each full sweep factors that are shrunk exactly to 0 are pruned, which means one can start with k larger than the true rank. The optimal best rank is then typically very well recovered.

The objective minimised after each full W-H sweep is $$ \underbrace{\frac{1}{\phi}\sum_{ij} D_{F}\!\bigl(D_{ij}\,\|\,\mu_{ij}\bigr)} _{-2\,\mathrm{logLik}} \;+\; \lambda_W \sum_{k}\mathbf 1[\|W_{\cdot k}\|_2>0] \;+\; \begin{cases} \lambda_H\sum_{k}\mathbf 1[\|H_{k\cdot}\|_2>0], &\text{if } grouped\_Hupdate=TRUE\\[6pt] \lambda_H\!\sum_{k,j}\mathbf 1[H_{kj}\neq 0], &\text{otherwise.} \end{cases} $$

The grouped penalty is a group-\(L_0\)/row-\(L_0\) penalty, equivalently an \(L_{2,0}\) penalty, counting whole latent components as on/off. The ungrouped penalty is the plain elementwise \(L_0\) penalty used by L0glm.fit.

Implementation notes

* Dispersion \(\phi\) is re-estimated after every sweep for Gaussian, Gamma and inverse-Gaussian models (or fixed to 1 for Poisson). * Rows of \(H\) are rescaled to unit \(L^2\) norm at the end; corresponding columns of \(W\) are scaled back to keep \(W\,H\) unchanged.

Examples

## Small nonnegative Gaussian matrix factorization example
if (requireNamespace("Matrix", quietly = TRUE)) {
  library(Matrix)
}
set.seed(123)
Wtrue <- matrix(rexp(30 * 3), 30, 3)
Htrue <- matrix(rexp(3 * 20), 3, 20)
D <- pmax(Wtrue %*% Htrue + matrix(rnorm(30 * 20, sd = 0.05), 30, 20), 0)
Wstart <- Matrix(Wtrue + matrix(runif(30 * 3, 0, 0.1), 30, 3), sparse = TRUE)

nmf <- L0BestRankGNMF(D,
                      Wstart   = Wstart,
                      IC       = "bic",
                      max_iter = 1,
                      family   = stats::gaussian("identity"))

cat("Rank recovered by L0BestRankGNMF:", nmf$k, "\n")

if (FALSE) { # \dontrun{
## ---------- Larger GC/MS deconvolution demo with Poisson noise ------------
if (requireNamespace("Matrix", quietly = TRUE)) {
  library(Matrix)
}
n  <- 200                 # time points
k0 <- 50                  # true rank/ number of compounds
set.seed(123)

peak_width <- 2.5         # unknown, but can be estimated from standards or tuned from chromatograms
psf_fun <- function(x,u,w = peak_width) exp(-0.5*((x-u)/w)^2)

sim <- simulate_spike_train(n       = n,
                            k       = k0,
                            family  = "poisson",
                            multichannel = TRUE,
                            mean_beta    = 1e4,
                            sparse       = TRUE,
                            psf          = psf_fun)

D <- sim$y                # observed GC/MS chromatogram

## build a design matrix with time-shifted copies of the PSF
Wband <- sim$X            # n x n banded matrix with time-shifted copies of PSF
## pre-select putative peaks using an unregularised nonnegative
## identity-link Poisson fit on the total ion chromatogram
sel   <- L0glm.fit(Wband, rowSums(D),
                   lambda = 0, nonnegative = TRUE,
                   family = poisson(link=identity),
                   solver = "osqp")$coefficients != 0
Winit <- as.matrix(Wband[, which(sel), drop = FALSE])

## sparse Poisson Best-Rank NMF with grouped-L0 penalty
nmf <- L0BestRankGNMF(D,
                      Wstart          = Matrix(Winit, sparse = TRUE),
                      grouped_Hupdate = TRUE,
                      IC              = "bic",
                      max_iter        = 5,
                      family          = poisson("identity"))

cat("True rank:", k0, "\n",
    "Rank recovered by Poisson L0BestRankGNMF:", nmf$k, "\n")
matplot(nmf$W, type="l", xlab="Time", ylab="Signal",
main="Elution profiles inferred by Poisson L0BestRankGNMF")

## ---------- Comparison with MCR-ALS (alternating non-negative least squares) -
if (requireNamespace("ALS", quietly = TRUE)) {
  library(ALS)
  als_fit <- als(CList   = list(Winit),
                 S       = matrix(1, ncol(D), ncol(Winit)),
                 PsiList = list(as.matrix(D)),
                 uniC    = TRUE,  # impose unimodality constraints
                 maxiter = 8,
                 normS   = 1)

  cat(" Rank recovered by MCR-ALS:", nrow(t(als_fit$S)), "\n")
}
} # }