Sparse log-concave NMF

splogconNMF() extends sparse deconvolution to compounds whose profiles are broader or more asymmetric than the point spread function. It combines an initial grouped L0glm.fit() with fast spectral fusion and local constrained updates. Every compound has a nonnegative spectrum and a log-concave elution profile. A low-dimensional spline basis absorbs the smooth background.

fit <- splogconNMF(X, Y, background = B, family = poisson("identity"),
                  precluster = 32L)
fit$rank                  # number of compounds
fit$H                     # sparse, unit-sum profiles
coef(fit)                 # spectra in every original channel
fitted(fit)               # profiles x spectra + background + offset
GIC(fit)

For experimental GC/MS intensities, family=quasipoisson("identity") allows extra variation through one pooled temporal dispersion estimate. Gaussian, negative-binomial and Gamma likelihoods are also available. The help page ?splogconNMF gives all arguments and examples. The GC/MS deconvolution tutorial shows complete experimental examples, including CDF input, MSP export and all four plot types.

Model and criterion

For scans i, channels j and compounds r,

\mu_{ij}=\sum_{r=1}^K H_{ir}S_{rj}+(BA)_{ij}+O_{ij},\qquad H,S,A\geq0.

The columns of H have unit sum, contiguous support and \Delta^2\log H_{\cdot r}\leq0. Spectra S therefore have the interpretation of integrated channel intensities. The objective is

J=-2\sum_{ij}w_{ij}\ell(y_{ij};\mu_{ij},\nu)+\lambda_R K, \qquad \lambda_R=m\kappa\lambda_{IC}.

Here m is the original channel count, \nu contains pooled noise parameters and \kappa=1 for independent observations. Numeric lambda specifies the total compound penalty before the optional covariance adjustment. Background coefficients and profile complexity have no additional penalty. The original channel count is retained in this penalty when the observations are aggregated, so changing the number of outcome groups does not change the compound penalty. Dividing by the channel count gives the equivalent criterion -2\ell/m+\kappa\lambda_{IC}K, with one rank penalty per compound. The reported IC methods use the compound-count convention mK, plus any estimated shared noise parameters.

Generalized matrix factorization

The broad statistical setting is generalized matrix factorization (GMF): an exponential-family mean is linked to a low-rank product, optionally with observed covariate effects and offsets. See Kidziński et al. (2022) and Castiglione et al. (2026). For example, the latter mean-model structure is

g(\mu)=X_{\rm row}B_{\rm row}^{\top} +B_{\rm col}X_{\rm col}^{\top}+UV^{\top}+O.

Our model uses the identity link, U=H, V^{\top}=S, and X_{\rm row}=B, B_{\rm row}^{\top}=A for the smooth background; the column-covariate term is omitted. It adds nonnegativity, unit-area log-concave profiles and a compound-count penalty to select the rank. Thus sparse log-concave nonnegative GMF describes the model, while the function retains the familiar splogconNMF() name.

Deviance matrix factorization and Bregman NMF

splogconNMF() is a rank-penalized, nonnegative deviance matrix factorization with log-concave profiles. It belongs to the deviance-based framework of Wang & Carvalho (2023), using an identity link for the additive mean \mu=HS+BA+O. The loss also connects directly to NMF with Bregman divergences.

For a regular exponential-dispersion model with fixed noise parameters,

\ell(y;\theta,\phi)=\frac{y\theta-b(\theta)}{\phi}+c(y,\phi), \qquad \mu=b'(\theta),

let \Psi=b^* be the convex conjugate of b. Its Bregman divergence is

D_\Psi(y,\mu)=\Psi(y)-\Psi(\mu)-\Psi'(\mu)(y-\mu).

Since \Psi'(\mu)=\theta(\mu) and \Psi(\mu)=\mu\theta(\mu)-b\{\theta(\mu)\}, substitution gives

\begin{aligned} D_\Psi(y,\mu) &=y\theta(y)-b\{\theta(y)\}-y\theta(\mu)+b\{\theta(\mu)\},\\ 2\{\ell(y;y,\phi)-\ell(y;\mu,\phi)\} &=\frac{2}{\phi}D_\Psi(y,\mu). \end{aligned}

Here the likelihood on the second line is written in mean coordinates. Boundary observations use the corresponding limits. Thus, conditional on the noise parameters, fitting is equivalent to minimizing

2\sum_{ij}\frac{w_{ij}}{\phi}D_\Psi(y_{ij},(HS+BA+O)_{ij}) +\lambda_RK

under the nonnegativity, profile-normalization and log-concavity constraints. The saturated likelihood contributes a constant in these conditional model comparisons. The corresponding divergences include:

Noise model	Generator \Psi(u)	Half unit deviance D_\Psi(y,\mu)
Gaussian	u^2/2	(y-\mu)^2/2
Poisson	u\log u-u	y\log(y/\mu)-y+\mu
Gamma, uncensored y>0	-\log u	y/\mu-\log(y/\mu)-1
NB2, fixed shape \theta	u\log u-(u+\theta)\log(u+\theta)	y\log(y/\mu)-(y+\theta)\log\{(y+\theta)/(\mu+\theta)\}

