2023 International Symposium on Nonlinear Theory and Its Applications
NOLTA2023, September 26-29, 2023, Catania and Online

A Novel HALS-Based Iterative Algorithm for
Randomized Nonnegative Matrix Factorization
Takao Masuda† , Kento Yamada‡ , Tsuyoshi Migita⋆ and Norikazu Takahashi⋆† Graduate School of Natural Science and Technology, Okayama University‡ Graduate School of Environmental, Life, Natural Science and Technology, Okayama University⋆ Faculty of Environmental, Life, Natural Science and Technology, Okayama University
3–1–1 Tsushima-naka, Kita-ku, Okayama 700-8530 Japan
Email: { pk396rxo@s. | pj3w250a@s. | migita-t@ | takahashi@ } okayama-u.ac.jp
Abstract— As an approach to eﬃciently perform largescale Nonnegative Matrix Factorization (NMF), randomized NMF was recently proposed. This approach reduces
the dimension of a given nonnegative matrix by multiplying it by a random matrix, and then performs matrix factorization. However, algorithms for standard NMF cannot bedirectly used in randomized NMF because linear inequality
constraints have to be considered instead of nonnegativity
ones. In this paper, we reformulate the optimization problem for randomized NMF from a diﬀerent perspective andpropose a novel iterative algorithm which is a combinationof the hierarchical alternating least squares algorithm andprojection onto the feasible region. We also prove that the
proposed algorithm has global convergence property.

An algorithm for randomized NMF was ﬁrst proposed by
Erichson et al. [5]. However, their algorithm does not exactly solve the optimization problem for randomized NMF.
In fact, it does not always return a solution that satisﬁes thelinear inequality constraints. Another algorithm for randomized NMF was proposed by the authors [6]. Although
this algorithm has been proved to be globally convergent
to the set of stationary points, it is slow because only one
variable is updated at a time.
In this paper, we reformulate the optimization problem
for randomized NMF from a diﬀerent perspective, and propose a novel iterative algorithm which is a combination of
the HALS algorithm and projection onto the feasible region. We also show through experiments using two real
datasets that the proposed algorithm is valid and eﬀective.

1. Introduction
2. Randomized NMF
Nonnegative Matrix Factorization (NMF) is the mathematical operation that decomposes a given nonnegative matrix into two low-rank nonnegative matrices [1, 2]. Since
NMF can extract latent features in the given nonnegative
matrix as a set of nonnegative vectors, it has been applied
to image processing, audio signal processing, text mining,recommender systems, network analysis, and so on. NMF
is formulated as an optimization problem with nonnegativity constraits [1]. Multiplicative Update Rule (MUR) [1]and Hierarchical Alternating Least Squares (HALS) algorithms [3] are widely used for solving such problems.
However, when the given nonnegative matrix is very large,
these algorithms require long computation time.
As an approach to eﬃciently perform large-scale NMF,
randomized NMF was recently proposed [5]. In this approach, a given nonnegative matrix is multiplied by a random matrix to obtain a new matrix of lower dimension, andthen this new matrix is decomposed into two low-rank matrices. However, algorithms for standard NMF cannot bedirectly used in randomized NMF because some of the nonnegativity constrains in the original optimization problem
are transformed into linear inequality constraints.
ORCID iDs Takao Masuda:
0009-0006-4965-1682, Kento Yamada:0009-0005-7269-0375, Tsuyoshi Migita:
0009-0006-3242-2090,
0000-0001-8222-5593
Norikazu Takahashi:

NMF is the operation that decomposes a given nonnegative matrix X ∈ R+M×N into two nonnegative matrices
W ∈ R+M×K and H ∈ R+N×K that satisfy X ≈ W HT , where
R+ denotes the set of nonnegative real numbers. NMF isformulated as an optimization problem in which the error
between X and W HT is minimized subject to the constraintthat all entries of W and H are nonnegative [1]. In the casewhere the error is measured in the Frobenius norm ∥ · ∥F ,the optimization problem for NMF is expressed as follows:
1
2
X − W HT F
2subject to W ≥ 0 M×K , H ≥ 0N×K
minimizef (W, H) =(1)

