An Iterative Algorithm for Regularized Non-negative Matrix Factorizations

The term non-negative matrix factorization refers to decomposing a given matrix $\mathsf{Y}$ with non-negative elements approximately as the product $\mathsf{Y}\approx\mathsf{L}\mathsf{R}$, for comformable non-negative matrices $\mathsf{L}$ and $\mathsf{R}$ of smaller rank than $\mathsf{Y}$. This factorization represents a kind of dimensionality reduction. The quality of the approximation is measured by some “divergence” between $\mathsf{Y}$ and $\mathsf{L}\mathsf{R}$, or some norm of their difference. Perhaps the most common is the Frobenius (or elementwise $\ell_{2}$) distance between them. Under this objective, non-negative matrix factorization can be expressed as the optimization problem:

where $\operatorname{tr}\left(\mathsf{M}\right)$ is the trace of the matrix $\mathsf{M}$, and $\mathsf{A}\geq 0$ means the elements of $\mathsf{A}$ are non-negative. [10]

Lee and Seung described an iterative approach to this problem, as well as for a related optimization with a different objective. [3] Their algorithm starts with some unspecified initial iterates ${\mathsf{L}}^{\left(0\right)}$ and ${\mathsf{R}}^{\left(0\right)}$, then iteratively improves ${\mathsf{L}}^{\left(k\right)}$ and ${\mathsf{R}}^{\left(k\right)}$ in turn. That is, given ${\mathsf{L}}^{\left(k\right)}$ and ${\mathsf{R}}^{\left(k\right)}$, an update is made to arrive at ${\mathsf{L}}^{\left(k+1\right)}$; then using ${\mathsf{L}}^{\left(k+1\right)}$ and ${\mathsf{R}}^{\left(k\right)}$, one computes ${\mathsf{R}}^{\left(k+1\right)}$, and the process repeats. The iterative update for ${\mathsf{L}}^{\left(k+1\right)}$ is

We use $\odot$ to mean the Hadamard, or ‘elementwise’, multiplication, and $\mathit{\oslash}$ to mean Hadamard division. The iterative update for ${\mathsf{R}}^{\left(k+1\right)}$ is defined mutatis mutandis. As should be apparent from inspection, this iterative update maintains non-negativity of estimates. That is, if $\mathsf{Y}$, ${\mathsf{L}}^{\left(k\right)}$, and ${\mathsf{R}}^{\left(k\right)}$ have only positive elements, then ${\mathsf{L}}^{\left(k+1\right)}$ is finite and has only positive elements. It is not clear what should be done if elements of the denominator term ${\mathsf{L}}^{\left(k\right)}$ ${\mathsf{R}}^{\left(k\right)}{{{\mathsf{R}}^{\left(k\right)}}^{\top}}$ are zero, and whether and how this can be avoided.

Gonzalez and Zhang describe a modification of the update to accelerate convergence. [2] Their accelerated algorithm updates ${\mathsf{L}}^{\left(k\right)}$ row by row, and ${\mathsf{R}}^{\left(k\right)}$ column by column, using the step length modification of Merritt and Zhang [5].

One attractive feature of the Lee and Seung iterative update is its sheer simplicity. In a high-level computer language with a rich set of matrix operations, such as Julia or octave, the iterative step can be expressed in a single line of code; the entire algorithm could be expressed in perhaps 20 lines of code. Conceivably this algorithm could even be implemented in challenged computing environments, like within a SQL query or a spreadsheet macro language. Simplicity makes this algorithm an obvious choice for experimentation without requiring an all-purpose constrained quadratic optimization solver, or other complicated software.

We generalize the optimization problem of 1 to include a weighted norm, regularization terms, and orthogonality penalties. For a given $\mathsf{Y}$, and non-negative weighting matrices ${\mathsf{W}}_{0,R}$, ${\mathsf{W}}_{0,C}$, ${\mathsf{W}}_{1,L}$, ${\mathsf{W}}_{1,R}$, ${\mathsf{W}}_{2,R,L,j}$, ${\mathsf{W}}_{2,C,L,j}$, ${\mathsf{W}}_{2,R,R,j}$, and ${\mathsf{W}}_{2,C,R,j}$, find matrices $\mathsf{L}$ and $\mathsf{R}$ to minimize

Essentially we are trying to model $\mathsf{Y}$ as $\mathsf{L}\mathsf{R}$. Our loss function is quadratic in the error $\mathsf{Y}-\mathsf{L}\mathsf{R}$, with weights ${\mathsf{W}}_{0,R}$ and ${\mathsf{W}}_{0,C}$, and $\ell_{1}$ and $\ell_{2}$ regularization via the terms with ${\mathsf{W}}_{1}$ and ${\mathsf{W}}_{2,R,j}$ and ${\mathsf{W}}_{2,C,j}$. This optimization problem includes that of 1 as a special case.

We note that a non-orthogonality penalty term on the $\mathsf{X}$ can be expressed as

