The Diagonalized Newton Algorithm for Nonnegative Matrix Factorization
Authors: Hugo Van hamme
ICLR 2013

Abstract: Non-negative matrix factorization (NMF) has become a popular machine learning approach to many problems in text mining, speech and image processing, bio-informatics and seismic data analysis to name a few. In NMF, a matrix of non-negative data is approximated by the low-rank product of two matrices with non-negative entries. In this paper, the approximation quality is measured by the Kullback-Leibler divergence between the data and its low-rank reconstruction. The existence of the simple multiplicative update (MU) algorithm for computing the matrix factors has contributed to the success of NMF. Despite the availability of algorithms showing faster convergence, MU remains popular due to its simplicity. In this paper, a diagonalized Newton algorithm (DNA) is proposed showing faster convergence while the implementation remains simple and suitable for high-rank problems. The DNA algorithm is applied to various publicly available data sets, showing a substantial speed-up on modern hardware.

 
The Diagonalized Newton Algorithm for Non-
negative Matrix Factorization 
 
 Hugo Van hamme 1 
 University of Leuven, dept. ESAT 2 
 Kasteelpark Arenberg 10 – bus 2441, 3001 Leuven, B elgium 3 
 hugo.vanhamme@esat.kuleuven.be 4 
Abstract 5 
Non-negative matrix factorization (NMF) has become a popular machine 6 
learning approach to many problems in text mining, speech and image 7 
processing, bio-informatics and seismic data analys is to name a few. In 8 
NMF, a matrix of non-negative data is approximated by the low-rank 9 
product of two matrices with non-negative entries. In this paper, the 10 
approximation quality is measured by the Kullback-L eibler divergence 11 
between the data and its low-rank reconstruction. T he existence of the 12 
simple multiplicative update (MU) algorithm for com puting the matrix 13 
factors has contributed to the success of NMF. Desp ite the availability of 14 
algorithms showing faster convergence, MU remains p opular due to its 15 
simplicity. In this paper, a diagonalized Newton al gorithm (DNA) is 16 
proposed showing faster convergence while the imple mentation remains 17 
simple and suitable for high-rank problems. The DNA  algorithm is applied 18 
to various publicly available data sets, showing a substantial speed-up on 19 
modern hardware. 20 
 21 
1 Introduction  22 
Non-negative matrix factorization (NMF) denotes the  process of factorizing a N×T  data 23 
matrix V of non-negative real numbers into the product of a  N×R  matrix W and a R×T 24 
matrix H, where both W and H contain only non-negative real numbers. Taking a c olumn-25 
wise view of the data, i.e. each of the T columns of V is a sample of N-dimensional vector 26 
data, the factorization expresses each sample as a (weighted) addition  of columns of W, 27 
which can hence be interpreted as the R parts that make up the data [1]. Hence, NMF can be  28 
used to learn data representations  from samples. In [2], speaker representations are learnt 29 
from spectral data using NMF and subsequently appli ed to separate their signals. Another 30 
example in speech processing is [3] and [4], where phone representations are learnt using a 31 
convolutional extention of NMF. In [5], time-freque ncy representations reminiscent of 32 
formant traces are learnt from speech using NMF. In  [6], NMF is used to learn acoustic 33 
representations for words in a vocabulary acquisiti on and recognition task. Applied to image 34 
processing, local features are learnt from examples  with NMF in order to represent human 35 
faces in a detection task [7].  36 
In this paper, the metric to measure the closeness of reconstruction Z = WH  to its target V is 37 
measured by their Kullback-Leibler divergence: 38 
 39 
qg18R⇠qg3112qg3113 qg4⇠⇠⇠VV VVquni112⊘quni1119ZZ ZZqg4⇠⇠∇=qg3R33qg18∇4qg3141qg314∇ qg18⇠4qg18⇠∇qg18R9 qg343⇠qg18∇4qg3141qg314∇ 
qg18∇8qg3141qg314∇ 
qg3441−qg3R33qg18∇4qg3141qg314∇ 
qg3141quni112⊘qg314∇
+qg3R33qg18∇8qg3141qg314∇ 
qg3141quni112⊘qg314∇qg3141quni112⊘qg314∇
 (1)
