Sparse log-concave nonnegative matrix factorization

Recover a parsimonious set of compounds from multichannel signals, with nonnegative spectra, hard log-concave elution profiles and a smooth background. Fast spectral fusion and local native updates accommodate broad, asymmetric and shark-fin peaks while retaining all original outcome channels.

Usage

splogconNMF(
  X,
  y,
  family = poisson("identity"),
  background = NULL,
  lambda = "gic",
  dispersion = NULL,
  theta = NULL,
  precluster = NULL,
  weights = NULL,
  offset = NULL,
  correlation = c("auto", "none", "lag1"),
  ncores = NULL,
  control = list(),
  start = NULL,
  verbose = FALSE
)

# S3 method for class 'splogconNMF'
fitted(object, ...)

# S3 method for class 'splogconNMF'
nobs(object, ...)

# S3 method for class 'splogconNMF'
coef(object, ...)

# S3 method for class 'splogconNMF'
logLik(object, ...)

# S3 method for class 'splogconNMF'
AIC(object, ..., k = 2)

# S3 method for class 'splogconNMF'
BIC(object, ...)

# S3 method for class 'splogconNMF'
GIC(object, ...)

# S3 method for class 'splogconNMF'
EBIC(object, ...)

# S3 method for class 'splogconNMF'
MBIC(object, ...)

# S3 method for class 'splogconNMF'
print(x, ...)

Arguments

X

Nonnegative scans-by-candidates matrix, typically shifted copies of a point-spread function from convmatrix() or L0glm_blind_deco(). Columns must have contiguous log-concave positive support. Dense and sparse inputs are accepted and stored sparsely; columns are normalized to unit sum.

y

Numeric scans-by-channels matrix or a single-channel vector.

family

Identity-link noise family: gaussian, poisson, quasipoisson, mgcv::nb, Gamma, nb1(), Gamma_scale() or Gamma_censored(). The corresponding names, including "nb", "nb1" and "Gamma_scale", are also accepted. Count likelihoods require integer counts; quasipoisson accepts nonnegative intensities. Gaussian observations may be negative; their fitted latent signals remain nonnegative.

background

Optional nonnegative scans-by-basis matrix, for example splines::bs(time, df=6, intercept=TRUE). Coefficients are nonnegative and unpenalized. NULL fits without a background.

lambda

"gic" (default), "aic", "bic", "ebic", "mbic", or a nonnegative numeric compound penalty. See Details for its scale.

dispersion

Known common Gaussian variance, quasipoisson scale, Gamma dispersion, NB1 variance-to-mean ratio or common Gamma scale. NULL estimates a pooled value from temporal differences and freezes it. Parameters supplied in a signal_families object are used when this argument is omitted. Poisson and NB2 use unit dispersion.

theta

NB2 shape, shared by all original channels. NULL estimates and freezes a temporal moment estimate; a positive value fixes it; a negative value initializes at its absolute value and enables pooled likelihood updates. A fixed shape in the supplied NB family is respected.

precluster

NULL or FALSE uses every channel throughout. TRUE uses up to 32 groups; an integer supplies another maximum. Size-adjusted Poisson clustering, aggregation, selection and fusion operate in that reduced space, followed by refinement on all original channels.

weights

Nonnegative likelihood weights: a scalar, row vector or matrix matching y. NULL uses one. Zero weights omit observations.

offset

Known nonnegative mean contribution, in the same supported formats as weights. NULL uses zero. It is included on the response scale because the link is identity.

correlation

"auto" or "lag1" estimates an optional lag-one working adjustment to the selection penalty. "none" uses independent observations; a number between zero and 0.49 supplies the lag-one correlation.

ncores

OpenMP workers. NULL uses the detected physical cores.

control

Named list: broad_repeats=1L, full_sweeps=1L, final_polishes=1L, native_sweeps=6L, block_sweeps=6L, psf_width=NULL, screen_factor=1, censor=0.5, linkage="ward", joint_trials=NULL, profile_cycles=10L, fusion_passes=5L, profile_fusion=TRUE. joint_trials uses direct Poisson fusion for Poisson data and joint local refinement for the other families. psf_width supplies the typical FWHM in scans. When omitted, local updates use the middle design column and broad-profile fusion uses the median of up to 21 representative design columns. screen_factor scales the initial selection penalty. censor is the Gamma detection limit: recorded zero intensities are left-censored below this value; positive values must be at least the limit. Zero requests an uncensored Gamma model. A supplied Gamma family constructor provides its own detection limit unless overridden here. linkage may be "ward" or "bregman". joint_trials allows overlapping compounds to adjust together during local fusion refinement. profile_cycles bounds the original-channel shape/background refinement; zero disables it. fusion_passes bounds broad-anchor and compatible-tail fusion after that refinement. profile_fusion enables a final conservative fusion pass on all original channels, followed by one local likelihood-based profile/spectrum refit.

start

