A fast two-stage algorithm for non-negative matrix factorization in streaming data

Given that the linear independence constraint qualification (LICQ) holds for NMF problem, the KKT conditions [17] for the problem can be written as

We denote $x=vec(X^{T})$, $y=vec(Y)$, $r=vec(R^{T})$, $s=vec(S)$, where $vec$ transforms a matrix to a vector by expanding it by columns. Meanwhile, we define the inverse operations $mat(x)=X^{T}$, $mat(y)=Y$, $mat(r)=R^{T}$, $mat(s)=S$.

Applying Newton’s method to the nonlinear system, in the variables $x$, $y$, $r$, $s$, we obtain

with $\mu=0$, where $graX=vec(Y(XY-M)^{T})$, $graY=vec(X^{T}(XY-M))$, $e$ is a vector of ones, $Diag(x)$ constructs a diagonal matrix with the diagonal elements given by $x$,

$X_{i}$ is the $i$th row of $X$, and $Y_{j}$ is the $j$th column of $Y$.

Let $\mu$ be strictly positive, then the variables $x$, $y$, $r$ and $s$ are forced to take positive values. The trajectory $(x(\mu),y(\mu),r(\mu),s(\mu))$ is called the primal-dual central path. The variables are updated by

where $\alpha_{1}\in(0,\alpha_{1}^{\max}]$, $\alpha_{2}\in(0,\alpha_{2}^{\max}]$, and

with $\tau\in(0,1)$. The condition (7) is called the fraction to the boundary rule [25], which is used to prevent the variables from approaching their lower bounds of 0 too quickly. In this work, we choose $\tau=0.9$.

A predictor or probing strategy [25] can also be used to determine the parameter $\mu$. We calculate a predictor (affine scaling) direction

by setting $\mu=0$. We probe this direction by letting $(\alpha_{1}^{aff},\alpha_{2}^{aff})$ be the longest step lengths that can be taken along the affine scaling direction before violating the nonnegativity conditions $(x,y,r,s)\geq 0$. Explicit formulas for these step lengths are given by (7) with $\tau=1$. We then define $\mu^{aff}$ to be the value of complementarity along the (shortened) affine scaling step, that is,

and a heuristic choice of $\mu$ is defined as following: $\mu=\sigma\mu$ and

We propose a two-level nested loop algorithm for the interior point search. In the inner loop, the parameter $\mu$ is fixed. In the outer loop, we gradually reduce $\mu$ to 0. We use the following error function to break the inner loop, which is based on the perturbed KKT system:

To guarantee the global convergence of the algorithm, we apply a line search approach. First we consider the second derivation of the interior-point method associated with the barrier problem

We use the exact merit function as same as the barrier function, which can be formed by

In the algorithm, we utilize Armijo line search to make the merit function sufficiently decrease.

Considering the convexity of the Hessian, we approximate

by a positive definite matrix to guarantee the direction $(\Delta x,\Delta y)$ is always a descent direction of $\phi(x,y)$, so that such a line search can be implemented. Notice that the original NMF problem is a least square problem. Thus we can utilize the Hessian in the traditional Guass-Newton algorithm. The Hessian matrix is

and it is guaranteed to be positive semi-definite. Here

which is the first item of $C$. For further safeguarding, we add a diagonal matrix $\rho I$ to this Hessian, where $\rho$ is a small positive constant. Then we obtain the primal-dual direction by solving

By eliminating $\Delta r$ and $\Delta s$ in (11), we have

where

and $Diag(x)^{-1}e$ is simply represented by $x^{-1}$. We simplify the above formula by

where $B$ represents the coefficient matrix and $p$ represents the direction $(\Delta x^{T},\Delta y^{T})^{T}$. Since $B$ is positive definite, the inner product $p^{T}\nabla\phi(x,y)>0$, which means that $p$ is a descent direction.