Given a data matrix V, the matrix factors W and H are then found by minimizing cost 40 
function (1), which yields the maximum likelihood e stimate if the data are drawn from a 41 
Poisson distribution. The multiplicative updates (M U) algorithm proposed in [1] solves 42 
exactly this problem in an iterative manner. Its si mplicity and the availability of many 43 
implementations make it a popular algorithm to date  to solve NMF problems. However, there 44 
are some drawbacks to the algorithm. Firstly, it on ly converges locally and is not guaranteed 45 
to yield the global minimum of the cost function. I t is hence sensitive to the choice of the 46 
initial guesses for W and H. Secondly, MU is very slow to converge. The goal  of this paper  47 
is to speed up  the convergence while the local convergence proper ty is retained. The 48 
resulting Diagonalized Newton Algorithm (DNA) uses only simple element-wise operations , 49 
such that its implementation requires only a few te ns of lines of code, while memory 50 
requirements and computational efforts for a single  iteration are about the double of an MU 51 
update. 52 
The faster convergence rate is obtained by applying  Newton’s method to minimize 53 
dKL (V,WH ) over W and H in alternation. Newton updates have been explored for the Frobe-54 
nius norm to measure the distance between V and Z in e.g. [8]-[13]. Specifically, in [11] a 55 
diagonal Newton method is applied Frobenius norms. For the Kullback-Leibler divergence, 56 
fewer studies are available. Since each optimizatio n problem is multivariate, Newton updates 57 
typically imply solving sets of linear equations in  each iteration. In [16], the Hessian is 58 
reduced by refraining from second order updates for  the parameters close to zero. In [17], 59 
Newton updates are applied per coordinate, but in a  cyclic order, which is troublesome for 60 
GPU implementations. In the proposed method, matrix  inversion is avoided by diagonalizing 61 
the Hessian matrix. The resulting updates resemble the ones derived in [18] to the extent that 62 
they involve second order derivatives. Important di fferences are that [18] involves the non-63 
negative k-residuals hence requiring flooring to zero. Of cou rse, the diagonal approximation 64 
may affect the convergence rate adversely. Also, Ne wton algorithms only show (quadratic) 65 
convergence when the estimate is sufficiently close  to the local minimum and therefore need 66 
damping, e.g. Levenberg-Marquardt as in [14], or st ep size control as in [15] and [16]. In 67 
DNA, these convergence issues are addressed by comp uting both the MU and Newton 68 
solutions and selecting the one leading to the grea test reduction in dKL (V,Z). Hence, since 69 
the cost is non-decreasing under MU, it will also b e under DNA updates. This robust safety 70 
net can be constructed fairly efficiently because t he quantities required to compute the MU 71 
have already been computed in the Newton update. Th e net result is that DNA iterations are 72 
only about two to three times as slow as MU iterati ons, both on a CPU and on a GPU. The 73 
experimental analysis shows that the increased conv ergence rate generally dominates over 74 
the increased cost per iteration such that overall balance is positive and can lead to speedups 75 
of up to a factor 6. 76 
2 NMF formulation 77 
To induce sparsity on the matrix factors, the KL-di vergence is often regularized, i.e. one 78 
seeks to minimize: 79 
quni1119qg18R⇠qg3112qg3113 qg4⇠⇠⇠VV VVquni112⊘quni1119WW WWquni1119qg1∇82qg4⇠⇠∇+ρqg3R33qg18∇Rqg3141qg314R 
qg3141quni112⊘qg314R
+λqg3R33ℎqg314Rqg314∇ 
qg314Rquni112⊘qg314∇
 (2)
subject to non-negativity constraints on all entrie s of W and H. Here, ρ and λ are non-negative re-80 
gularization parameters. 81 
Minimizing the regularized KL-divergence (2) can be  achieved by alternating updates of W and H 82 
for which the cost is non-increasing. The updates for this form of block coordinate descent are: 83 
qg1∇82quni1119quni1119quni1119quni1119quni2191quni1119quni11⇠1quni11∇2quni11⇠∇quni1119quni11⇠⊙quni11⇠9quni11⇠E 
qg1∇82qg4R93≥qg2∇∇∇quni1119
quni1119qg343∇qg18R⇠qg3112qg3113 qg4⇠⇠⇠VV VVquni112⊘quni1119WW WWquni1119qg1∇82qg4R93qg4⇠⇠∇+λqg3R33ℎqg314Rqg314∇ 
qg4R93
qg314Rquni112⊘qg314∇
qg3441 (3)
qg1∇9∇quni1119quni1119quni1119quni1119quni2191quni1119quni11⇠1quni11∇2quni11⇠∇quni1119quni11⇠⊙quni11⇠9quni11⇠E 
qg1∇9∇qg4R93≥qg2∇∇∇quni1119
quni1119qg343∇qg18R⇠qg3112qg3113 qg4⇠⇠⇠VV VVquni112⊘quni1119qg1∇9∇qg4R93qg1∇82qg4⇠⇠∇+ρqg3R33qg18∇Rqg3141qg314R qg4R93
qg3141quni112⊘qg314R
qg3441 (4)
Because of the symmetry property dKL (V,WH ) = dKL (Vt,HtWt), where superscript- t denotes 84 
matrix transpose,  it suffices to consider only the  update on H. Furthermore, because of the 85 
summation over all columns in (1), minimization (3)  splits up into T independent optimiza-86 
tion problems. Let v denote any column of V and let h denote the corresponding column of H, 87 
then the following is the core minimization problem to be considered: 88 
qg1818quni1119quni1119quni1119quni1119quni2191quni1119quni11⇠1quni11∇2quni11⇠∇quni1119quni11⇠⊙quni11⇠9quni11⇠E 
qg1818qg4R93≥qg2∇∇∇quni1119
quni1119qg4⇠⇠⇠qg18R⇠qg3112qg3113 qg4⇠⇠⇠VV VVquni112⊘quni1119WW WWquni1119qg1818qg4R93qg4⇠⇠∇+λqg2∇∇8qg314∇qg1818qg4R93qg4⇠⇠∇ (5)
where 1 denotes a vector of ones of appropriate length. Th e solution of (5) should satisfy the KKT 89 
conditions, i.e. for all r with hr > 0 90 
qg2134qg4⇠⇠⇠qg18R⇠qg3112qg3113 qg4⇠⇠⇠quni11∇⇠quni1119quni11∇⇠quni1119 quni11∇⇠quni1119quni11∇⇠quni1119quni112⊘quni1119WW WWquni1119qg1818qg4⇠⇠∇+λqg2∇∇8qg314∇qg1818qg4⇠⇠∇
qg2134ℎqg314R
=−qg3R33qg18∇4qg3141
qg18∇Rqg3141qg314R 
qg4⇠⇠⇠Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 qg4⇠⇠∇qg3141qg3141
+qg3R33qg18∇Rqg3141qg314R 
qg3141
+λ =0 
where hr denotes the r-th component of h. If hr = 0, the partial derivative is positive. Hence the  91 
product of hr and the partial derivative is always zero for a solution of (5), i.e. for r = 1 … R: 92 
qg3R33qg18∇4qg3141
qg18∇Rqg3141qg314R ℎqg314R
qg4⇠⇠⇠Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 qg4⇠⇠∇qg3141qg3141
−ℎqg314Rqg343∇qg3R33qg18∇Rqg3141qg314R 
qg3141
+λqg3441=0 (6)
Since W-columns with all-zeros do not contribute to Z, it can be assumed that column sums of W 93 
are non-zero, so the above can be recast as: 94 
ℎqg314R
qg4⇠⇠⇠WW WWqg314∇qq qqqg4⇠⇠∇qg314R
qg4⇠⇠⇠WW WWqg314∇11 11qg4⇠⇠∇qg314R+λ −ℎqg314R=0 
where qn = vn / (Wh )n. To facilitate the derivations below, the following notations are introduced:  95 
qg18R3qg314R=
qg4⇠⇠⇠WW WWqg314∇qq qqqg4⇠⇠∇qg314R
qg4⇠⇠⇠WW WWqg314∇11 11qg4⇠⇠∇qg314R+λ −1 (7)
which are functions of h via q. The KKT conditions are hence recast as [20] 96 
quni1119qg18R3qg314Rquni1119ℎqg314R=0quni1119quni1119quni1119quni1119quni1119quni11⇠⇠quni11⇠Fquni11∇2quni1119quni11∇2quni1119=quni11191quni1119quni212⇠quni1119R (8)
Finally, summing (6) over r yields 97 
qg3R33qg4⇠⇠⇠qg4⇠⇠⇠qg1∇9∇qg314∇qg2∇∇8qg4⇠⇠∇qg314R+λqg4⇠⇠∇ℎqg314R
qg314R
=qg3R33qg18∇4qg3141
qg3141
 (9)