Optional list H, S and alpha, with original dimensions. H has scans in rows and compounds in columns; S has compounds in rows and all original channels in columns. alpha contains background coefficients and defaults to zero. Profiles are normalized with the inverse change applied to spectra. A restart bypasses initial screening and preclustering; nonnegative spectral entries can reopen during refitting.

verbose

Print stage timings, objective changes and compound counts.

object, x

A fitted splogconNMF object.

...

Further method arguments.

k

Complexity multiplier for stats::AIC().

Value

An object of class splogconNMF, containing sparse unit-sum profiles H, spectra S, background coefficients alpha, the background basis, fitted means mu, rank, resolved noise settings spec, original data, offsets and weights, objective, loglik, stage ledger history, ledger_error, aggregation and profile_fusion diagnostics, controls and timings seconds. coef() returns spectra and fitted() the full mean matrix. logLik(), AIC(), BIC(), GIC(), EBIC() and MBIC() score the returned model. Columns of H and the matching rows of S are ordered by increasing apex scan, also returned in apex_scan. Their labels identify the same compounds.

Details

Model and selection criterion. With scans in rows, write $$\mu=HS+B A+O,\qquad H,S,A\geq0.$$ Each profile has unit sum, contiguous positive support and \(\Delta^2\log H_{\cdot r}\leq0\) on that support. The objective is $$J=-2\sum_{ij}w_{ij}\ell(y_{ij};\mu_{ij},\nu)+\lambda_R K,$$ where \(K\) is the number of compounds and \(\nu\) contains the shared noise parameters. An IC preset uses \(\lambda_R=m\kappa\lambda_{IC}\), with \(m\) original channels. A numeric lambda supplies the total compound penalty before multiplication by \(\kappa\). The initial grouped fit uses \(\lambda_R/m\). Each compound contributes one selection count; spectral entries and profile complexity have no separate penalty. This is a compound-count model-selection convention. The IC methods use \(mK\) plus the number of estimated shared noise parameters as their complexity count; it is distinct from counting every free profile ordinate.

Initialization and local updates. A variance-matched grouped Poisson working fit supplies candidate components, using L0glm.fit() with relaxed refitting and compulsory background columns. The subsequent updates use the specified family's exact score and likelihood. For a spectral coefficient or profile ordinate, the conditional mean is \(\mu_i=r_i+x_i b\), with the other components fixed in \(r_i\). If \(g_i\) and \(I_i>0\) are the negative log-likelihood score and a positive working curvature, a nonnegative quadratic proposal is $$b^+=\max\left\{0,b-\frac{\sum_i w_i x_i g_i} {\sum_i w_i x_i^2 I_i}\right\}.$$ Backtracking checks the actual loss, with an explicit zero-boundary check. Poisson and quasipoisson use a specialized bracketed scalar likelihood solve. Background terms use small constrained quadratic solves, and all channels run independently through OpenMP. Poisson profile proposals use four warm native IRLS iterations. During a full profile sweep, an adaptive ridge tether centred on the preceding profile stabilizes weak tails, with local penalty weights proportional to the inverse squared profile ordinate. The final acceptance test uses the complete likelihood-and-rank criterion.

Spectral clustering proposes fusions without fitting a model. It compares normalized spectra through Jensen–Shannon information loss and joins overlapping compatible profiles. A least concave majorant of the log profile enforces hard log concavity. Local refinement reallocates signal between overlapping compounds and checks the actual criterion. Every complete stage is accepted against its full working data, and final refinement and scoring use all original observations. history records each accepted or rejected stage; ledger_error checks the accumulated criterion change.

Broad-profile fusion. After refitting the original channels, an additional algebraic pass compares elution profiles. Candidate pairs must have cosine similarity at least 0.95 and at least one profile with FWHM at least four times the typical PSF FWHM. This can reunite complementary ion sets assigned to nearly proportional broad profiles. The area-weighted profile is projected to hard log concavity and the spectra are added. Each merge must improve the original-channel likelihood-and-compound criterion with the noise scale, weights, offsets and background fixed. Candidate fusion requires no refits. One final local likelihood pass updates the profiles and spectra in affected elution regions. It uses Fisher scoring in log-profile coordinates with piecewise-linear interpolation and linear concavity constraints; the noise likelihood is evaluated on the original intensity scale. Backtracking preserves likelihood improvement. Neighboring profile modes stay within one quarter of their starting FWHM, while merged profiles can move freely. In log-profile coordinates \(\eta=\log H_{\cdot r}\), the score and working information are \(g_t=H_{tr}\sum_c w_{tc}S_{rc}g_{\mu,tc}\) and \(I_t=H_{tr}^2\sum_c w_{tc}S_{rc}^2I_{\mu,tc}\). With interpolation matrix \(C\), the quadratic step minimizes \(\sum_t I_t\{(Cb)_t-\eta_t+g_t/I_t\}^2\), subject to nonincreasing slopes of \(Cb\). Zero-information rows receive numerical stabilization. profile_fusion records the merge audit and source component groups, numbered in temporal order as they entered the final fusion pass. Its source_components list follows the returned temporal ordering.

