Conic-Optimization Based Algorithms for
Nonnegative Matrix Factorization

Nonnegative matrix factorization (NMF) is the problem of approximating a given nonnegative matrix, $V\in\mathbb{R}_{+}^{F\times N}$, as the product of two smaller nonnegative matrices, $W\in\mathbb{R}_{+}^{F\times K}$ and $H\in\mathbb{R}_{+}^{K\times N}$, where $K$ is a given parameter known as the factorization rank. One aims at finding the best approximation, that is, the one that minimizes the discrepancy between $V$ and the product $WH$, often measured by the Frobenius norm of their difference, $\left\|V-WH\right\|_{F}$. Despite the fact that NMF is NP-hard in general [30] (see also [27]), it has been used successfully in many domains such as probability, geoscience, medical imaging, computational geometry, combinatorial optimization, analytical chemistry, and machine learning; see [11, 13] and the references therein. Many local optimization schemes have been developed to compute NMFs. They aim to identify local minima or stationary points of optimization problems that minimize the discrepancy between $V$ and the approximation $WH$. Most of these iterative algorithms rely on a two-block coordinate descent scheme that consists in (approximatively) optimizing alternatively over $W$ with $H$ fixed, and vice-versa; see [5, 13] and the references therein. In this paper, we are interested in computing high-quality local minima for the NMF optimization problems without relying on the block coordinate descent (BCD) framework. We will perform the optimization over $W$ and $H$ jointly. Moreover, our focus is on finding exact NMFs, that is, computing nonnegative factors $W$ and $H$ such that $V=WH$, although our approaches can be used to find approximate NMFs as well.

The minimum factorization rank $K$ for which an exact NMF exists is called the nonnegative rank of $V$ and is denoted $\operatorname{rank}_{+}(V)$, we have

The computation of the nonnegative rank is NP-hard [30], and is a research topic on its own; see [13, Chapter 3][9, 8] and the references therein for recent progress on this question.

In this section we propose two new formulations for NMF, where the feasible set is represented using the exponential cone (Section 2.1) and the second-order cone (Section 2.2).

In this section, we present an iterative algorithm to tackle problems (5) and (9). Both problems can be written as the minimization of a concave function $\Phi$ over a convex set denoted by $Q$. Note that $Q$ designates either the feasible set of (5) or the feasible set of (9). We perform this minimization by solving a sequence of simpler problems in which the objective function is replaced by its linearization constructed at the current solution $(U,T)$. Let us denote $Z^{(i)}=(U^{(i)},T^{(i)})$ the $i$th iterate of our algorithm. At each iteration $i$, we update $Z$ as follows:

where $\Phi$ is the objective function of (5) or (9). Since the objective of (14) is linear in $Z$, the subproblems become convex. Moreover they are particular structured conic optimization problems. In this paper, we use the MOSEK software [25] to solve each successive problem (14) with an IPM. Algorithm 1 summarizes our proposed method to tackle (5) and (9).

To initialize $U$ and $T$, we chose to randomly initialize $W$ and $H$ (using the uniform distribution in the interval $[0,1]$ for each entry of $W$ and $H$) and apply the two changes of variables, $G(.)$, to compute the initializations for $U$ and $T$.

In this paper, we use a tolerance for the relative error equal to $10^{-6}$, that is, we assume that an exact NMF $(W,H)$ is found for an input matrix $V$ as soon as $\frac{\left\|V-WH\right\|_{F}}{\left\|V\right\|_{F}}\leq 10^{-6}$, as done in [29].

The main algorithm integrates a procedure that automatically updates the optimization problems in the case subsets of entries of the solution tend to zero. Indeed, due to numerical limitations of the solver, the required level of accuracy cannot be reached in some numerical tests even if the solution is close to convergence. This procedure is referred to as Sparsity Patterns Integration (SPI) and is detailed in Appendix B. Note that the update of the optimization problem is computationally costly, in particular the update of the matrix of coefficients defining the constraints. Hence, SPI is triggered twice; at 80% and 95% of the maximum number of iterations. This practical choice has been motivated by numerical experiments that showed that two activations in the final iterations are sufficient to reach the tolerance error when the current solution is close enough to a high-accuracy local optimum. However, for exponential cones, in some numerical tests, the relative error can be stuck in the interval $[10^{-4},10^{-5}]$. In this context only, a final refinement step further improves the output of the main algorithm using the state-of-the-art accelerated HALS algorithm, an exact BCD method for NMF, from [14] to go below the $10^{-6}$ tolerance.

In Section 4 , we discuss the convergence guarantees for Algorithm 1.

Let us focus on the following optimization problem

where $\Phi$ is a concave continuously differentiable function over the domain $Q$ which is assumed to be convex and compact. Let us first describe the convergence of the sequence of objective function values $\{\Phi(Z^{(i)})\}$ obtained with Algorithm 1. Since $\Phi(Z)$ is concave, its linearization around the current iterate $Z^{(i)}$ provides an upper approximation, that is,

This upper bound is tight at the current iterate and is exactly minimized over the feasible set $Q$ at each iteration. Hence Algorithm 1 is a majorization-minimization algorithm. This implies that $\Phi(Z)$ is nonincreasing under the updates of Algorithm 1 and since $\Phi(Z)$ is bounded below on $Q$, by construction, the sequence of objective function values $\{\Phi(Z^{(i)})\}$ converges.