which is satisfied for any guess h by renormalizing: 98 
ℎqg314R⇠ℎqg314R
quni11∇⇠quni11∇⇠ quni11∇⇠quni11∇⇠qg314∇11 11
qg1818qg314∇qg4⇠⇠⇠qg1∇9∇qg314∇qg2∇∇8+λqg4⇠⇠∇ (10)
2.1  Multiplicative updates 99 
For the more generic class of Bregman divergences, it was shown in a.o. [20] that multipli-100 
cative updates (MU) are non-decreasing at each upda te of W and H. For KL-divergence, MU 101 
are identical to a fixed point update of (6), i.e. 102 
ℎqg314R⇠ℎqg314Rqg4⇠⇠⇠1+qg18R3qg314Rqg4⇠⇠∇=ℎqg314R
qg4⇠⇠⇠WW WWqg314∇qq qqqg4⇠⇠∇qg314R
qg4⇠⇠⇠WW WWqg314∇11 11qg4⇠⇠∇qg314R+λ (11)
Update (11) has two fixed points: hr = 0 and ar = 0. In the former case, the KKT conditions imply 103 
that ar is negative. 104 
 105 
2.2  Newton updates 106 
To find the stationary points of (2), R equations (8) need to be solved for h.  In general, let g(h) be 107 
an R-dimensional vector function of an R-dimensional variable h. Newton’s update then states: 108 
qg1818⇠qg1818−quni1119qg4⇠⇠⇠∇qg181∇qg4⇠⇠∇qg28∇9qg28⇠9qg181∇qg4⇠⇠⇠qg1818qg4⇠⇠∇quni1119quni1119quni1119wquni11⇠9tquni11⇠8 quni1119quni1119quni1119qg4⇠⇠⇠∇qg181∇qg4⇠⇠∇qg314Rqg3139=∂qg18R9qg314Rqg4⇠⇠⇠qg1818qg4⇠⇠∇
∂ℎqg3139
 (12)
Applied to equations (8): 109 
qg4⇠⇠⇠∇qg181∇qg4⇠⇠∇qg314Rqg3139=qg18R3qg3139qg2112qg314Rqg3139− ℎqg314R
qg4⇠⇠⇠WW WWqg314∇11 11qg4⇠⇠∇qg314R+λqg3R33qg18∇4qg313⇠qg18∇Rqg313⇠qg314R qg18∇Rqg313⇠qg3139
qg4⇠⇠⇠qg1∇9∇qg1818 qg4⇠⇠∇qg313⇠
qg28∇1
qg313⇠
 (13)
where δ rl  is Kronecker’s delta. To avoid the matrix inversio n in update (12), the last term in 110 
equation (13) is diagonalized, which is equivalent to solving the r-th equation in (8) for hr with all 111 
other components fixed. With 112 
qg18R4qg314R= 1
qg4⇠⇠⇠WW WWqg314∇11 11qg4⇠⇠∇qg314R+λqg3R33qg18∇4qg3141
qg18∇Rqg3141qg314R qg28∇1
qg4⇠⇠⇠Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 Wquni1119quni11⇠8 qg4⇠⇠∇qg3141qg28∇1
qg3141
 (14)
which is always positive, an element-wise Newton update for h is obtained: 113 
ℎqg314R⇠ℎqg314R
ℎqg314Rqg18R4qg314R
ℎqg314Rqg18R4qg314R−qg18R3qg314R
 (15)
