Alternating Nonnegative Least Squares

Why Alternate Between NNLS Problems?

For nonnegative data, non-negative matrix factorization (NMF) represents an $m\times n$ matrix $X$ approximately as the product of two nonnegative factors: an $m\times r$ matrix $W$ and an $r\times n$ matrix $H$. With the Frobenius-norm loss, the goal is to minimize the squared reconstruction error while requiring every factor entry to be nonnegative. Optimizing $W$ and $H$ simultaneously is a nonconvex problem, but fixing either factor makes the objective a convex quadratic in the other. Alternating nonnegative least squares (ANLS) uses this fact: it solves one factor block at a time using nonnegative least squares (NNLS).

The method is useful when accurate block solutions are important and the NNLS subproblems can be solved efficiently. It also makes the source of each update explicit: the current $W$ defines NNLS problems for the columns of $H$, and the current $H$ defines NNLS problems for the rows of $W$. An exact solution of a block subproblem cannot increase the reconstruction error, although it does not guarantee that the overall nonconvex problem has been solved globally.

Definition: Given a nonnegative data matrix $X$, ANLS seeks nonnegative factors $W$ and $H$ by repeatedly minimizing the Frobenius-norm objective over one factor while holding the other fixed: $$ \min_{W\geq 0,\ H\geq 0} \frac{1}{2}\lVert X-WH\rVert_F^2. $$ Here $W\geq 0$ and $H\geq 0$ mean entrywise nonnegativity, and $\lVert\cdot\rVert_F$ is the square root of the sum of squared matrix entries.

The Two Blockwise Subproblems

First hold $W$ fixed. Write $x_j$ for column $j$ of $X$, and $h_j$ for column $j$ of $H$. The squared Frobenius norm is the sum of squared column residuals, so the minimization over $H$ separates into $n$ NNLS problems:

\min_{H\geq 0}\frac{1}{2}\lVert X-WH\rVert_F^2 = \sum_{j=1}^{n}\min_{h_j\geq 0}\frac{1}{2}\lVert x_j-Wh_j\rVert_2^2.

Each column $h_j$ gives the nonnegative coefficients that approximate sample $j$ using the columns of $W$. Next hold $H$ fixed. Write $x_{i:}$ for row $i$ of $X$, and $w_i$ for the column vector whose transpose is row $i$ of $W$. The row residuals separate, giving:

\min_{W\geq 0}\frac{1}{2}\lVert X-WH\rVert_F^2 = \sum_{i=1}^{m}\min_{w_i\geq 0}\frac{1}{2}\lVert x_{i:}^{\mathsf T}-H^{\mathsf T}w_i\rVert_2^2.

Thus, the first block update solves NNLS with design matrix $W$, once for each column of $X$. The second solves NNLS with design matrix $H^{\mathsf T}$, once for each row of $X$. Because the subproblems are independent within a block, they can be solved separately or in parallel.

Result: For an NNLS problem $\min_{z\geq 0}\frac{1}{2}\lVert y-Az\rVert_2^2$, let $g=A^{\mathsf T}(Az-y)$ be the gradient. A feasible $z$ is optimal exactly when $g\geq 0$ and $z_jg_j=0$ for every coordinate $j$. These conditions say that an active zero coordinate cannot have a negative gradient pointing toward a feasible decrease, while each positive coordinate has zero gradient.

The NNLS objective is convex, so these conditions characterize a global minimizer of that subproblem. If $A$ has full column rank, the quadratic is strictly convex and the minimizer is unique. If the columns are linearly dependent, multiple minimizers may exist; the minimum residual is still well-defined, but a solver may return different coefficient vectors.

How NNLS Solvers Enter ANLS

An active-set solver divides coordinates into a free set, whose coefficients are allowed to be positive, and an active set, fixed at zero. It solves least squares on the free columns, checks whether any free coefficient violates nonnegativity, and checks the gradient signs at active coordinates. A violating coefficient may be moved to the active set; an active coordinate with a negative gradient may be released. The process continues until the NNLS optimality conditions hold. Particular implementations differ in how they choose and update the sets.