where 0I×J denotes the I × J zero matrix, and the inequality
between matrices is applied entry-wise.
Since the objective function f (W, H) of the problem (1)is not jointly convex, it is diﬃcult in general to ﬁnd an optimal solution. However, if we consider W or H as constant, the objective function becomes convex. For this reason, many existing algorithms for solving (1) such as the
MUR [1] and the HALS algorithm [3] do not update entries of W and H at the same time.
In the case where the matrix X ∈ R+M×N is very large,
the above-mentioned algorithms require long computation

This work is licensed under a Creative Commons
Attribution-NonCommercial-NoDerivatives 4.0 International.

– 704 –

time. In order to overcome this issue, Erichson et al. proposed a framework called randomized NMF [5], which
is based on approximate matrix decompositions. In theﬁrst step of the framework, one generates a random matrix Ω ∈ RN×K whose entries are drawn i.i.d. from thestandard normal distribution or the uniform distribution on[0, 1], and obtains Y = XΩ ∈ R M×K , where K is an integer
smaller than min{M, N} and R denotes the set of real numbers. The matrix Y is considered as a collection of K points
randomly chosen from the column space of X, and hence
the column space of Y is considered as a K-dimensional
subspace of that of X. Throughout this paper, we impose
the following assumption on Y.
Assumption 1 Y is a column full-rank matrix.
In the second step, one performs a QR decomposition
on Y to obtain a matrix Q ∈ R M×K and an upper triangularmatrix R ∈ RK×K that satisfy both Y = QR and QT Q = IK ,where IK denotes the K × K identity matrix. Since the setof K columns of Q is an orthonormal basis for the column
space of Y, we have X ≈ QQT X = QB, where B = QT X ∈RK×N is called a surrogate matrix of X.
In the third step, one transforms the optimization problem (1) into another with a lower dimension. Multiplying
both sides of X ≈ W HT by QT from the left and lettingW̃ = QTW, we have B ≈ W̃ HT . Since W ≈ QQTW = QW̃
follows from QQT ≈ I M , the nonnegativity constraint W ≥0 M×K is approximately equivalent to the linear inequality
constraint QW̃ ≥ 0 M×K . Therefore, one obtains the following optimization problem in W̃ ∈ RK×K and H ∈ RN×K :minimizesubject to1
2
B − W̃ HT F
2
QW̃ ≥ 0 M×K , H ≥ 0N×K .
f1 (W̃, H) =In order to avoid this ambiguity, we propose to transform(1) using the change of variable W = QW̃. In other words,
we restrict the columns of W to those obtained by linear
combination of the columns of Q. Then, the problem (1) istransformed into the optimization problem:
minimizesubject to

3. Reformulation of Randomized NMF
Erichson et al. [5] transformed (1) into (2) using thechange of variable W̃ = QTW. In this approach, however,the nonnegativity constraint W ≥ 0 M×K in (1) is not exactly
equivalent to the linear inequality constraint QW̃ ≥ 0 M×Kin (2) because QW̃ = QQTW is not exactly equal to W.

(3)

An advantage of our approach is that, for any feasible solution (W̃, H) of (3), we can immediately obtain a feasible
solution of (1), which is given by (QW̃, H).
The problems (2) and (3) have the same constraints but
diﬀerent objective functions. However, these two problems
are equivalent to each other as shown in the next theorem.
Theorem 1 The objective functions f1 in (2) and f2 in (3)diﬀer only by a constant.
Proof: Let us deﬁne an orthonormal matrix P ∈ R M×M
and a matrix V ∈ R M×K by!
W̃
P = Q Q̄ , V =0where Q̄ is an M × (M − K) matrix satisfying Q̄T Q̄ = I
and QT Q̄ = 0, that is, the set of M − K columns of Q̄ isan orthonormal basis for the orthogonal complement of thecolumn space of Q. Using the equalities PT P = PPT = I,
PV = QW̃ and V T V = W̃ TW̃, we have
1 T
2
Q X − W̃ HT F
2
!
!
2
1
1 T 2
QT
W̃
T
=
Q̄ X F
X
−
H
−
Q̄T
0
2
2
F
1 T
1 T 2
2
=
Q̄ X FP X − VHT F −
2
2
N
1X T
1 T 2
2
=
P xn − V(HT )n 2 −Q̄ X F2 n=1
2