Notice that this update does not automatically sati sfy (9), so updates should be followed by a 114 
renormalization (10). One needs to pay attention to  the fact that Newton updates will attract 115 
towards both local minima and local maxima. Like fo r the EM-update, hr = 0 and ar = 0 are the 116 
only fixed points of update (15), which are now sho wn to be locally  stable. In case the optimizer is 117 
at hr = 0, ar is negative by the KKT conditions, and update (15)  will indeed decrease hr. In a 118 
sufficiently small neighborhood of a point where th e gradient vanishes, i.e. ar = 0, update (15) will 119 
increase (decrease) hr if and only if (11) increases (decreases) its esti mate. Since if (11) never 120 
increases the cost, update (15) attracts to a minimum. 121 
However, this only guarantees local convergence for  per-element updates and Newton methods 122 
are known to suffer from potentially small converge nce regions. This also applies to update (15), 123 
which can indeed result in limit cycles in some cas es. In the next subsections, two measures are 124 
taken to respectively increase the convergence region and to make the update non-increasing. 125 
2.3  Step size li mitation 126 
When ar is positive, update (15) may not be well-behaved i n the sense that its denominator can 127 
become negative or zero. To respect nonnegativity a nd to avoid the singularity, it is bounded 128 
below by a function with the same local behavior around zero:  129 
ℎqg314Rqg18R4qg314R
ℎqg314Rqg18R4qg314R−qg18R3qg314R
= 1
1− qg18R3qg314R
ℎqg314Rqg18R4qg314R
≥1+ qg18R3qg314R
ℎqg314Rqg18R4qg314R
 
Hence, if ar ≥ 0, the following update is used: 
(16)
ℎqg314R⇠ℎqg314Rqg343⇠1+ qg18R3qg314R
ℎqg314Rqg18R4qg314R
qg3441=ℎqg314R+qg18R3qg314R
qg18R4qg314R
 (17)
Finally, step sizes are further limited by flooring  resp. ceiling the multiplicative gain applied to hr 130 
in update (15) and (17) (see Algorithm 1, steps 11 and 24 for details). 131 
2.4  Non-increase of the cost 132 
Despite the measures taken in section 2.3, the dive rgence can still increase under the Newton 133 
update. A very safe option is to compute the EM upd ate additionally and compare the cost 134 
function value for both updates. If the EM update i s be better, the Newton update is rejected and 135 
the EM update is taken instead. This will guarantee  non-increase of the cost function. The 136 
computational cost of this operation is dominated b y evaluating the KL-divergence, not in 137 
computing the update itself. 138 
3 The Diagonalized Newton Algorithm for KLD-NMF 139 
In Algorithm 1, the arguments given above are joined to form the Diagonalized Newton Algorithm 140 
(DNA) for NMF with Kullback-Leibler divergence cost . Matlab TM  code is available from 141 
www.esat.kuleuven.be/psi/spraak/downloads  both for the case when V is sparse or dense. 142 
 143 
Algorithm 1: pseudocode for the DNA KLD-NMF algorithm. ⊘ and ⊙ are element-wise division 144 
and multiplication respectively and [x] ε = max(x, ε). Steps not labelsed with “MU” is the additional 145 
code required for DNA. 146 
Input: data V , initial guess for W  and H , regularization weights ρ  and λ . 147 
MU - Step 1: divide the r-th column of W  by quni2211qg18∇Rqg3141qg314R qg3141+λ. Multiply the r-th row of H  by the same number. 148 
MU - Step 2: Z =  WH  149 
MU - Step 3: qg1∇91=qg1∇9⇠⊘qg1811 150 
Repeat until convergence 151 
MU - Step 4: precompute qg1∇9∇⨀qg1∇9∇ 152 
MU - Step 5: qg1∇∇R=qg1∇9∇qg314∇qg1∇91−1 153 
MU - Step 6: qg1∇82qg3114qg3122 =qg1∇82+qg1∇∇R⨀qg1∇82 154 
MU - Step 7: qg1811qg3114qg3122 =qg1∇9∇qg1∇82qg3114qg3122  155 
MU - Step 8: qg1∇91qg21⇠9qg21∇∇ =qg1∇9⇠⊘qg1811qg3114qg3122  156 
MU - Step 9: qg1814qg3114qg3122 =qg2∇∇8qg314∇qg343Rqg1∇9⇠⨀lquni11⇠Fquni11⇠∇quni1119qg4⇠⇠⇠qg1∇91qg3114qg3122 qg4⇠⇠∇qg3439 157 
Step 10: qg1∇91qg33⇠R=qg1∇9⇠⊘qg4⇠⇠⇠qg2182⨀qg2182qg4⇠⇠∇;quni1119qg1∇∇⇠=qg4⇠⇠⇠qg1∇9∇⨀qg1∇9∇qg4⇠⇠∇qg2212qg1∇91qg33⇠R  158 
Step 11:  qg1∇82qg311Rqg311Rqg3112 =quni1119qg1∇82⊙qg4⇠∇1qg1∇∇⇠⊘qg4⇠⇠⇠qg1∇∇⇠−qg1∇∇Rqg4⇠⇠∇qg4⇠∇1qg2239for the entries for which A < 0 159 
qg1∇82qg311Rqg311Rqg3112 =quni1119qg1∇82+qg1813qg1819qg1814 qg4⇠⇠⇠qg1∇∇R⊘qg4⇠⇠⇠qg1∇82⊙qg1∇∇⇠qg4⇠⇠∇quni112⊘qg2119qg1∇82qg4⇠⇠∇ for the entries for which A ≥  0 160 
  multiply t-th column of H DNA  with the t-th entry of qg4⇠⇠⇠qg2∇∇8qg314∇qg1∇9⇠qg4⇠⇠∇⊘qg4⇠⇠⇠qg2∇∇8qg314∇qg1∇82qg311Rqg311Rqg3112 qg4⇠⇠∇ 161 