To sum up all the approaches above, we present the whole algorithm of interior point method in (2). It contains two loops. The parameter $\mu$ is fixed in the inner loop. Due to the Armijo line search (13) in the inner loop, the stopping criterion of the inner loop can be satisfied in finite iterations. Further, the parameter $\mu$ and $\epsilon_{\mu}$ are reduced gradually, and by the definition of error function $E$, the solutions of the two-loop algorithm satisfy the KKT system (3) within the error $\epsilon_{TOL}$. In practice, the barrier stop tolerance can be defined as

The complete convergence theorem and its proof are given in (1).

In general, the computational cost in each iteration of the interior point method is usually the cubic power of its size, which is impractical for large scale problems. However, for the special NMF problem under consideration here, the amount of computation can be greatly reduced, so that the proposed method can be applied to practical problems involving streaming data. In the following, we analyze the computational cost of algorithm 2.

First of all, to compute the gradient in (1), $O(nmk)$ flops are needed.

The main computational cost of (2) is computing the primal dual direction (11). One can solve (12) to obtain $(\Delta x,\Delta y)$ first, and then compute $(\Delta r,\Delta s)$ within a low cost.

We rewrite (12) as

In order to minimize the computational cost, we first decompose $\bar{Q}_{1}$.

can be obtained by Cholesky factorization or eigenvalue decomposition. Since the matrix $\bar{Q}_{1}$ is composed of $n$ positive definite diagonal blocks with the size of k times k, one can obtain $P$ and $P^{-1}$ within $O(nk^{3})$ flops. (19) is equivalent to

Then we solve $\Delta y$ from

and compute $\Delta x$ by

We need $O(nmk^{2})$ flops for constructing $\bar{C}$ and $O(nmk^{3})$ flops for computing $\bar{C}P^{-1}$. By considering that each block of $\bar{C}$ is rank one, we can compute $\bar{C}P^{-1}$ within $O(nmk^{2})$ flops. The dominant computation is the computation of $(\bar{C}P^{-1})(\bar{C}P^{-1})^{T}$ which costs $O(nm^{2}k^{3})$ flops. If we consider that each block of $\bar{C}$ is rank one, it can be reduced to $O(nm^{2}k^{2})$ flops. When solving (20) by Cholesky factorization, since the size of the coefficient matrix is $mk$ by $mk$, the computational cost is $O(m^{3}k^{3})$. Other computations, like computing the right side of (20) and computing $\Delta x$, are $O(nmk^{2})$ flops.

To sum up, the computation cost of algorithm 2 in each iteration is $O((n+mk)m^{2}k^{2})$. Compared with the computation of gradient $O(nmk)$, the cost is no more than $O(mk^{2})$ times. Since $n\gg m\gg k$ and $k$ is usually very small in our streaming data, the computational complexity is completely acceptable.

The KKT conditions of (1) are given in (3). The definition of the error function $E(x,y,r,s;0)$ (10) measures the violation of the KKT conditions. Therefore, we set

to be the stopping criterion, where

The ANLS-based algorithms do not generate the dual variables $r$ and $s$. Here we give a reasonable definition that

In some cases, we also limit the maximum CPU time, in cases where some algorithms can not reach the given accuracy.

Zhao et al. [32] measured the X-ray diffraction data during the nucleation and growth of zeolite-supported Ag nanoparticles through reduction of Ag-exchanged mordenite (MOR), and processed the data with pair distribution function (PDF) measurement. In the field of chemistry, more and more people use mathematical tools to analyze their measured data. Chapman et al. [4] used principle component analysis (PCA) and some post-processing to analyze the data given in [32], and obtained 3 principle components. Since PDF is a distribution-type function, NMF may be intuitively more applicable.

We simply remove the negativity of the raw data by shifting each PDF up by the opposite of its original minimum value. Then we perform the NMF algorithms on the data. The size of the data is

and in the algorithms, $k$ is set to be 3 based on the analysis of [4].