We now focus on the convergence analysis of the sequence of iterates $\{Z^{(i)}\}$ generated by Algorithm 1, in particular convergence to a stationary point. To achieve this goal, we first recall some basics about the Frank-Wolfe (FW) algorithm. The FW algorithm [10] is a popular first-order method to solve (15) that relies on the ability to compute efficiently the so-called Linear Minimization Oracle (LMO), that is, $LMO(D):=\underset{Z\in Q}{\text{argmin }}\langle D,Z\rangle$ where $D$ denotes some search direction. The FW algorithm with adaptive step size is given in Algorithm 2. A step of this algorithm can be briefly summarized as follows: at a current iterate $Z^{(i)}$, the algorithm considers the first-order model of the objective function (its linearization), and moves towards a minimizer of this linear function, computed on the same domain $Q$.

Algorithm 1 is a particular case of Algorithm 2 for which $\tau^{(i)}=1$ for all $i$. In the last decade, FW-based methods have regained interest in many fields, mainly driven by their good scalability and the crucial property that Algorithm 2 maintains its iterates as a convex combination of a few extreme points. This results in the generation of sparse and low-rank solutions since at most one additional extreme point of the set $Q$ is added to the convex combination at each iteration. More details and insights about the later observations can be found in [6, 17]. FW algorithms have been recently studied in terms of convergence guarantees for the minimization of various classes of functions over convex sets, such as convex functions with Lipschitz continous gradient [18], and non-convex differentiable functions [22]. However, we are not aware of any convergence rates proven for Algorithm 2 when solving (15) assuming only concavity of the objective. To derive rates in the concave setting, we need to define a measure of stationarity for our iterates. In this paper, we consider the so-called FW gap of $\Phi$ at $Z^{(i)}$ defined as follows:

This quantity is standard in the analysis of FW algorithms, see [18, 22] and the references therein. A point $Z^{(i)}$ is a stationary point for the constrained optimization problem (15) if and only if $\mu^{(i)}=0$. Moreover, the FW gap

• •

  provides a lower bound on the accuracy: $0\leq\mu^{(i)}\leq\Phi(Z^{(i)})-\Phi^{*}$ for all $i$, where $\Phi^{*}:=\min_{Z\in Q}\Phi(Z)$,

• •

  is affine invariant, that is, it is invariant with respect to an affine transformation of the domain $Q$ in problem (15) [18], and

• •

  is not tied to any specific choice of norms, unlike criteria such as $\left\|\nabla\Phi(Z^{(i)})\right\|$.

Let us provide a convergence rate for the FW algorithm (Algorithm 2).

Theorem 4.1 shows that it takes at most $\mathcal{O}(\frac{1}{\epsilon})$ iterations to find an approximate stationary point with gap smaller than $\epsilon$. Note that Theorem 4.1 requires $\min_{i}\tau^{(i)}>0$ and, for example, the standard choice $\tau^{(i)}=\frac{2}{i+1}$ does not satisfy this assumption.

In this section, Algorithm 1, using both formulations (5) and (9), is tested for the computation of exact NMF for particular classes of matrices usually considered in the literature: (1) $10$-by-$10$ matrices randomly generated with nonnegative rank $5$ (each matrix is generated by multiplying two random rank-5 nonnegative matrices), (2) four $6$-by-$6$ infinitesimally rigid factorizations with nonnegative rank $5$ [20], denoted $V_{\text{inf}i}$ for $i=1,2,3,4$, and (3) four $5$-by-$5$ slack matrices corresponding to nested hexagons, denoted $V_{a=x}$, with nonnegative ranks $3,4,5,5$ depending on a parameter $x=2,3,4,+\infty$, respectively. These matrices are described in more details in the Appendix A.

In order to make Algorithm 1 practically more effective, we incorporate two improvements in the context of Exact NMF:

1. 1.

   Sparsity patterns integration: in NMF, entries of $W$ and $H$ are expected to be equal to zero at optimality. Hence, when some entries of $W$ or $H$ are sufficiently close to zero, fixing them to zero for all remaining iterations reduces the number of variables and hence accelerates the subsequent iterations of Algorithm 1; see Appendix B for the details.

2. 2.

   Final refinement: once a solution is generated by Algorithm 1 for the formulation based on exponential cones, it can typically be slightly improved by applying a standard NMF algorithm (we use a few iterations of A-HALS [14]). This is due to the bound constraints on $W$ and $H$; see Section 4.1.

For each of the matrices, we run our Algorithm 1 for 750 iterations for nested hexagons and random matrices and 3000 iterations for infinitesimally rigid matrices, and compare with the state-of-the-art algorithms from [29] with Multi-Start 1 heuristic ”ms1” and the Rank-by-rank heuristic ”rbr”. For each method, we run 100 initializations with SPARSE10, as recommended in [29], with target precision $10^{-6}$. Note that a different random matrix is generated each time for the experiments on random matrices, following the procedure described in the Appendix A. All tests are preformed using Matlab R2021b on a laptop Intel CORE i7-11800H CPU @2.3GHz 16GB RAM. The code is available from https://bit.ly/3FqMqhD.

Table 1 reports the number of successes over 100 attempts to compute the exact NMF of the input matrices, where the success is defined as obtaining a solution where $\frac{\left\|V-WH\right\|_{F}}{\left\|V\right\|_{F}}$ is below the target precision, namely $10^{-6}$.