Step 12: qg1811qg311Rqg311Rqg3112 =qg1∇9∇qg1∇82qg311Rqg311Rqg3112  162 
Step 13: qg1∇91qg21⇠1qg21∇1qg21R∇ =qg1∇9⇠⊘qg1811qg21⇠1qg21∇1qg21R∇  163 
Step 14: qg1814qg311Rqg311Rqg3112 =qg2∇∇8qg314∇qg343Rqg1∇9⇠⨀lquni11⇠Fquni11⇠∇quni1119qg4⇠⇠⇠qg1∇91qg311Rqg311Rqg3112 qg4⇠⇠∇qg3439 164 
Step 15:  copy H , Z  and Q  from: 165 
H DNA , Z DNA  and Q DNA  for the columns for which d DNA  < d MU  166 
H EM , Z EM  and Q EM  for the columns for which d DNA  ≥ d MU 167 
MU - Step 16: divide (multiply) the r-th row (column) of H  (W ) by quni2211ℎqg3141qg314∇ qg314∇+ρ.  168 
Step 17: precompute qg1∇82⨀qg1∇82  169 
MU - Step 18: qg1∇∇R=qg1∇91qg1∇82qg2212−1 170 
MU - Step 19: qg1∇9∇qg3114qg3122 =qg1∇9∇+qg1∇∇R⨀qg1∇9∇ 171 
MU - Step 20: qg1811qg3114qg3122 =qg1∇9∇qg3114qg3122 qg1∇82 172 
MU - Step 21: qg1∇91qg21⇠9qg21∇∇ =qg1∇9⇠⊘qg1811qg3114qg3122  173 
MU - Step 22: qg1814qg3114qg3122 =qg343Rqg1∇9⇠⨀lquni11⇠Fquni11⇠∇quni1119qg4⇠⇠⇠qg1∇91qg3114qg3122 qg4⇠⇠∇qg3439qg2∇∇8 174 
Step 23: qg1∇91qg33⇠R=qg1∇9⇠⊘qg4⇠⇠⇠qg2182⨀qg2182qg4⇠⇠∇;quni1119qg1∇∇⇠=qg1∇91qg33⇠Rqg4⇠⇠⇠qg1∇82⨀qg1∇82qg4⇠⇠∇qg2212  175 
Step 24:  qg1∇9∇qg311Rqg311Rqg3112 =quni1119qg1∇9∇⊙qg4⇠∇1qg1∇∇⇠⊘qg4⇠⇠⇠qg1∇∇⇠−qg1∇∇Rqg4⇠⇠∇qg4⇠∇1qg2239for the entries for which A < 0 176 
qg1∇9∇qg311Rqg311Rqg3112 =quni1119qg1∇9∇+qg1813qg1819qg1814 qg4⇠⇠⇠qg1∇∇R⊘qg4⇠⇠⇠qg1∇82⊙qg1∇∇⇠qg4⇠⇠∇quni112⊘qg2119qg1∇9∇qg4⇠⇠∇ for the entries for which A ≥  0 177 
  multiply the n -th row of W DNA  with the n -th entry of qg4⇠⇠⇠qg1∇9⇠qg2∇∇8 qg4⇠⇠∇⊘qg4⇠⇠⇠qg1∇9∇qg311Rqg311Rqg3112 qg2∇∇8qg4⇠⇠∇ 178 
Step 25: qg1811qg311Rqg311Rqg3112 =qg1∇9∇qg311Rqg311Rqg3112 qg1∇82 179 
Step 26: qg1∇91qg21⇠1qg21∇1qg21R∇ =qg1∇9⇠⊘qg1811qg21⇠1qg21∇1qg21R∇  180 
Step 27: qg1814qg311Rqg311Rqg3112 =qg343Rqg1∇9⇠⨀lquni11⇠Fquni11⇠∇quni1119qg4⇠⇠⇠qg1∇91qg311Rqg311Rqg3112 qg4⇠⇠∇qg3439qg2∇∇8 181 
Step 28:  copy W , Z  and Q  from: 182 
W DNA , Z DNA  and Q DNA  for the rows for which d DNA  < d MU  183 
W EM , Z EM  and Q EM  for the rows for which d DNA  ≥ d MU  184 
MU - Step 29: divide(multiply) the r-th column (row) of W  (H ) by quni2211qg18∇Rqg3141qg314R qg3141+λ.  185 
Notice that step 9, 14 22 and 27 require some care for the zeros in V, which should not contribute 186 
to the cost. In terms of complexity, the most expen sive steps are the computation of A, B, ZMU  and 187 
ZDNA , which require O( NRT ) operations. All other steps require O( NR ), O( RT ) or O( NR ) 188 
operations. Hence, it is expected that a DNA iterat ion is about twice as slow as MU iteration. On 189 
modern hardware, parallelization may however distor t this picture and hence experimental 190 
verification is requied. 191 
4 Experiments 192 
DNA and MU are run on several publicly available 1 data sets. In all cases, W is initialized with a 193 
random matrix with uniform distribution, normalized  column-wise. Then H is initialized as WtV 194 
and one MU iteration is performed. The same initial  values are used for both algorithms. Sparsity 195 
is not included in this study, so ρ = λ = 0. The algorithm parameters are set to ε = 0.01 and α =4. 196 
CPU timing measurements are obtained on a quad-core  AMD TM  Opteron 8356 processor running 197 
the MATLAB TM  code available at www.esat.kuleuven.be/psi/spraak/downloads  which uses the 198 
built-in parallelization capability. Timing measure ments on the graphical processing unit (GPU) 199 
are obtained on a TESLA C2070 running MATLAB and Accelereyes Jacket v2.3. 200 
4.1  Dense data matrices 201 
The first dataset considered is a set of 400 fronta l face greyscale 64 × 64 images of 40 people 202 
showing 10 different expressions.  The resulting 40 96 × 165 dense matrix is decomposed with 203 
factors of a common dimension R of 10, 20, 40 and 80. Figure 1 shows the KL diverg ence as a 204 
function of iteration number and CPU time as measur ed on the CPU. The superiority of DNA is 205 
obvious: for instance, at R = 40, DNA reaches the same divergence after 33 ite rations as MU 206 
obtains after 500 iterations. This implies a speed- up of a factor 15 in terms of iterations or 6.3 in 207 
terms of CPU time.  208 
 209 