f1 (W̃, H) =(2)

Note that the number of variables in (2) is smaller than thatin (1) because K < M. However, the feasible region of (2)is more complicated than that of (1).
In the last step, one ﬁnds an approximate local optimal
solution (W̃ ∗ , H∗ ) of (2) in some way, and then returns(QW̃ ∗ , H∗ ) as an approximate local optimal solution of (1).
Erichson et al. [5] proposed an iterative algorithm basedon the HALS algorithm for solving (2). However, theiralgorithm does not always return a solution that satisﬁes
the linear inequality constraint QW̃ ≥ 0 M×K . The authorsof the present paper proposed another iterative algorithm
for solving (2) [6]. This algorithm converges to the set of
stationary points, but is slow because only one variable isupdated at a time.

f2 (W̃, H) = 12 ∥X − QW̃ HT ∥2F
QW̃ ≥ 0 M×K , H ≥ 0N×K .

1 X Tx xn − 2xTn QW̃(HT )n2 n=1 n 1
2
− ((HT )n )TW̃ TW̃(HT )n −Q̄T X F
2
N
1X
1 T 2
2
=
xn − QW̃(HT )n 2 −Q̄ X F
2 n=1
2
N

=

1 T 2
Q̄ X F .
2
Since the second term of the last formula is a constant, thefunctions f1 and f2 diﬀer only by a constant.
□= f2 (W̃, H) −In the ﬁrst formulation of the randomized NMF [5],Erichson et al. assumed that the number of columns of Ω isequal to the target rank K. However, they also considered
the case where the number of columns of Ω is greater than
K. We thus hereafter assume that Ω is an N × L randommatrix, where L is an integer satisfying K ≤ L ≤ M − 1. In
this case, Q is an M × L matrix and W̃ is an L × K matrix.

– 705 –

Algorithm 1 Update Rule for Randomized NMF

4. Proposed Algorithm
The optimization problem (3) is equivalent to (2) asshown in Theorem 1. Also, the objective function of (2)is simpler than that of (3). Therefore, we focus our attention on (2) in this paper, and design an algorithm for ﬁndinga stationary point.
It was reported by Erichson et al. [5] that Ω based on theuniform distribution on [0, 1] gives better results than thestandard normal distribution. We thus impose the following
assumptions on the matrices Ω and Q.
Assumption 2 The entries of Ω ∈ RN×L are drawn
i.i.d. from the uniform distribution on [0, 1]. All entries
in the ﬁrst column of Q ∈ R M×L are nonnegative, and there
exists at least one nonzero entry in each row of Q.
Note that XΩ is nonnegative when Ω is nonnegative. So
we can assume without loss of generality that the ﬁrst column of Q is nonnegative. The last condition in Assumption 2 is also very mild, and usually satisﬁed.
We now present an overview of our iterative algorithm
for ﬁnding a stationary point of the optimization problem(2). Given an initial feasible solution (W̃, H) such that≥ 0, if l = 1,w̃lk (4)
= 0, if l ∈ {2, 3, . . . , L},
our algorithm updates the 2K columns of W̃ and H one by
one in the order h1 → w̃1 → h2 → w̃2 → · · · → hK →w̃K , where w̃k and hk denotes the k-th column of W̃ and H,
respectively. Note that the initial W̃ given by (4) satisﬁes
QW̃ ≥ 0 because of Assumption 2.
The update rule for (W̃, H) in our algorithm is shown in
Algorithm 1. In Step 3 of Algorithm 1, hk is updated in
the same way as the one presented in the literature [2, 4],
where δ is a positive constant and [t]+ represents the vector
obtained from t by replacing all negative entries with zeros. In Step 4, hk is normalized to guarantee the boundedness of the sequence of solutions generated. Note that thisnormalization process does not change the value of w̃k hTk
and hence the value of the objective function f1 (W̃, H). In
Steps 5–16, the system of equations:

Input: X ∈ R+M×N , Q ∈ R M×L , W̃ ∈ RL×K , H ∈ R+N×K ,Λ ∈ R+M×K , δ, τg ∈ R++ , J ∈ Z++ .
Output: Updated W̃, H and Λ.
1: Set E ← QT X − W̃ HT and k ← 1.
2: Set Sk ← E + w̃k hT
k.
3: Set hk ← [(STw̃+δhk )/(∥w̃k ∥22 + δ)]+ .
kk
4: If hk , 0N×1 then w̃k ← w̃k ∥hk ∥2 and hk ← hk /∥hk ∥2 .
Otherwise set w̃k ← 0K×1 and hk ← uk , where uk ∈ R+N
is any nonnegative unit vector.
5: Set w̃min← Sk hk .
k
6: Set Ipos ← {m ∈ {1, 2, . . . , M} | λmk > 0} and P ← 0L×L .
7: For each m ∈ Ipos , set vm ← (QT )m − P(QT )m , and thenset P ← P + vm vTm /∥vm ∥22 if vm , 0L×1 .
8: Set w̃∗k ← w̃min− Pw̃minkk .
∗9: If Qw̃k ≥ 0 M×1 then set w̃k ← w̃∗k and go to Step 12.
Otherwise set θmax ← max{θ ∈ [0, 1] | Q(w̃k + θ(w̃∗k −w̃k )) ≥ 0 M×1 }, w̃k ← w̃k + θmax (w̃∗k − w̃k ) and Inew ←{m ∈ {1, 2, . . . , M} | (Qw̃k )m = 0} \ Ipos .
10: For each m ∈ Inew , set vm ← (QT )m − P(QT )m and thenset P ← P + vm vTm /∥vm ∥22 if vm , 0L×1 .
11: Set w̃∗k ← w̃min− Pw̃minand return to Step 9.
kk
12: Set Iact ← {m ∈ {1, 2, . . . , M} | (Qw̃k )m = 0}.
13: Set j ← 1.
14: Update the values of λmk (m ∈ Iact ) one by one as


T
T
T

 −w̃min − P ′′ k (Q )m′
λ
(Q
)
m
m
m
∈I
\{m}
act
k

 .
λmk ← 
∥(QT )m ∥22+

15: If j <J then set j ← j + 1 and return to Step 14.
Otherwise return to Step 6.
− QT λk ∥22 ≤ τg then go to Step 17.
16: If 12 ∥w̃k − w̃mink
17: Set E ← Sk − w̃k hT
k.
18: If k = K then return W̃, H, Λ and stop. Otherwise setk ← k + 1 and return to Step 2.
Theorem 2 Suppose that τg is suﬃciently small. Then, for
any initial solution (W̃ (0) , H(0) ), the sequence of solutions{(W̃ (t) , H(t) )}∞t=0 generated by Algorithm 1 has at least one
convergent subsequence, and the limit of any convergent
subsequence is a stationary point of (2).

w̃k − Sk hk − QT λk = 0L×1 ,5. Experiments

λk ⊙ (Qw̃k ) = 0 M×1 ,λk ≥ 0 M×1 ,which is the Karush-Kuhn-Tucker condition for the optimization problem:
minimizesubject to12T
2 ∥w̃k ∥2 − (Sk hk ) w̃k
Qw̃k ≥ 0 M×K ,(5)

is solved approximately for w̃k and λk .
The global convergence of the proposed algorithm isshown as follows. The proof is omitted due to lack of space.

In order to evaluate the validity and eﬀectiveness of the
proposed algorithm, the authors conducted experiments using two image datasets COIL-20 [7] and Olivetti* . From
the COIL-20 dataset, a 1024 × 1440 matrix X with values between 0 and 1 was constructed. From the Olivetti
dataset, a 4096 × 400 matrix X with values between 0 and
1 was constructed. For these two nonnegative matrices, thetarget rank K was set to 10 and L, the number of columns of* https://scikit-learn.org/0.19/datasets/olivetti_
faces.html– 706 –

