Why Penalize NMF Factors?
Nonnegative matrix factorization represents a nonnegative data matrix as \(X\approx WH\). With \(p\) features, \(n\) samples, and rank \(r\), use \(X\in\mathbb{R}_{\geq 0}^{p\times n}\), \(W\in\mathbb{R}_{\geq 0}^{p\times r}\), and \(H\in\mathbb{R}_{\geq 0}^{r\times n}\). A column of \(W\) describes a nonnegative feature profile for a component; a column of \(H\) gives that sample’s component scores. The Frobenius loss rewards accurate reconstruction, but by itself it does not favor factors with many zero entries.
An L1 penalty adds a cost proportional to the sum of factor entries. Because NMF factors are nonnegative, this sum is also the L1 norm. Penalizing \(W\) can make component profiles selective across features; penalizing \(H\) can make sample representations use fewer components. These are different scientific aims, not interchangeable descriptions of “sparse NMF.”
Gradients and Nonnegative Optimality
First hold \(H\) fixed and differentiate with respect to \(W\), then hold \(W\) fixed and differentiate with respect to \(H\). On the nonnegative domain, the penalty terms are linear sums of entries, so their gradients are constant. The resulting gradients are
where each \(\mathbf{1}\) is an all-ones matrix of the same dimensions as the gradient it appears in. The gradients describe the smooth objective on the nonnegative orthant, including its one-sided derivative at a zero entry. At a constrained minimum, each entry must satisfy a boundary condition: a positive entry has zero gradient, while a zero entry has nonnegative gradient. In compact form, for each factor \(A\) and corresponding gradient \(G\), these conditions are \(A\geq 0\), \(G\geq 0\), and \(A\odot G=0\), entrywise. They are necessary for a block minimum; jointly they do not make the factorization problem convex.
A single-coordinate calculation makes the sparsity mechanism explicit. If all other terms are held fixed, let \(z\) be the residual data for one coefficient \(w\), and let \(h\) be its fixed multiplier. The part of the objective involving \(w\geq 0\) is \(\frac12(z-wh)^2+\lambda_Ww\). For \(h\ne 0\), differentiation gives \(h^2w-hz+\lambda_W=0\). Enforcing nonnegativity yields a thresholded solution.
If \(h=0\), the objective is \(\frac12z^2+\lambda_Ww\), which does not depend on \(w\) when \(\lambda_W=0\), so every \(w\geq 0\) minimizes it. When \(\lambda_W>0\), it is minimized at \(w=0\). This separate case matters: the threshold formula has a zero denominator at \(h=0\). The corresponding calculation for a coefficient in \(H\) has the same form, with the roles of the factor and its multiplier exchanged.
Worked Example: A Single Penalized Coefficient
Suppose a fixed multiplier is \(h=2\), the residual data value is \(z=8\), and \(\lambda_W=1\). The coefficient minimizes \(\frac12(8-2w)^2+w\) over \(w\geq 0\). Since \(h\ne 0\), use the threshold formula:
\(w^*=\max(0,(2\cdot 8-1)/2^2)=\max(0,15/4)=3.75.\) The unpenalized coefficient would be \(hz/h^2=16/4=4\), so the penalty reduces it. If instead \(hz\leq\lambda_W\), the constrained solution would be zero. The zero is a consequence of the penalty overcoming the residual’s support for that coefficient, not evidence that the corresponding feature or sample is absent in reality.
Multiplicative Updates and Their Implications
For fixed \(H\), the gradient for \(W\) can be separated into a nonnegative part \(WHH^T+\lambda_W\mathbf{1}\) and a subtractive part \(XH^T\). A standard multiplicative update scales each current entry by the ratio of these parts. Alternating the same construction for \(H\) gives
Products, ratios, and \(\odot\) are entrywise where appropriate. Each update uses the current value of the other factor; in an alternating algorithm, the \(H\) update typically uses the newly updated \(W\). With positive denominators and nonnegative initialization, the updates preserve nonnegativity. They are commonly derived through auxiliary-function or majorization arguments for the Frobenius objective. They do not turn the joint optimization into a convex problem or guarantee a global minimum.
A zero entry remains zero under multiplication, even if its gradient would favor becoming positive. This is zero locking. Positive initialization avoids locking entries at the start, but does not prevent later entries from reaching zero. If a denominator is zero, the ratio is undefined: the update needs an explicit convention, such as retaining the entry or setting it to zero when the corresponding numerator is also zero. In implementations, check denominators and objective values rather than assuming every ratio is well-defined.
Worked Example: Penalizing Feature Profiles
Consider the hypothetical nonnegative matrix \(X=\begin{pmatrix}2&3\\2&2\end{pmatrix}\), a rank-one factor \(H=\begin{pmatrix}2&1\end{pmatrix}\), and initial \(W=\begin{pmatrix}1\\1\end{pmatrix}\). Use \(\lambda_W=0.5\) and \(\lambda_H=0\), so only feature profiles are penalized. The \(W\) update needs \(HH^T=2^2+1^2=5\) and
\(XH^T=\begin{pmatrix}2(2)+3(1)\\2(2)+2(1)\end{pmatrix}=\begin{pmatrix}7\\6\end{pmatrix}.\)
Thus the denominator is \(WHH^T+0.5\mathbf{1}=\begin{pmatrix}5.5\\5.5\end{pmatrix}\), and
\(W_{\mathrm{new}}=\begin{pmatrix}1\\1\end{pmatrix}\odot\begin{pmatrix}7/5.5\\6/5.5\end{pmatrix}\approx\begin{pmatrix}1.273\\1.091\end{pmatrix}.\)
Without the penalty, the denominators would be 5, giving \(W_{\mathrm{new}}=(7/5,6/5)^T=(1.400,1.200)^T\). The penalized update reduces both feature weights. This one iteration does not make either weight zero; sparsity means that sufficiently weak entries can be driven to zero over optimization, not that every positive penalty immediately creates zeros.
Which Factor Should Be Sparse?
A penalty on \(W\) targets feature profiles: within a component, many feature weights may be zero, so that component is associated with a narrower subset of measured features. A penalty on \(H\) targets sample representations: within a sample, many component scores may be zero, so that sample is represented using fewer components. Penalizing both can encourage both patterns. The factor being penalized determines the interpretation of the induced zeros.
Worked Example: Feature Sparsity Versus Sample Sparsity
For a coordinate update, hold all other components fixed and write \(z\) for the residual vector for one sample. A score \(h\geq 0\) for a fixed profile \(w\) minimizes \(\frac12\|z-wh\|_2^2+\lambda_Hh\). When \(\|w\|_2^2>0\), differentiating and applying the nonnegative constraint gives \(h^*=\max(0,(w^Tz-\lambda_H)/\|w\|_2^2)\). Take \(w=(2,1)^T\) and \(z=(5,2)^T\). With \(\lambda_H=2\), the score is \(\max(0,(2\cdot5+1\cdot2-2)/(2^2+1^2))=10/5=2\). For another sample with residual \(z=(1,0)^T\) and \(\lambda_H=3\), the score is \(\max(0,(2\cdot1+1\cdot0-3)/5)=0\). This penalty has made that sample’s representation sparse for this component.
For a feature-profile coordinate, hold \(H\) fixed and let \(y\) be the residual row across samples. A coefficient \(w\) for profile scores \(h\) minimizes \(\frac12\|y-wh\|_2^2+\lambda_Ww\), so \(w^*=\max(0,(h y^T-\lambda_W)/\|h\|_2^2)\) when \(\|h\|_2^2>0\). With \(h=(2,1)\), \(y=(3,1)\), and \(\lambda_W=2\), this is \(\max(0,(3\cdot2+1\cdot1-2)/5)=1\). If another feature has residual row \(y=(1,0)\), its coefficient is \(\max(0,(1\cdot2+0\cdot1-2)/5)=0\). Here the penalty has induced sparsity across feature weights in a profile, not across component scores for a sample.
Assumptions, Scale, and Practical Checks
These formulas concern nonnegative data modeled by the Frobenius reconstruction loss, a fixed rank, and factor entries constrained to be nonnegative. The gradients and coordinate solutions are exact for their stated objectives. Alternating optimization solves one factor block at a time; each block problem is convex when the other factor is fixed, but the joint problem is nonconvex. Results can depend on initialization and may be stationary solutions rather than global optima.
L1 penalties also interact with NMF’s scale indeterminacy. For any \(c>0\), replacing \(W\) by \(cW\) and \(H\) by \(H/c\) leaves \(WH\) unchanged, but changes the total penalty to \(\lambda_W c\|W\|_1+\lambda_H\|H\|_1/c\). When both penalties and factor norms are positive, this expression is minimized over \(c\) at \(c=\sqrt{\lambda_H\|H\|_1/(\lambda_W\|W\|_1)}\). Thus penalty magnitudes cannot be interpreted independently of factor scaling. If only one factor is penalized and the other can be rescaled freely, the penalty can be driven toward zero along this scaling path without changing the reconstruction. Normalization or other constraints can alter that issue, but also alter the model being fit.
A useful check is to compare the reconstruction loss and each penalty contribution separately, rather than reporting only their sum. If a larger \(\lambda_W\) makes nearly all feature weights vanish while reconstruction error rises sharply, the profiles may be over-penalized. If a larger \(\lambda_H\) makes most sample scores zero, inspect whether the resulting reconstruction remains adequate. Repeat fits from different initializations, check for zero locking, and verify the objective after updates. The unpenalized Frobenius NMF objective is the main reference when sparsity is not required; hard constraints on the number of nonzero entries are another option, but they require combinatorial choices and do not provide the same continuous shrinkage as L1 penalties.
Check Your Understanding
Try each question, then check the solutions below.
- In \(X\approx WH\), which factor should be penalized to encourage fewer nonzero component scores per sample?
- For \(h=3\), \(z=2\), and \(\lambda_W=1\), find the constrained scalar optimum.
- Why can a zero entry stay zero under a multiplicative update even if a nonzero value could improve the objective?
- What happens to the L1 penalty terms under \(W\mapsto cW,\ H\mapsto H/c\), and why does this matter?
Solutions. 1. Penalize \(H\), whose columns are sample representations. 2. Since \(h\ne0\), \(w^*=\max(0,(3\cdot2-1)/3^2)=5/9\approx0.556\). 3. Multiplication by the current zero entry keeps the updated entry at zero; this is zero locking. 4. The reconstruction is unchanged, but the penalty becomes \(\lambda_Wc\|W\|_1+\lambda_H\|H\|_1/c\), which varies with \(c\). Penalty interpretation therefore depends on how factor scale is set.