Figure 1: convergence of DNA and MU on the ORL image dataset as a function of the number of 210 
iterations (left) and CPU time (right) for different ranks R. 211 
The second test case is the CMU PIE dataset which c onsists of 11554 greyscale images of 32 × 32 212 
pixels showing human faces under different illumina tion conditions and poses. The data are 213 
shaped to a dense 1024 × 11554 matrix and a decomposition of rank R = 10, 20, 40 and 80 are 214 
attempted with the MU and DNA algorithms.  As obser ved in Figure 2, the proposed DNA still 215 
outperforms MU, but by a smaller margin. 216 
                                                           
1 www.cad.zju.edu.cn/home/dengcai/Data/data.html 
0 100 200 300 400 5000
1
2
3
4
5
6
7 x 10
6
iteration
KLD
 
 
R=10
R=20
R=40
R=80
MU
DNA
0 10 20 30 40 500
1
2
3
4
5
6
7 x 10
6
time (s)
KLD
 
 
R=10
R=20
R=40
R=80
MU
DNA
 217 
Figure 2: convergence of DNA and MU on the CMU PIE image dataset as a function of the 218 
number of iterations (left) and CPU time (right). 219 
An overview of the time required for a single itera tion on both data sets is given in Table 1. For 220 
MU, the first row lists the time if the KL divergen ce is not computed as this is not required if the 221 
number of iterations is fixed in advance instead of  stopping the algorithm based on a decrease in 222 
KLD. The table shows that the computational cost of  MU can be reduced by about a third by not 223 
computing KLD. Compared to MU with cost calculation , DNA requires typically about 2.5 to 3 224 
times more time per iteration on the CPU. On the GPU, the ratio is rather 2 to 2.5.   225 
4.2  Sparse data matrices 226 
The third matrix considered originates from the NIS T Topic Detection and Tracking Corpus 227 
(TDT2). For 10212 documents (columns of V), the frequency of 36771 terms (rows of V) was 228 
counted leading to a sparse 36771 × 10212  matrix with only 0.35% non-zero entries. The  fourth 229 
matrix originates from the Newsgroup corpus results  in a 61188 × 18774 sparse frequency matrix 230 
with 0.2% non-zeros. Both for MU and DNA a MATLAB i mplementation using the sparse matrix 231 
class was made. In this case, an iteration of DNA i s twice as slow a MU iteration. Again, the 232 
convergence of both algorithms is shown in Figure 3 . In this case, DNA is only marginally faster 233 
than MU in terms of CPU time. 234 
 235 
Table 1: time per iteration in milliseconds as measured on the CPU and GPU implementations for 236 
different ranks ( R) and dense matrices (ORL/PIE).   237 
dataset 
R
processor CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU 
MU without cost 78 3.7 85 4.1 96 5.0 115 6.8 310 18 310 20 330 27 400 38 
MU with cost 114 6.4 118 7.0 130 7.9 161 9.4 480 35 490 37 510 44 580 55 
DNA 269 15.5 280 15.9 319 17.9 425 23.6 1180 71 1430 76 1720 95 19 60 125 
ORL PIE 
10 20 40 80 10 20 40 80 
238 
 239 
 240 