Noise models. Gaussian variance is \(\phi\); quasipoisson variance is \(\phi\mu\); NB2 variance is \(\mu+\mu^2/\theta\); standard Gamma variance is \(\phi\mu^2\). NB1 has variance \(\phi\mu\), while common-scale Gamma has shape \(\mu/b\) and scale \(b\), hence variance \(b\mu\). Their exact scores and densities are given in signal_families. A censored Gamma zero contributes a log CDF at the detection limit, and a positive observation contributes its log density. For quasipoisson, \(\ell_Q=\{y\log\mu-\mu-\log\Gamma(y+1)\}/\phi\); logLik() marks this scaled kernel with quasi=TRUE, and the regular IC method names return quasi-information criteria. Compare candidate ranks under the same family and frozen scale.

For an updated NB2 shape, the pooled score is $$U_\theta=\sum_{ij}w_{ij}\left\{\psi(y_{ij}+\theta)-\psi(\theta) +\log\theta+1-\log(\theta+\mu_{ij}) -\frac{y_{ij}+\theta}{\theta+\mu_{ij}}\right\}.$$ Bounded Newton proposals in \(\log\theta\) are backtracked against the pooled likelihood. One common shape is retained across original channels.

Temporal calibration and aggregation. temporal_dispersion() uses fourth and eighth differences in quieter parts of a smoothed pilot signal. This estimates the marginal noise scale; an optional lag-one working correction multiplies the rank penalty by \(\kappa=1+2\widehat\rho_1\). It leaves the fitted variance parameter on its marginal scale. This local covariance approximation is recorded in spec; use correlation="none" for independent simulated data.

Preclustering uses glm_cluster() with size_fac=colSums(pmax(y,0)), removing differences in total ion intensity when comparing channel shapes. Poisson and common-parameter NB1 latent signals add exactly; independent common-scale Gamma latent signals also add. Standard NB2 and Gamma use moment-matched aggregate parameters. Censored Gamma sums and differing channel weights give working approximations in the reduced space. The returned model always contains all original channels, their offsets and their original-channel likelihood. The compound penalty uses the original channel count throughout aggregation. Selection and initial fusion run in the reduced space. Lifting the spectra to the original channels is followed by a spectral/background refit and local shape refinement, with no repetition of the reduced-space search.

Smooth local profiles. Original-channel refinement fits \(\eta(t)=a+bt-\tfrac12\sum_j c_j\{(t-u_j)_+^2-(t-v_j)_+^2\}\), with \(c_j\geq0\). Thus \(\eta''(t)\leq0\) and \(H(t)=\exp\{\eta(t)\}\) is log-concave. Fisher scoring uses the original observation model and the chain-rule score/information given above. Small constrained quadratic solves and likelihood backtracking update each local profile. Joint spectral updates separate overlapping compounds; compatible broad-profile and tail fusions are accepted under the same original-channel likelihood and compound-count penalty.

See also

GCMSdeco(), L0glm_blind_deco(), L0glm.fit(), glm_cluster()

Examples

set.seed(12)
n <- 300L
psf <- convkernel("gaussian", 21L, list(sigma = 2))
X <- convmatrix(psf, n)
H <- cbind(X[, c(45, 95, 160, 235)],
           pmax(0, 1 - abs(seq_len(n) - 185) / 22))
H <- H %*% Matrix::Diagonal(x = 1 / Matrix::colSums(H))
S <- matrix(rexp(5 * 8, rate = 1 / 3000), 5, 8)
B <- splines::bs(seq_len(n), df = 4, intercept = TRUE)
mu <- as.matrix(H %*% S) + B %*% matrix(20, 4, 8)
Y <- matrix(rpois(length(mu), mu), n, 8)
fit <- splogconNMF(X, Y, background = B, family = poisson("identity"),
  precluster = 4L, correlation = "none", ncores = 2L)
fit
c(compounds = fit$rank, GIC = GIC(fit), seconds = fit$seconds[["total"]])
# Broad profiles with nearly identical shapes can share complementary ion sets.
# A conservative fusion pass runs on all original channels after preclustering.
fit$profile_fusion$merges
matplot(fit$H, type = "l", lty = 1, xlab = "Scan", ylab = "Elution profile")

# For overdispersed intensities, estimate one common noise scale from temporal
# differences. Selection and fusion use four groups; the result has all channels.
Yq <- 5 * matrix(rpois(length(mu), mu / 5), n, 8)
fitq <- splogconNMF(X, Yq, background = B, family = quasipoisson("identity"),
  precluster = 4L, correlation = "none", ncores = 2L)
GIC(fitq)

# Other identity-link choices include gaussian(), mgcv::nb(), Gamma(),
# nb1(dispersion=5, link="identity"), and
# Gamma_scale(scale=5, link="identity", censor=0.5).