For NB2, \phi=1 and the shape \theta enters the generator. Gaussian, Poisson and Gamma therefore give squared-error, generalized KL and Itakura–Saito NMF losses. Quasipoisson gives the KL kernel scaled by its dispersion. Censored Gamma and the mean-dependent-shape extensions use their full likelihood kernels described below; estimated noise parameters are compared using the full likelihood as well.

From a chromatogram to NIST spectra

GCMSdeco() connects the factorization to ANDI-MS .cdf files and NIST .msp spectral libraries:

fit <- GCMSdeco("beer.cdf", "beer.msp", alkanes_file="alkanes.cdf",
                time_range=c(6, Inf))
fit$compounds
coef(fit)                 # spectra in all original mass channels
GIC(fit)

Choose time_range to match the acquisition’s solvent delay or the window of interest. The default workflow uses ten overlapping blind-PSF windows, size-adjusted aggregation to at most 32 outcome groups, and a low-df spline background. Alkane peaks provide an optional calibrated PSF start. Final refinement and exported spectra retain every original mass channel.

For instrumental intensities the default quasipoisson model estimates one shared dispersion. Its variance scale also sets the working precision of the total-ion PSF pilot. The final likelihood uses the original channels and the same frozen noise parameters. The reference page details the working-weight formula, automatic spline dimension, retention-time units and exported metadata.

Local quadratic updates

For a coefficient b\geq0, keep the rest of the signal in an offset r_i and write \mu_i=r_i+x_i b. With negative log-likelihood score g_i and positive working curvature I_i, propose

b^+=\max\left\{0,b-\frac{\sum_i w_i x_i g_i} {\sum_i w_i x_i^2 I_i}\right\}.

Actual-loss backtracking and a zero-boundary check complete each scalar solve. The Poisson branch uses a faster bracketed likelihood solve. Its score \sum_iw_ix_i(1-y_i/\mu_i) is monotone, and its nonnegative root is bounded above by \sum_iw_iy_i/\sum_iw_ix_i. Channels are independent conditional on the profiles, allowing deterministic OpenMP parallelism.

Local profile updates use the original observation likelihood. For \eta_i=\log h_i, \mu_{ij}=r_{ij}+h_i s_j, the chain-rule score and working information are

g_{\eta,i}=h_i\sum_jw_{ij}s_jg_{\mu,ij},\qquad I_{\eta,i}=h_i^2\sum_jw_{ij}s_j^2I_{\mu,ij}.

The working target is \eta_i-g_{\eta,i}/I_{\eta,i}, with tiny stabilization where the information vanishes. This quadratic approximation updates the log profile parameters; the observed intensities remain on their original scale and acceptance uses their actual likelihood.

A smooth concave log profile can be written as

\eta(t)=a+bt-\frac12\sum_j c_j \{(t-u_j)_+^2-(t-v_j)_+^2\},\qquad c_j\geq0,

where the intervals [u_j,v_j) partition the local time window. Its second derivative is -\sum_jc_j\mathbf1_{[u_j,v_j)}(t)\leq0, giving hard log-concavity. The intercept and slope are unrestricted, and a small nonnegative quadratic solve estimates the curvature coefficients. For more detailed local adjustments, piecewise-linear log-profile increments preserve the existing fine-scale curvature; a concave-envelope projection removes numerical constraint residuals. Backtracking between feasible log profiles preserves concavity and accepts only likelihood improvements.

Fusion starts with spectral compatibility and profile overlap. Broad-anchor and contained-tail proposals are refined only within the affected windows; other compounds remain in an offset. The full original-channel criterion decides whether each proposed rank change is retained. A final conservative pass can reunite broad profiles with complementary ion sets when their profile cosine is at least 0.95. Neighbouring modes are guarded during its local refit. control=list(profile_cycles=0L) disables the main shape-refinement stage; control=list(profile_fusion=FALSE) disables the final conservative pass.

Noise model	Variance	Negative log-likelihood score in the mean
Gaussian	\phi	(\mu-y)/\phi
Poisson	\mu	1-y/\mu
Quasipoisson	\phi\mu	(1-y/\mu)/\phi
NB2	\mu+\mu^2/\theta	\theta(\mu-y)/\{\mu(\mu+\theta)\}
Gamma	\phi\mu^2	(\mu-y)/(\phi\mu^2)
NB1	\phi\mu	\{\psi(a)-\psi(a+y)+\log\phi\}/(\phi-1), a=\mu/(\phi-1)
Common-scale Gamma	b_0\mu	\{\psi(\mu/b_0)+\log b_0-\log y\}/b_0

The final two models use the exact distributions documented in signal_families. Their observed positive curvatures are differences of trigamma functions for NB1 and \psi_1(\mu/b_0)/b_0^2 for common-scale Gamma. The latter is convex in the mean; the NB1 curvature is nonnegative because trigamma decreases on the positive real line.