Thus a non-orthogonality penalty term can be expressed with the ${\mathsf{W}}_{2,R,j}$ and ${\mathsf{W}}_{2,C,j}$. Indeed the motivation for including a sum of the $\ell_{2}$ terms is to allow for both straightforward regularization, where ${\mathsf{W}}_{2,R,j}=\mathsf{I}$ and ${\mathsf{W}}_{2,C,j}=c\mathsf{I}$, and orthogonality terms. For simplicity of notation, in the discussion that follows we will often treat the $\ell_{2}$ terms as if there were only one of them, i.e., $J=1$, and omit the $j$ subscripts. Restoring the more general case consists only of expanding to a full sum of $J$ terms.

This problem formulation is appropriate for power users who can meaningfully select the various weighting matrices. And while it is easy to conceive of a use for general forms for ${\mathsf{W}}_{0,R}$ and ${\mathsf{W}}_{0,C}$ (the rows of $\mathsf{Y}$ have different importances, or the columns of $\mathsf{Y}$ are observed with different accuracies, for example), most users would probably prefer simpler formulations for the ${\mathsf{W}}_{1,\cdot}$ and ${\mathsf{W}}_{2,\cdot,\cdot,j}$ matrices. To achieve something like an elasticnet factorization with non-orthogonality penalties, one should set $J=2$ and take

where here $\mathsf{1}$ stands in for an all-ones matrix of the appropriate size. All the $\lambda$ and $\gamma$ parameters must be non-negative, and the identity matrices are all appropriately sized. This is the reduced form of the problem that depends mostly on a few scalars. A simple form for the weighted error is to take ${\mathsf{W}}_{0,R}$ and ${\mathsf{W}}_{0,C}$ to be diagonal matrices,

Depending on the regularization terms, the solution to the problem of 3 is unlikely to be unique. Certain the solution to 1 is not unique, since if $\mathsf{Q}$ is a permutation matrix of the appropriate size, then $\mathsf{L}\mathsf{Q}$ and ${{\mathsf{Q}}^{-1}}\mathsf{R}$ are another set of non-negative matrices with the same objective value. Perhaps the regularization terms can help avoid these kinds of ambiguity in the solution.

Theorem 2.3 only tells us that the sequence $\phi\left({\boldsymbol{x}}^{\left(k\right)}\right)$ is non-decreasing. It does not guarantee convergence to a global, or even a local, minimum. In fact, if we only restrict $\boldsymbol{d}$ to be non-positive, it is easy to construct cases where ${\boldsymbol{x}}^{\left(k\right)}$ will not converge to the constrained optimum. For if an element of $\boldsymbol{d}$ is zero, then that same element of ${\boldsymbol{x}}^{\left(k\right)}$ will be zero for all $k>0$. However, the optimum may occur for ${\boldsymbol{x}}$ with a non-zero value for that element. (Certainly the unconstrained solution $-{{\mathsf{G}}^{-1}}{\boldsymbol{d}}$ may have a non-zero value for that element.)

Here we analyze the Lee and Seung update as a more traditional additive step, rather than a multiplicative step. Expressing the iteration as an additive step, we have

where

This may not be the optimal length step in the direction of ${\boldsymbol{h}}^{\left(k+1\right)}$. If we instead take

then $\phi\left({\boldsymbol{x}}^{\left(k+1\right)}\right)$ is a quadratic function of $c$ with optimum at

However, if $c^{*}>1$ elements of ${\boldsymbol{x}}^{\left(k\right)}+c^{*}{\boldsymbol{h}}^{\left(k+1\right)}$ could be negative. One would instead have to select a $c$ that results in a strictly non-negative ${\boldsymbol{x}}^{\left(k+1\right)}$. However, if $c^{*}<1$, then the original algorithm would have overshot the optimum.

We note that if $c^{*}\leq 1$ then $c^{\prime}=c^{*}$, since non-negativity is sustained for smaller step sizes.

This iterative update was originally described by Merritt and Zhang for solution of “Totally Nonnegative Least Squares” problems. [5] Their algorithm has guaranteed convergence for minimization of $\left\|{\mathsf{A}{\boldsymbol{x}}-\boldsymbol{b}}\right\|_{2}^{2}$, without the weighting or regularization terms, under certain conditions. They assume strict positivity of $\boldsymbol{b}$, which we cannot easily express in terms of $\mathsf{G}$ and $\boldsymbol{d}$, and it is not obvious that their proof can be directly used to guarantee convergence. Without guarantees of convergence, we express their method in our nomenclature as Algorithm 2, computing the step as ${\boldsymbol{h}}^{\left(k+1\right)}\leftarrow-\nabla\phi\left({\boldsymbol{x}}^{\left(k\right)}\right)\odot{\boldsymbol{x}}^{\left(k\right)}\mathit{\oslash}\left(\mathsf{G}{\boldsymbol{x}}^{\left(k\right)}\right)$. The vague language around selecting the step length $\tau_{k}$ is due to Merritt and Zhang; presumably one can choose it randomly, or always take $\tau_{k}=\left(\tau+1\right)/2$.