Figure 3: convergence of DNA and MU on the sparse TDT2 (left) and Newsgroup (right) data. 241 
5 Conclusions 242 
The DNA algorithm is based on Newton’s method for s olving the stationarity conditions of the 243 
constrained optimization problem implied by NMF. Th is paper only addresses the Kullback-244 
Leibler divergence as a cost function. To avoid mat rix inversion, a diagonal approximation is 245 
made, resulting in element-wise updates. Experiment al verification on publicly available matrices 246 
with a CPU and GPU MATLAB implementation for dense data matrices and a CPU MATLAB 247 
implementation for sparse data matrices show that, depending on the case and matrix sizes, DNA 248 
iterations are 2 to 3 times slower than MU iteratio ns. In most cases, the diagonal approximation is 249 
good enough such that faster convergence is observed and a net gain results. 250 
Since Newton updates can in general not ensure mono tonic decrease of the cost function, the step 251 
size was controlled with a brute force strategy of falling back to MU in case the cost is increased. 252 
More refined step damping methods could speed up DN A by avoiding evaluations of the cost 253 
function, which is next on the research agenda. 254 
A c kno w l e dg me nt s  255 
This work is supported by IWT-SBO project 100049 (A LADIN) and by KU Leuven research 256 
grant OT/09/028(VASI).   257 
0 20 40 60 80 1004
4.5
5
5.5
6
6.5
7 x 10
6
iteration
KLD
 
 
R=10
R=20
R=40
R=80
MU
DNA
0 20 40 60 80 100
8
8.5
9
9.5
10
10.5
11 x 10
6
iteration
KLD
 
 
R=10
R=20
R=40
R=80
MU
DNA
R e f e re nc e s  258 
[1]  D. Lee, and H. Seung, “Algorithms for non-negative matrix factorization,” Advances in Neural 259 
Information Processing Systems , vol. 13, pp. 556–562, 2001. 260 
[2]  B. Raj, R. Singh and P. Smaragdis, “Recognizing Spe ech from Simultaneous Speakers”, in proceedings 261 
of Eurospeech , pp. 3317-3320, Lisbon, Portugal, September 2005 262 
[3]  P. Smaragdis, “Convolutive Speech Bases and Their A pplication to Supervised Speaker Separation,” 263 
IEEE Transactions on Audio, Speech and Language Processing , vol. 15, pp. 1-12, January 2007  264 
[4]  P. D. O'Grady and B. A. Pearlmutter, “Discovering S peech Phones Using Convolutive Non-negative 265 
Matrix Factorisation with a Sparseness Constraint.,” Neurocomputing , vol. 72, no. 1-3, pp. 88-101, December 266 
2008, ISSN 0925-2312. 267 
[5]  M. Van Segbroeck and H. Van hamme, “Unsupervised le arning of time-frequency patches as a noise-268 
robust representation of speech,” Speech Communication , volume 51, no. 11, pp. 1124-1138, November 269 
2009. 270 
[6]  H. Van hamme, “HAC-models: a Novel Approach to Cont inuous Speech Recognition,” In Proc. 271 
International Conference on Spoken Language Process ing, pp. 2554-2557, Brisbane, Australia, September 272 
2008. 273 
[7]  X. Chen, L. Gu, S. Z. Li and H.-J. Zhang, “Learning  representative local features for face detection,”  in 274 
proceedings of the IEEE Computer Society Conference  on Computer Vision and Pattern Recognition , pp. 275 
1126-1131, Kauai, HI, USA, December 2001. 276 
[8]  D. Kim, S. Sra and I. S. Dhillon, “Fast Projection- Based Methods for the Least Squares Nonnegative 277 
Matrix Approximation Problem,” Statistical Analy Data Mining , vol. 1, 2008 278 
[9]  C.-J. Lin, “Projected gradient methods for non-nega tive matrix factorization,” Neural Computation , vol. 279 
19, pp. 2756-2779, 2007 280 
[10]  R. Zdunek, A. H. Phan and A. Cichocki, “Damped Newt on Iterations for Nonnegative Matrix 281 
Factorization,” Australian Journal of Intelligent Information Processing Systems , 12(1), pp. 16-22, 2010 282 
[11]  Y. Zheng, and Q. Zhang, “Damped Newton based Iterat ive Non-negative Matrix Factorization for 283 
Intelligent Wood Defects Detection,” Journal of software , vol. 5, no. 8, pp. 899-906, August 2010. 284 
[12]  P. Gong, and C. Zhang, “Efficient Nonnegative Matri x Factorization via projected Newton method”, 285 
Pattern Recognition , vol. 45, no. 9, pp. 3557-3565, September 2012. 286 
[13]  S. Bellavia, M. Macconi, and B. Morini, “An interio r point Newton-like method for nonnegative least-287 
squares problems with degenerate solution,” Numerical Linear Algebra with Applications , vol. 13, no. 10, pp. 288 
825-846, December 2006. 289 
[14]  R. Zdunek and A. Cichocki, “Non-Negative Matrix Fac torization with Quasi-Newton Optimization,”  290 
Lecture Notes in Computer Science, Artificial Intelligence and Soft Computing  4029, pp. 870-879, 2006 291 
[15]  R. Zdunek and A. Cichocki, "Nonnegative Matrix Fact orization with Constrained Second-Order 292 
Optimization", Signal Processing , vol. 87, pp. 1904-1916, 2007 293 
[16]  G. Landi and E. Loli Piccolomini, “A projected Newt on-CG method for nonnegative astronomical image 294 
deblurring,” Numerical Algorithms , no. 48, pp. 279–300, 2008 295 
[17]  C.-J. Hsieh and I. S. Dhillon, “Fast Coordinate Des cent Methods with Variable Selection for Non-296 
negative Matrix Factorization,” in proceedings of the 17th ACM SIGKDD International  Conference on 297 
Knowledge Discovery & Data Mining  (KDD), San Diego, CA, USA, August 2011 298 
[18]  L. Li, G. Lebanon and H. Park, “Fast Bregman Diverg ence NMF using Taylor Expansion and 299 
Coordinate Descent,” in proceedings of the 18 th  ACM SIGKDD Conference on Knowledge Discovery and 300 
Data Mining,  2012 301 
[19]  A. Cichocki, S. Cruces, and S.-I. Amari, “Generaliz ed Alpha-Beta Divergences and Their Application to 302 
Robust Nonnegative Matrix Factorization,” Entropy , vol. 13, pp. 134-170, 2011; doi:10.3390/e13010134 303 
[20]  I. S. Dhillon and S. Sra, “Generalized Nonnegative Matrix Approximations with Bregman Divergences,” 304 
Neural Information Proc. Systems , pp. 283-290, 2005 305 


Peer review:
[Review]:
I would like to thank the reviewers for their investment of time and effort to formulate their valued comments. The paper was updated according to your comments. Below I address your concerns:

A common remark is the lack of comparison with state-of-the-art NMF solvers for Kullback-Leibler divergence (KLD). I compared the performance of the diagonalized Newton algorithm (DNA) with the wide-spread multiplicative updates (MU) exactly because it is the most common baseline and almost every algorithm has been compared against it. As you suggested, I did run comparison tests and I will present the results here. I need to find a method to post some figures to make the point clear. First, I compared against the Cyclic Coordinate Descent (CCD) by Hsieh & Dhillon using the software they provide on their website. I ran the synthetic 1000x500 example (rank 10). The KLD as a function of iteration number for DNA and CCD are very close (I did not find a way to post a plot on this forum). However, in terms of CPU (ran on the machine I mention in the paper) DNA is a lot faster with about 200ms per iteration for CCD and about 50ms for DNA. Note that CCD is completely implemented in C++ (embedded in a mex-file) while DNA is implemented in matlab (with one routine in mex - see the download page mentioned in the paper). As for the comparison with SBCD (scalar block coordinate descent), I also ran their code on the same example, but unfortunately, one of the matrix factors is projected to an all-zero matrix in the first iteration. I have not found the cause yet.
What definitely needs investigation is that I observe CCD to be 4 times slower than DNA. Using my implementation for MU, 1200 MU iterations are actually as fast as the 100 CCD iteration. (My matlab MU implementation is 10 times faster than the one provided by Hsieh&Dhillon). For these reasons, I am not too keen on quickly including a comparison in terms of CPU time (which is really the bottom line), as implementation issues seem not so trivial. Even more so for a comparison on a GPU, where the picture could be different from the CPU for the cyclic updates in CCD. A thorough comparison on these two architectures seems like a substantial amount of future work.  But I hope the data above data convince you the present paper and public code are significant work. 

Reply to Anonymous 57f3
' it's not clear that matrix factorization is a problem for which optimization speed is a primary concern (all of the experiments in the paper terminate after only a few minutes)'

>> There are practical problems where NMF takes hours, e.g. the problems of [6], which is essentially learning a speech recognizer model from data. We are now applying NMF-based speech recognition in learning paradigms that learn from user interaction examples. In such cases, you want to wait seconds, not minutes. Also, there is an increased interest in 'large-sccale NMF problems'.

'Using a KL-divergence objective seems strange to me since there aren't any distributions involved, just matrices, whose entries, while positive, need not sum to 1 along any row or column. Are the entries of the matrices supposed to represent probabilities? '

>> Notice that the second and third term in the expression for KLD (Eq. 1) are normalization terms such that we don't require V or Z to sum to unity. This very common in the NMF literature, and was motivated in a.o. [1]. KLD is appropriate if the data follow a (mixture of) Poisson distribution. While this is realistic for counts data (like in the Newsgroup corpus), the KLD is also applied on Fourier spectra, e.g. for speaker separation or speech enhancement, with success. Imho, the relevance of KLD does not need to be motivated in a paper on algorithms, see also [18] and [20] ( numbering in the new paper).

'I understand that this is a formulation used in previous work ([1]), but it should be briefly explained. '
>> Added a sentence about the Poisson hypothesis after Eq. 1.

'You should explain the connection between your work and [17] more carefully. Exactly how is it similar/different? '
>> Reformulated. [17] (now [18]) uses a totally different motivation, but also involves the second order derivatives, like a Newton method.

'Has a diagonal Newton-type approach ever been used for the squared error objective? '
>> A reference is given now. Note however that KLD behaves substantially different.
'the smallest cost' -> 'leading to the greatest reduction in d_{KL}(V,Z)'
 'the variables required to compute' -> 'the quantities required to compute' 
>> corrected

You should avoid using two meanings of the word 'regularized' as this can lead to confusion. Maybe 'damped' would work better to refer to the modifications made to the Newton updates that prevent divergence? 
>> Yes. A lot better. Corrected.

'Have you compared to using damped/'regularized' Newton updates instead of your method of selecting the best between the Newton and MU updates? In my experience, damping, along the lines of the LM algorithm or something similar, can help a great deal. '
>> yes. I initially tried to control the damping by adding lambda*I to the Hessian, where lambda is decreased on success and increased if the KLD increases. I found it difficult to find a setting that worked well on a variety of problems. 

I would recommend using '	op' to denote matrix transposition instead of what you are doing. Section 2 needs to be reorganized. It's hard for me to follow what you are trying to say here. First, you introduce some regularization terms. Then, you derive a particular fixed-point update scheme. When you say 'Minimizing [...] is achieved by alternative updates...' surely you mean that this is just one particular way it might be done. 
>> That's indeed what I meant to say. 'is' => 'can be'

You say you are applying the KKT conditions, but your derivation is strange and you seem to skip a bunch of steps and neglect to use explicit KKT multipliers (although the result seems correct based on my independent derivation). But when you say: 'If h_r = 0, the partial derivative is positive. Hence the product of h_r and the partial derivative is always zero', I don't see how this is a correct logical implication. Rather, the product is zero for any solution satisfying complementary slackness. 
>> I meant this holds for any solution of (5). This is corrected.

And I don't understand why it is particularly important that the sum over equation (6) is zero (which is how the normalization in eqn 10 is justified). Surely this is only a (weak) necessary condition, but not a sufficient one, for a valid optimal solution. Or is there some reason why this is sufficient (if so, please state it in the paper!). 
>> A Newton update may yield a guess that does not satisfy this (weak) necessary condition. We can satisfy this condition easily with the renormalization (10), which is reflected in steps 16 and 29.

I don't understand how the sentence on line 122 'Therefor...' is not a valid logical implication. Did you actually mean to use the word 'therefor' here? The lower bound is, however, correct. 'floor resp. ceiling'??
>> 'Therefore' => 'To respect the nonnegativity and to avoid the singularity”

Reply to Anonymous 4322
See comparison described above.
I added more about the differences with the prior work you mention.

Reply to Anonymous 482c
See also comparison data detailed above.
You are right there is a lot of generic work on Hessian preconditioning. I refer to papers that work on damping and line search in the context of NMF ([10], [11], [12], [14] ...). Diagonalization is only related in the sense that it ensures the Hessian to be positive definite (not in general, but here is does).