Coordinate descent instead updates one coefficient at a time, holding the others fixed. For a nonzero column $a_j$ of $A$, a full coordinate minimization is:

z_j\leftarrow \max\left(0,\ z_j+\frac{a_j^{\mathsf T}(y-Az)}{a_j^{\mathsf T}a_j}\right).

The residual in this expression is computed before the update, using the current $z$. The formula follows by writing the residual with all other coordinates fixed: changing $z_j$ adds a multiple of $a_j$, so the unconstrained one-dimensional quadratic has its minimum at the current coefficient plus the residual projection divided by $a_j^{\mathsf T}a_j$. Restricting the coefficient to be nonnegative replaces that minimum by its projection onto $[0,\infty)$. A zero column has no effect on the fitted values and must be handled separately rather than used as a denominator.

Worked Example: Solving Both NMF Blocks for Fixed Factors

Suppose $$ X=\begin{pmatrix}1&1\\1&2\end{pmatrix},\qquad W=\begin{pmatrix}1&0\\1&1\end{pmatrix}, $$ and update $H$ while keeping $W$ fixed. For the first column $x_1=(1,1)^{\mathsf T}$, solve $Wh_1=x_1$. The equations are $h_{11}=1$ and $h_{11}+h_{21}=1$, so $h_1=(1,0)^{\mathsf T}$, which is feasible and has zero residual.

For the second column $x_2=(1,2)^{\mathsf T}$, the equations are $h_{12}=1$ and $h_{12}+h_{22}=2$. Hence $h_2=(1,1)^{\mathsf T}$, also feasible with zero residual. Combining the columns gives $$ H=\begin{pmatrix}1&1\\0&1\end{pmatrix},\qquad WH=X. $$

Now hold this $H$ fixed and update the rows of $W$. For the first row, $x_{1:}^{\mathsf T}=(1,1)^{\mathsf T}$, and $H^{\mathsf T}=\begin{pmatrix}1&0\\1&1\end{pmatrix}$. Solving $H^{\mathsf T}w_1=x_{1:}^{\mathsf T}$ gives $w_1=(1,0)^{\mathsf T}$. For the second row, the target is $(1,2)^{\mathsf T}$, giving $w_2=(1,1)^{\mathsf T}$. Both row updates are nonnegative and exact. This example has an exact factorization, so both block losses are zero; in general, the NNLS fits need not be exact.

Comparing Active-Set and Coordinate-Descent Updates

The following small NNLS problem shows both the active-set solution and a sequence of coordinate-descent updates. It also illustrates why an intermediate iterate should not be mistaken for the final constrained optimum.

Worked Example: Active Set and Coordinate Descent on the Same NNLS Problem

Let $$ A=\begin{pmatrix}1&0\\1&1\end{pmatrix},\qquad y=\begin{pmatrix}1\\0.2\end{pmatrix}. $$ We minimize $\frac{1}{2}\lVert y-Az\rVert_2^2$ subject to $z\geq0$. If both coordinates were free, the equations $Az=y$ would give $z_1=1$ and $z_2=-0.8$. That unconstrained solution is infeasible, so an active-set solution places coordinate 2 at zero. With $z_2=0$, least squares on the first column $a_1=(1,1)^{\mathsf T}$ gives $$ z_1=\max\left(0,\frac{a_1^{\mathsf T}y}{a_1^{\mathsf T}a_1}\right) =\max\left(0,\frac{1(1)+1(0.2)}{1^2+1^2}\right)=0.6. $$

The resulting residual is $r=y-Az=(0.4,-0.4)^{\mathsf T}$. The gradient is $A^{\mathsf T}(Az-y)=-A^{\mathsf T}r=(-0,0.4)^{\mathsf T}$, or $(0,0.4)^{\mathsf T}$. Thus the positive first coordinate has gradient zero and the active second coordinate has a nonnegative gradient. The NNLS conditions hold, confirming $z=(0.6,0)^{\mathsf T}$ is optimal. Its residual sum of squares is $0.4^2+(-0.4)^2=0.32$, and the objective with the factor $1/2$ is $0.16$.