We use the same initial point for each algorithm. The relationship between KKT violation $E$ and CPU time of the four algorithms are shown in Fig. 1. It can be seen that the performances of NeNMF, ANLS-BPP and our 2-STAGE algorithm are similar in the early stage, and the decline of error function value for QRPBB is the most obvious. After that, there is a rise of error function value in 2-STAGE, because our algorithm begins to enter the second stage and the variables change, and some shocks occur later, which is caused by the instability of the early iterations of the interior point method. After several iterations, 2-STAGE becomes stable and converges quickly to meet the termination criterion. Compared with NeNMF and ANLS-BPP, QRPBB is always more accurate. Even so, it can not achieve the given criterion in a short time.

Next, we generate 10 different random initial points. For each of them, all algorithms are implemented. The results are shown in Tab. I. Because none of the three algorithms compared with the 2-STAGE algorithm can meet the stopping criterion in a short time, we set the maximum running CPU time as 20 seconds. Tab. I presents the average, minimum and maximum CPU time, the KKT violation $E$ and the objective function value of each algorithm respectively. It can be seen that our algorithm always has the minimum KKT violation and the minimum objective function value. It has a great advantage in finding a high-precision solution within a short amount of time.

The Yale face database is widely tested for face recognition [3]. It contains 165 grayscale images of 15 individuals. There are 11 images per subject, one per different facial expression or configuration: lighting (center-light, left-light and right-light), with/without glasses, facial expressions (normal, happy, sad, sleepy, surprised, and wink). The original image size is $320\times 243$ pixels. To reduce the computational burden, the image size was reduced to $64\times 64$ pixels. Because the image of face also has some degree of continuity, we regard the data as streaming data to test our algorithm. The data size is $n=4096,m=165$.

We gradually increase $m$ to explore the influence of $m$ variation on the algorithm. We selected the face images of the first $i$ individuals, $i=1,2,4,8,15$, i.e. $m=11,22,44,88,165$. In the following tests, we fix $k=3$. We limit a maximum CPU time of 60 seconds for each algorithm and generate 10 different random initial points for each sample. We see from Fig. 2, where we take an example of $m=44$ and an example of $m=165$, that the first stage of 2-STAGE is not as efficient as other algorithms. The main reason is that the continuity of face data is not strong enough compared with streaming data, thus the active set changes frequently in the algorithm resulting large computational cost. In Fig. 2 (top), the second stage of 2-STAGE, interior point method, shows a fast local convergence, since the slope is significantly steeper than the others. However, in Fig. 2 (bottom), the slope of the second stage of 2-STAGE is similar with others. The reason is that when $m$ increases, the computational cost of interior point method increases faster than that of other algorithms. Therefore, we see from Tab. II that, as anticipated, the larger $m$ is, the less efficient 2-STAGE is compared with other algorithms.

In order to test more problems on different scales, we artificially synthesized some data for testing.

We note that the focus of this paper is on streaming data. Chapman et al.[4] use PCA method to study AgMor data. They found that there were three dominant components (singular values) of the data, which means that the original data is a rank-3 matrix plus noise. According to this characteristics, we construct the artificial data by the following methods. First, we decide the problem size $(n,m,k)$. Then, we generate a random matrix $X$ whose size is $n\times k$, and each element is uniformly distributed from 0 to 1. In the same way, we generate a random matrix $Y$ whose size is $k\times m$. Next, we compute $M=XY$ and add Gaussian noise, whose expectation is 0 and standard deviation is $0.1$, to each element of $M$. Finally, we change the negative elements in $M$ to zeros.

In this set of tests, each algorithm is limited to a maximum CPU time of 60 seconds. We generate 10 different random initial points for sample. The results are shown in Tab. III. From an average point of view, the performance of 2-STAGE is the best of the four algorithms, followed by QRPBB. Except for M = 50, k = 6, ANLS-BPP exceeds QRPBB. In some cases, other algorithms may be better than 2-STAGE. For example, when $m=100$ and $k=6$, the minimum time of QRPBB is much less than the maximum time of 2-STAGE. Generally speaking, the performance of 2Stage in getting high-precision solution truly stands out.