Censored Gamma intensities

Recorded zeros may represent intensities below a detection limit d, whose contribution is \log P(Y<d;\mu). Positive observations contribute their ordinary density. For standard Gamma, set a=1/\phi, z=d/(\phi\mu); its negative log-CDF derivative is

g(\mu)=\frac{z^a e^{-z}}{\Gamma(a)P(a,z)\mu}.

For common-scale Gamma, z=d/b_0 stays fixed and a=\mu/b_0 changes. The convergent incomplete-Gamma series

P(a,z)=\frac{e^{-z}z^a}{\Gamma(a+1)} \left\{1+\sum_{k\geq1}\frac{z^k}{(a+1)\cdots(a+k)}\right\}

provides its shape derivative when z is small. The general branch uses the Gamma CDF and a numerical shape derivative. Both branches are checked against the actual censored likelihood. The default detection limit is 0.5 intensity units and can be changed through control$censor.

Fusion, aggregation and rank selection

1. Obtain a variance-matched grouped deconvolution start, retaining the spline background as compulsory unpenalized terms.
2. Compare normalized spectra by Jensen–Shannon information loss. Propose fusions when their elution profiles overlap and their combined profile can be made log-concave. Clustering itself performs no regression fits.
3. Refine each proposal locally, holding all outside contributions fixed. Overlapping spectra and the candidate profile can adjust together.
4. Accept a rank change when -2(\ell_{new}-\ell_{old})+\lambda_R(K_{new}-K_{old})\leq0.
5. Update broad profiles, smooth background coefficients and spectra, followed by final boundary checks and rank pruning.

Algebraic fusion uses a least concave majorant of the log profile, followed by normalization. Local likelihood updates impose hard concavity constraints in their quadratic solves. Every complete accepted stage decreases the actual criterion, within numerical tolerance. history and ledger_error make this property directly checkable.

With precluster=32L, the initial workflow uses at most 32 aggregated outcome groups. Size factors colSums(pmax(Y,0)) account for different total ion intensities. Final refinement returns to every original channel. Independent NB1 variables sharing p=1/\phi are closed under sums, because their probability-generating functions multiply as

\prod_j\left\{\frac{p}{1-(1-p)t}\right\}^{\mu_j/(\phi-1)} =\left\{\frac{p}{1-(1-p)t}\right\}^{\sum_j\mu_j/(\phi-1)}.

The same property holds for common-scale Gamma latent signals, through \prod_j(1+b_0t)^{-\mu_j/b_0}=(1+b_0t)^{-\sum_j\mu_j/b_0}. Standard NB2 and Gamma use moment matching in the reduced problem. Censoring or different channel weights also give working aggregate approximations; final scoring always uses the original responses and their specified model.

Temporal variation

High-order differences suppress the slowly varying mean. With a_k=(-1)^k{r\choose k} and stationary covariance c_h,

\frac{\operatorname{Var}(\Delta^r e)}{{2r\choose r}} =c_0+2\sum_{h=1}^r(-1)^h \frac{{2r\choose r+h}}{{2r\choose r}}c_h.

The identity follows by expanding the variance and using \sum_k a_k a_{k+h}=(-1)^h{2r\choose r+h}. With lag-one covariance this reduces to c_0-2r c_1/(r+1). Fourth and eighth differences estimate the marginal variance and lag-one covariance. The optional penalty multiplier \kappa=1+2\widehat\rho_1 approximates the local information loss for smooth directions, while the marginal scale stays in the component likelihood. This is a scalar working covariance adjustment. Composite-likelihood theory relates the broader correction to \operatorname{tr}(JH^{-1}), where H is sensitivity and J is score covariance (Varin, Reid & Firth, 2011).

References

Kidziński, Ł., Hui, F. K. C., Warton, D. I. & Hastie, T. J. (2022). Generalized matrix factorization: efficient algorithms for fitting generalized linear latent variable models to large data arrays.

Castiglione, C., Segers, A., Clement, L. & Risso, D. (2026). Stochastic gradient descent estimation of generalized matrix factorization models with application to single-cell RNA sequencing data.

Wang, L. & Carvalho, L. E. (2023). Deviance matrix factorization.

Li, X. & Tyagi, A. (2023). Block-active ADMM to minimize NMF with Bregman divergences.

Banerjee, A. et al. (2005). Clustering with Bregman divergences.

Witten, D. M. (2011). Classification and clustering of sequencing data using a Poisson model.

Lee, D. D. & Seung, H. S. (2000). Algorithms for non-negative matrix factorization.

Eilers, P. H. C. (2005). Unimodal smoothing.

Varin, C., Reid, N. & Firth, D. (2011). An overview of composite likelihood methods.

The likelihood distributions follow the R negative-binomial and Gamma parameterizations. The compound-count criterion, local log-concave fusion and their integration with L0glm’s grouped initialization are developments of this package.