For coordinate descent, start at $z=(0,0.2)^{\mathsf T}$, updating coordinate 1 first. Holding $z_2=0.2$ fixed gives the residual before the update $$ y-a_2z_2=(1,0)^{\mathsf T}, \qquad z_1\leftarrow\max\left(0,0+\frac{(1,1)(1,0)^{\mathsf T}}{2}\right)=0.5. $$ The residual after this update is $y-A(0.5,0.2)^{\mathsf T}=(0.5,-0.5)^{\mathsf T}$. Now update coordinate 2, holding $z_1=0.5$ fixed. Its residual excluding coordinate 2 is $y-a_1(0.5)=(0.5,-0.3)^{\mathsf T}$, so $$ z_2\leftarrow\max\left(0,0.2+\frac{(0,1)(0.5,-0.3)^{\mathsf T}}{1}\right)=0. $$ On the next update of coordinate 1, the residual excluding coordinate 1 is $y=(1,0.2)^{\mathsf T}$, and $$ z_1\leftarrow\max\left(0,0+\frac{(1,1)(1,0.2)^{\mathsf T}}{2}\right)=0.6. $$ This reaches the active-set optimum. The example depends on using $y=(1,0.2)^{\mathsf T}$: with the stated starting value $z_2=0.2$, subtracting $a_2z_2=(0,0.2)^{\mathsf T}$ gives the first coordinate update's residual $(1,0)^{\mathsf T}$.

Comparison: An active-set method can identify a boundary solution directly by setting a coefficient to zero and checking the KKT conditions. Coordinate descent is simple and uses inexpensive one-coordinate updates, but its progress can depend on update order and the correlation among design columns. Both solve convex NNLS subproblems; neither makes the full NMF objective convex.

One Alternating Cycle and Its Guarantee

An ANLS cycle first minimizes over $H$ with $W$ fixed, then minimizes over $W$ with the updated $H$ fixed. The guarantee is a direct consequence of the definition of a minimum, not a guarantee that the algorithm reaches the best possible factorization.

Update $H$.
For each column $x_j$, solve $\min_{h_j\geq0}\frac{1}{2}\lVert x_j-Wh_j\rVert_2^2$. The total objective cannot increase because the old $H$ was feasible for this minimization.

Update $W$.
For each row $x_{i:}$, solve $\min_{w_i\geq0}\frac{1}{2}\lVert x_{i:}^{\mathsf T}-H^{\mathsf T}w_i\rVert_2^2$, using the updated $H$. Again the old $W$ is feasible, so the objective cannot increase.

Check progress and stopping criteria.
Record the objective after a full cycle and stop according to a prespecified criterion, such as a sufficiently small relative decrease over successive cycles, subject to an iteration limit.

More formally, let $F(W,H)=\frac{1}{2}\lVert X-WH\rVert_F^2$. If $H^{(t+1)}$ minimizes the $H$-block for $W^{(t)}$, and $W^{(t+1)}$ minimizes the $W$-block for $H^{(t+1)}$, then feasibility of the previous block values gives

F(W^{(t)},H^{(t+1)})\leq F(W^{(t)},H^{(t)}), \qquad F(W^{(t+1)},H^{(t+1)})\leq F(W^{(t)},H^{(t+1)}).

Since $F\geq0$, the sequence of objective values is nonincreasing and bounded below, and therefore has a limit. This argument does not imply that the factors converge, that a limit is globally optimal, or even that every stationary point is reached. The product $WH$ couples the two blocks, so alternating minimization can settle at a suboptimal solution.

Worked Example: An Exact Alternating Cycle

Take $$ X=\begin{pmatrix}2&4\\1&2\end{pmatrix},\qquad W^{(0)}=\begin{pmatrix}1\\1\end{pmatrix}. $$ This is a rank-one setting. With $W^{(0)}$ fixed, each column of $H$ has one nonnegative coefficient. For $x_1=(2,1)^{\mathsf T}$, the least-squares coefficient is $(1,1)x_1/((1,1)(1,1)^{\mathsf T})=3/2$. For $x_2=(4,2)^{\mathsf T}$, it is $6/2=3$. Both are nonnegative, so $$ H^{(1)}=\begin{pmatrix}1.5&3\end{pmatrix}. $$