Table 1: Computation time for the ﬁrst 50 iterations.
Algorithm
Proposed (L = 10)
Proposed (L = 20)
Proposed (L = 30)
Proposed (L = 40)
Proposed (L = 50)
HALS [4]proposed (L=10)
proposed (L=20)
proposed (L=30)
proposed (L=40)
proposed (L=50)
HALS

objective function value

105

COIL-20 (sec)
3.65
6.75
7.56
8.04
10.84
12.16

Olivetti (sec)
3.65
6.29
8.60
9.86
14.73
13.40

6. Conclusions
0

10

20
30
number of iterations

40

50

We have reformulated randomized NMF using thechange of variable W = QW̃, and proposed a new iterativealgorithm for solving the optimization problem of randomized NMF. We have also shown through experiments using
two real datasets that the proposed algorithm is eﬀective
in terms of computation time when L is small. However,when L is large, the proposed algorithm has no signiﬁcantadvantage over the HALS algorithm. Developing faster algorithms for randomized NMF is a future challenge.

objective function value(a)proposed (L=10)
proposed (L=20)
proposed (L=30)
proposed (L=40)
proposed (L=50)
HALS

105

Acknowledgments
This work was supported by JSPS KAKENHI Grant
Number JP21H03510.

104

0

10

20
30
number of iterations

40

50

References(b)
Figure 1: Evolution of the value of f (W, H): (a) COIL-20
and (b) Olivetti.

[1] D. D. Lee and H. S. Seung, “Learning the parts ofobjects by non-negative matrix factorization,” Nature,
vol. 401, pp. 788–791, 1999.
[2] N. Gillis, Nonnegative Matrix Factorization, SIAM,
2020.

Ω, was set to 10, 20, 30, 40 and 50. The random matrix Ω
is generated using pseudo-random numbers from the uniform distribution on [0, 1]. The (l, k)-th entry of W̃ was
initialized to 1 if l = 1 and 0 otherwise so that the condition (4) is satisﬁed. H was initialized using pseudo-random
numbers from the uniform distribution on [0, 0.1], and Λ
was initialized to 0 M×K . The constants were set as follows:δ = 10−12 , τg = 10−2 and J = 20. The proposed algorithm
and the HALS algorithm [4] were implemented in Python
3.9.7 and executed on a computer with Intel Core i5-9500
CPU, 8GB RAM and Windows 10 Pro.
Figure 1 shows the evolution of the objective functionf (W, H) of (1) for the ﬁrst 50 iterations. Since W is restricted to those given by W = QW̃ in the proposed algorithm, the value of f (W, H) is higher than the HALS algorithm. It is easily seen from the ﬁgures that the valueof f (W, H) decreases monotonically in all cases, and thatthe larger the value of L, the smaller the value of f (W, H).
Table 1 shows the computation time of the proposed and
HALS algorithms for the ﬁrst 50 iterations. The proposed algorithm requires shorter computation time than the
HALS algorithm when L ∈ {10, 20, 30, 40}, but there is no
signiﬁcant diﬀerence when L = 50.

[3] A. Cichocki, R. Zdunek and S. I. Amari, “Hierarchical ALS algorithms for nonnegative matrix and 3D
tensor factorization,” Proc. ICA2007, pp. 169–176,
2007.
[4] T. Sano, T. Migita and N. Takahashi, “A novel update
rule of HALS algorithm for nonnegative matrix factorization and Zangwill’s global convergence,” Journal of Global Optimization, vol. 84, no. 3, pp.755–
781, 2022.
[5] N. B. Erichson, A. Mendible, S. Wihiborn and J.
N. Kutz, “Randomized nonnegative matrix factorization,” Pattern Recognition Letters, vol. 104, pp. 1–7,
2018.
[6] T. Masuda, T. Migita and N. Takahashi, “An algorithm for randomized nonnegative matrix factorization and its global convergence,” Proc. IEEE
SSCI2021, pp.1–7, 2021.
[7] S. Nene, S. Nayar and H. Murase, “Columbia object
image library (COIL-20),” Technical Report CUCS005-96, 1996.

– 707 –