Now each row of $W$ is a scalar NNLS problem with design row $(1.5,3)$. Its squared norm is $1.5^2+3^2=11.25$. For the first data row $(2,4)$, the coefficient is $(2(1.5)+4(3))/11.25=15/11.25=4/3$. For the second row $(1,2)$, it is $(1(1.5)+2(3))/11.25=7.5/11.25=2/3$. Both are nonnegative, yielding $W^{(1)}=(4/3,2/3)^{\mathsf T}$. Multiplication gives $W^{(1)}H^{(1)}=X$, so the final objective is zero. The exact fit follows from the hypothetical rank-one matrix used here; it is not a general ANLS guarantee.

Diagnostics, Alternatives, and Limitations

A falling objective is necessary evidence of progress, but it is not enough to diagnose a successful fit. Check that the reported $W$ and $H$ are nonnegative, calculate the reconstruction error independently, and inspect whether the NNLS solver met its stopping tolerance. For active-set solutions, verify that free-coordinate gradients are near zero and active-coordinate gradients are nonnegative. For coordinate descent, check whether repeated sweeps still change the objective and whether coordinate updates are small. A solver that stops its inner NNLS iterations early may make the outer objective decrease slowly or appear to stall.

Poor conditioning or strongly correlated columns of $W$ (or rows of $H$) can make NNLS coefficients unstable and coordinate descent slow. Rank deficiency can produce different coefficient vectors with nearly identical fitted values. Check for non-finite values, zero columns in the design matrix, sensitivity to initialization, and large discrepancies between the stated and independently recomputed residual norm. Scaling can also affect numerical behavior, so compare objective values on the same data scale and document any scaling used.

Active-set NNLS is often attractive when the number of factors is modest and the active constraints can be identified efficiently; it may solve a subproblem in relatively few steps, though each step can require a least-squares solve. Coordinate descent has simple updates and can be effective for many modest, separable problems, but correlated predictors can slow convergence. Multiplicative updates, covered in the earlier material on multiplicative updates, offer a different update rule but should not be treated as equivalent to exactly solving each NNLS block. Projected gradient methods are another alternative; their updates are based on gradient steps and projection rather than NNLS subproblem solutions.

ANLS remains limited by the nonconvex coupling of $W$ and $H$: different initial factors can lead to different local solutions. Exact block solutions may also be computationally costly for large matrices, while approximate NNLS solves require care because the simple monotonicity argument assumes each block objective is actually minimized (or at least decreased). ANLS is therefore a method for improving a chosen factorization, not a proof that its final factorization is uniquely determined or globally best.

Common Mistake

Wrong reasoning: “Each NNLS subproblem is convex, so ANLS must find the global NMF solution.”
Correction: Convexity applies only when one factor is held fixed. Jointly, $WH$ depends on both factors, making the full optimization nonconvex. Exact block minimization ensures the objective does not increase at each update; it does not rule out a suboptimal stopping point.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

When $W$ is fixed, how many NNLS problems are needed to update $H$, and what is the design matrix for each?
For an NNLS solution $z$, state the gradient conditions for a positive coordinate and for a zero coordinate.
For $A=(1,2)^{\mathsf T}$, $y=(2,4)^{\mathsf T}$, and a current scalar $z=0.5$, calculate one coordinate-descent update.
What does monotonic decrease of the ANLS objective establish, and what does it not establish?

Solutions. (1) There are $n$ NNLS problems, one for each column $x_j$ of $X$; each uses design matrix $W$. (2) With $g=A^{\mathsf T}(Az-y)$, a positive coordinate must have $g_j=0$, while a zero coordinate must have $g_j\geq0$. (3) The residual is $y-Az=(1,3)^{\mathsf T}$, and $a^{\mathsf T}a=1^2+2^2=5$. Thus $z\leftarrow\max(0,0.5+((1,2)(1,3)^{\mathsf T})/5)=\max(0,0.5+7/5)=1.9$. (4) With exact block minimization, the objective values are nonincreasing and bounded below by zero, so they have a limit. This does not establish global optimality or convergence of the factors themselves.