arXiv:1612.04112v5 [] 1 Oct 2017

Upper Bound of Bayesian Generalization Error in
Non-Negative Matrix Factorization
Naoki HAYASHI, Sumio WATANABE
October 3, 2017
Abstract
Non-negative matrix factorization ( NMF ) is a new knowledge discovery method that is used for text mining, signal processing, bioinformatics,
and consumer analysis. However, its basic property as a learning machine
is not yet clarified, as it is not a regular statistical model, resulting that
theoretical optimization method of NMF has not yet established. In this
paper, we study the real log canonical threshold of NMF and give an upper bound of the generalization error in Bayesian learning. The results
show that the generalization error of the matrix factorization can be made
smaller than regular statistical models if Bayesian learning is applied.

1

Introduction

Recently, non-negative matrix factorization ( NMF ) [1, 2] has been applied
to text mining [3], signal processing [4, 5, 6], bioinformatics [7], and consumer
analysis [8]. Experiments has shown that a new knowledge discovery method is
derived by NMF, however, its mathematical property as a learning machine is
not yet clarified, since it is not a regular statistical model. A statistical model
is called regular if a function from a parameter to a probability density function
is one-to-one and if the likelihood function can be approximated by a Gaussian
function.
It is proved that, if a statistical model is regular and if a true distribution is
realizable by a statistical model, then the generalization error is asymptotically
equal to d/(2n), where d, n, and the generalization error are the dimension of the
parameter, the sample size, and the expected Kullback-Leibler divergence of the
true distribution and the estimated learning machine, respectively. However, the
statistical model used in NMF is not regular because the map from a parameter
to a probability density function is not injective. As a result, its generalization
error is still unknown. If there is a theoretical value of generalization error, we
can confirm whether the numerical calculation result is correct or not.
There are many non-regular statistical models in machine learning. For example, a neural network, a reduced rank regression, a normal mixture, a hidden

1

Markov model, and Boltzmann machine are such examples. From the theoretical point of view, the generalization error of a non-regular learning machine in
Bayesian learning was proved to be asymptotically equal to λ/n, where λ is a
real log canonical threshold ( RLCT ) [9, 10]. The RLCTs for several learning
machines, have been clarified. In fact, a 3-layer neural network [10], a reduced
rank regression [11], a normal mixture [12], and a hidden Markov model [13],
they are clarified by using resolution of singularities. A statistical model selection method using RLCT has also been proposed [14].
In this paper, we theoretically derive the upper bound of the RLCT of NMF,
by which we can derive the upper bound of the Bayesian generalization error
of NMF. This is the first derivation of the theoretical generalization error in
Bayesian NMF. This paper consists of five parts. In the second section, we
describe the upper bound of the RLCT in NMF in the Main Theorem. In the
third section, we mathematically prepare the proof of the Main Theorem. In
the fourth section, we prove the Main Theorem. And in the fifth section, we
describe a theoretical application of the Main Theorem to Bayesian learning.

2

Framework and Main Result

In this section, we explain the framework of Bayesian learning and introduce
the main result of this paper.

2.1

Framework of Bayesian Learning

First, we explain the general theory of Bayesian learning.
Let q(x) and p(x|θ) be probability density functions on a finite dimensional
real Euclidean space, where θ is a parameter. In learning theory, q(x) and
p(x|θ) represent a true distribution and a learning machine with θ respectively.
A probability density function on a set of parameters ϕ(θ) is called a prior.
Let X n = (X1 , X2 , ..., Xn ) be a set of random variables that are independently
subject to q(x), where n and X n are referred to as the sample size and training
data. The posterior distribution of w is defined by
p(θ|X n ) =

n
Y
1
ϕ(θ)
p(Xi |θ),
Z
i=1

where Z is the normalizing constant that is determined by the condition
1. The Bayesian predictive distribution is also defined by
Z
n
p(x|X ) = p(x|θ)p(θ|X n )dθ.

R

p(θ|X n )dθ =

The generalization error Gn is defined by the Kullback-Leibler divergence from
the true distribution q(x) and the predictive one p(x|X n ),
Z
q(x)
dx.
G = q(x) log
p(x|X n )
2

Note that G is a function of X n hence it is also a random variable. Its expected
value overall training data E[G] is called the expected generalization error. Assume there exists a parameter θ0 that satisfies q(x) = p(x|θ0 ). By the singular
learning theory [9, 10], it was proven that
 
λ
1
E[G] = + o
n
n
holds when n tends to infinity, even if the posterior distribution can not be
approximated by any normal distribution. The constant λ is the RLCT which is
an important birational invariant in algebraic geometry. From the mathematical
point of view, RLCT is characterized by the following property. We defined a
zeta function by
Z
ζ(z) = K(θ)z ϕ(θ)dθ,
(1)

where

K(θ) =

Z

q(x) log

q(x)
dx.
p(x|θ)

Then this is holomorphic in Re(z) > 0 which can be analytically continued to
a unique meromorphic function on the entire complex plane. The poles of this
extended function are all negative rational numbers. Let (−λ) be the nearest
pole to the origin; λ is then equal to the RLCT. If p(x|θ) is regular then λ = d/2;
however, it is not usually general. In this paper we show an upper bound of the
RLCT of the NMF.

2.2

Main Theorem

Second, we introduce the main result of this paper. In the followings, θ = (X, Y )
is a parameter and x = W is an observed random variable.
Let M(M, N, C) be a set of M × N matrices whose elements are in C, where
C is a subset of R. Let K be a compact subset of R≧0 = {x ∈ R|x ≧ 0}
and let K0 be a compact of subset R>0 = {x ∈ R|x > 0}. We denote that
X ∈ M(M, H, K), Y ∈ M(H, N, K) and A ∈ M(M, H0 , K0 ), B ∈ M(H0 , N, K0 )
are NMFs of AB such that they give the minimal H0 , where H ≧ H0 and
{(x, y, a, b) ∈ K 2 × K02 |xy = ab} 6= ∅. We set that k · k is Frobenius norm.
Definition 2.1 (An RLCT of NMF). Assume that the largest pole of the
function of one complex variable z,
Z
Z
z

dY kXY − ABk2
dX
ζ(z) =
M(M,H,K)

M(H,N,K)

is equal to (−λ). Then λ is said to be the RLCT of the NMF.
In this paper, we prove the following theorem.

3

Theorem 2.1 (Main Theorem). The RLCT of NMF λ satisfies the following
inequality:
1
λ ≦ [(H − H0 ) min{M, N } + H0 (M + N − 1)] .
2
If H = H0 = 1, then the above bound gives the exact value.
We prove this theorem in the next section. As an application of this theorem,
we obtain an upper bound of the Bayesian generalization error of NMF. Sometimes NMF is studied in the case when only one target matrix is decomposed,
however, in general, decomposition of a set of independent matrices should be
studied because target matrices are often obtained daily, monthly, or different
places [8]. In such cases, decomposition of a set of matrices results in statistical inference. The following theorem shows a statistical bound of Bayesian
estimation of NMF.
Theorem 2.2. Let the probability density functions of W ∈ M(M, N, K) be
q(W ) and p(W |X, Y ), which represent a true distribution and a learning machine respectively defined by


1
q(W ) ∝ exp − kW − ABk2 ,
2


1
p(W |X, Y ) ∝ exp − kW − XY k2 .
2
Also let ϕ(X, Y ) be a probability density function such that it is positive on a
compact subset of M(M, H, K) × M(H, N, K) including 0. Then, the expected
generalization error E[G] satisfies the following inequality:
 
1
1
[(H − H0 ) min{M, N } + H0 (M + N − 1)] + o
.
E[G] ≦
2n
n
If H = H0 = 1, then the above bound gives the exact value.
In this theorem, we study a case when a set of random matrices W1 , W2 , ..., Wn
are observed and the true decomposition A and B are statistically estimated.
A statistical model p(W |X, Y ) which has parameters (X, Y ) are employed for
estimation. Then the generalization error of Bayesian estimation is given by
this theorem. If Theorem2.1 is proved, Theorem2.2 is immediately derived.
Therefore, we prove Main Theorem2.1 in the followings.

3

Preparation

We need the following five lemmas in order to prove the Main Theorem.
Lemma 3.1. Put θ = (X, Y ), ϕ = 1, and K(X, Y ) = kXY k2 in equlity
(1),where X ∈ M(M, H, K) and Y ∈ M(H, N, K). The RLCT λ satisfies the
following equality :
H min{M, N }
λ=
.
2
4

H,N
Proof. We set X = (xik )M,H
i,k=1,1 , Y = (ykj )k,j=1,1 , then

2

kXY k =

H
M X
N
X
X
i=1 j=1

xik ykj

k=1

!2

.

Then, since elements of matrices is nonnegative
H
X

2
x2ik ykj

≦

k=1

H
X

xik ykj

k=1

!2

≦H

H
X

2
.
x2ik ykj

k=1

As a log canonical threshold is not changed by any constant factor and it is
order isomorphic, all we have to do is calculating an RLCT of
!
! N
H
M X
N X
H
M
X
X
X
X
2
2
2 2
ykj .
xik
xik ykj =
i=1 j=1 k=1

k=1

i=1

j=1

The RLCT λ becomes a sum of each ones about k. For each k, we consider
blow-ups of variables {xik }, {ykj } each other.

!
!
H
M
N
 X
X
X
1
1 
λ =
min
,

2
2 
i=1
j=1
k=1

=

H
X
min{M, N }

2

k=1

=

H min{M, N }
.
2

Remark 3.1. Lemma 3.1 means that the bound in Main Theorem formally
gives the exact value if H0 = 0 i.e. AB = O.
Lemma 3.2. If H0 = H = 1, the equal sign of the Main Theorem holds.
T
T
N
N
.
, A = (ai )M
We set X = (xi )M
i=1 , and B = (bj )j=1
i=1 . Y = (yj )j=1
We prove the next lemma for the proof of Lemma 3.2.
Lemma 3.3. We put fkl := xk yl − ak bl (k, l ∈ N). Then ∀M, N ∈ N≧2 = {n ∈
N|n ≧ 2}, ∃C > 0 s.t.
!
N
M
X
X
2
2
2
2
.
(2)
f1j
+ f11
fi1
+
fMN
≦C
j=2

i=2

Proof. We arbitrarily take i, j(2 ≦ i ≦ M, 2 ≦ j ≦ N ) and fix them. We
consider the case where fij = 0. The elements of the matrices are positive,

5

which causes
xi yj − ai bj =

xi
(xi−1 yj − ai−1 bj )
xi−1
ai−1 bj
+
(xi yj−1 − ai bj−1 )
xi−1 yj−1
ai b j
(xi−1 yj−1 − ai−1 bj−1 ).
−
xi−1 yj−1

Generally,
(α + β + γ)2 ≦ 2(α2 + β 2 + γ 2 ), for ∀α, β, γ ∈ R
is attained. We square both side
fij2

≦

!
2
 x 2
 a b 2
 ab
i
i−1 j
i j
2
2
2
2
f(i−1)j +
fi(j−1) +
f(i−1)(j−1)
xi−1
xi−1 yj−1
xi−1 yj−1

2

i
Because of X ∈ M(M, 1, K), Y ∈ M(1, N, K) and what K is compact, 2 xxi−1
2
2


ai−1 bj
a b
, 2 xi−1
and 2 xi−1i yjj−1
have a maximum value m(i−1)j , mi(j−1) and
yj−1
m(i−1)(j−1) ,respectively. Thus
2
2
2
.
+ m(i−1)(j−1) f(i−1)(j−1)
+ mi(j−1) fi(j−1)
fij2 ≦ m(i−1)j f(i−1)j

Besides, we set C(i, j) = max{m(i−1)j , mi(j−1) , m(i−1)(j−1) } and get
2
2
2
).
+ f(i−1)(j−1)
+ fi(j−1)
fij2 ≦ C(i, j)(f(i−1)j

(3)

We prove inequality(2) by mathematical induction of M, N ∈ N≧2 , using inequality(3).
(Step 1) In the case of (M, N ) ∈ (N≧2 )2 s.t. M + N = 4, 5.
M = N = 2 ⇒ clear. In case of M = 2 ∧ N = 3, using inequality(3), we give
2
2
2
2
)
≦ C(2, 3)(f13
+ f22
+ f12
f23

and
2
2
2
2
+ f21
+ f11
).
f22
≦ C(2, 2)(f12

We show the following inequality, using them.
 2
2
2
2
2
2
.
f23
≦ C(2, 3) f13
+ C(2, 2)(f12
+ f21
+ f11
) + f12

Thus, if we set C := C(2, 3) max{1, C(2, 2)}, the inequality(2) is attained. The
case of M = 3 ∧ N = 2 can be proven in the same way.
(Step 2) We assume that inequality(2) is attained if (M, N ) ∈ (N≧2 )2 s.t. M +
N = i + j, i + j + 1.
(Case 2-1) In case of i − k ≧ 2 ∧ j + k ≧ 2 i.e. 2 − j ≦ k ≦ i − 2.

6

Using inequality(3),

2
2
2
2
.
+ f(i−k)(j−k)
+ f(i−k+1)(j+k)
≦ C(i − k + 1, j + k + 1) f(i−k)(j+k+1)
f(i−k+1)(j+k+1)
According to the above assumption, ∃C1 , C2 , C3 > 0 s.t.
2
f(i−k)(j+k+1)

≦ C1

i−k
X

2
fi1

j+k+1
X

+

≦ C2

i−k+1
X

2
fi1

+

i−k
X

≦ C3

2
fi1

+

2
f(i−k+1)(j+k+1)

j+k
X

2
f1j

j+k
X

2
f1j

+

2
f11

!

!

+

2
f11

+

j+k
X

j=2

i=2

Using the above three inequalities,

!

j=2

i=2

2
f(i−k)(j+k)

+

2
f11

j=2

i=2

2
f(i−k+1)(j+k)

2
f1j

.

(

≦ C(i−k+1, j +k+1)

(C1 +C2 +C3 )

i−k
X

2
fi1

2
f1j

+

2
f11

j=2

i=2

!

)

2
2
.
+C1 f1j+k+1
+ C2 f(i−k+1)1

Hence we set C := C(i−k+1, j+k+1)(C1+C2+C3 ), the inequality(2) is attained
because of C1 , C2 , C3 > 0.
(Case 2-2) In case of k = 1 − j.
Using inequality(3),
2
2
2
2
).
+ f(i+j−1)1
+ f(i+j)1
≦ C(i + j, 2)(f(i+j−1)2
f(i+j)2

According to the above assumption, ∃C ′ > 0 s.t.
2
≦ C′
f(i+j−1)(2)

i+j−1
X

!

2
2
2
fi1
+ f12
+ f11
.

i=2

Using the above inequality,
(
2
fi+j2

≦ C(i + j, 2) C

′

i+j−1
X

2
fi1

+

i=2

2
f12

+

2
f11

!

+

2
f(i+j)1

+

2
f(i+j−1)1

)

.

Hence, we set C := C(i + j, 2) max{1, C ′ }, and inequality(2) is attained.
(Case 2-3) In case of k = i − 1
This can be derived in the same way as Case 2-2 and Lemma 3.3, as follows.

7

According to Lemma 3.3, the next corollary is easily verified.
Corollary 3.1. Under the same assumption as Lemma 3.3, ∃D > 0 s.t.
!
N
M
N
M X
N
M
X
X
X
X
X
2
2
2
2
2
2
2
f1j + f11 .
fi1 +
fij ≦ D
f1j + f11 ≦
fi1 +
i=2

j=2

i=2

i=1 j=1

j=2

We immediately get the next lemma using Corollary 3.1.
Lemma 3.4. When an RLCT of F is equal to one of G, it is denoted by
F ∼ G. Under the same assumption as Lemma 3.3, the following relation is
attained:
M X
N
X

fij2 ∼

M
X

2
fi1
+

2
2
f1j
+ f11
.

(4)

j=2

i=2

i=1 j=1

N
X

Using Lemma 3.4, Lemma 3.2 is proven.
N
M
Proof of Lemma 3.2. We set X = (xi )M
i=1 , Y = (yj )j=1 , A = (ai )i=1 , B =
(bj )N
j=1 , fkl = xk yl − ak bl , k, l ∈ N. On account of Lemma 3.4,

kXY − ABk2 =

M X
N
X

fij2 ∼

M
X

2
fi1
+

2
2
f1j
+ f11
.

j=2

i=2

i=1 j=1

N
X

Thus, all we have to do is calculate an RLCT of the right side. We arbitrarily
take i, j(1 ≦ i ≦ M, 1 ≦ j ≦ N, i, j ∈ N) and fix them. Let g : (xi , yj ) 7→ (ξi , fij )
be the following coordinate transformation:
ξi

=

xi

fij

=

xi yj − ai bj .

The Jacobi matrix of the above transformation is equal to
! 

∂fij
∂ξi
∂(ξi , fij )
1 yj
∂xi
∂xi
.
=
=
∂fij
∂ξi
0 xi
∂(xi , yj )
∂yj
∂yj
Because of
∂(ξi , fij )
= xi > 0,
∂(xi , yj )
g is an analytic isomorphism. Hence the zeta function of (xi yj − ai bj )2 is equal
to
ZZ
ZZ
fij2z ξi dξi dfij .
(xi yj − ai bj )2z dxi dyj =
g(K 2 )

K2

8

Owing to ξi = xi > 0, ξi does not contribute to the maximum pole of the zeta
function. Therefore we can consider only
ZZ
fij2z dξi dfij .
g(K 2 )

Since g is an analytic isomorphism and K 2 is compact, ∃c1 , c2 , c3 , c4 ∈ R (c1 <
c2 , c3 < c4 ) s.t.

g(K 2 ) = (ξi , fij ) ∈ R2>0 |c1 ≦ ξi ≦ c2 ∧ c3 ≦ fij ≦ c4
This is clearly integrable, therefore we calculate an iterated integral. Thus
Z c4
Z c2
ZZ
fij2z dfij
dξi
fij2z dξi dfij =
g(K 2 )

c1

=

c3

C
for ∃C 6= 0(constant).
2z + 1

Therefore an RLCT of fij2 = (xi yj − ai bj )2 is equal to 1/2.
According to Lemma 3.4, f11 , f21 , ..., fM1 , f12 , ...,P
f1N are independent
and
PN
2
2
2
the number of them is M + N − 1. An RLCT λ of M
i=2 fi1 +
j=2 f1j + f11
is equal to
N
M
X
M +N −1
1 X1 1
+
+ =
.
λ=
2 j=2 2 2
2
i=2
Let X ∈ M(M, H, K), Y ∈ M(H, N, K) be
X = (x1 , . . . , xH ), xk = (xik )M
i=1 ,
Y = (y1 , . . . , yH )T , yk = (ykj )N
j=1 ,
and A ∈ M(M, H0 , K0 ), B ∈ M(H0 , N, K0 ) be
A = (a1 , . . . , aH0 ), ak = (aik )M
i=1 ,
B = (b1 , . . . , bH0 )T , bk = (bkj )N
j=1 ,
respectively.
Lemma 3.5. In the case of H = H0 , the Main Theorem is attained.

9

Proof.
kXY − ABk2
=

M X
N
X

(xi1 y1j + ... + xiH yHj − ai1 b1j − aiH bHj )2

i=1 j=1

=

H
M X
N
X
X
i=1 j=1

≦

C

!2

(xik ykj − aik bkj )

k=1

M X
N X
H
X

(xik ykj − aik bkj )2 for ∃C > 0(const.)

i=1 j=1 k=1

∼

M X
N X
H
X

(xik ykj − aik bkj )2

i=1 j=1 k=1

=

H
M X
N
X
X
k=1

=

H
X

2

(xik ykj − aik bkj )

i=1 j=1

xk ykT − ak bTk

2

!

.

k=1

P
The variables in the
of the last term are independent for k. That causes,
using Lemma 3.2, an RLCT λ̃ of the term is equal to
λ̃ =


H 
X
M +N −1
M +N −1
=H
.
2
2
k=1

An RLCT is order isomorphic therefore we get
λ ≦ λ̃ = H

M +N −1
.
2

10

4

Proof of Main Theorem

Proof of Main Theorem.
kXY − ABk2
=

M X
N
X

(xi1 y1j + ... + xiH yHj − ai1 b1j − ... − aiH0 bH0 j )2

i=1 j=1

=

H
M X
N
X
X
i=1 j=1

=

H0
M X
N
X
X
i=1 j=1

≦ C

=

=

2

(xik ykj − aik bkj ) +

k=1

xik ykj
H
X

(xik ykj − aik bkj )2 +

2
x2ik ykj

H
X

2
x2ik ykj

k=H0 +1

(xik ykj − aik bkj )2

i=1 j=1

xk ykT − ak bTk

k=1

!2

k=H0 +1

k=1

H0
M X
N
X
X

H0
X

H
X

k=H0 +1

H0
M X
N
X
X

k=1

!2

(xik ykj − aik bkj ) +

H0
M X
N
X
X

i=1 j=1

aik bkj

k=1

k=1

i=1 j=1

∼

xik ykj −

k=1

H0
X

!

+

H
X

k=H0 +1



x1(H0 +1)
2

..
+ 
.
xM(H0 +1)

!

for ∃C > 0(const.)

!

M X
N
X
i=1 j=1

2
x2ik ykj

!


y(H0 +1)1
. . . x1H
..
..  
..
.
.
. 
yH1
. . . xMH


. . . y(H0 +1)N

..
..

.
.
...
yHN

Let λ̃ be an RLCT of the right-most side, λ1 be an RLCT of 1st term in
the right-most side, and λ2 be an RLCT of the 2nd one. Since variables are
independent and RLCTs are order isomorphic,
λ ≦ λ̃ = λ1 + λ2 .

(5)

Since the 1st term corresponds to the proof of Lemma 3.5 in the case of H ← H0 ,
M +N −1
.
2
In contrast, the 2nd term corresponds to Lemma 3.1 in the case of H ← H −H0 .
That causes
(H − H0 ) min{M, N }
.
λ2 =
2
Using the above two equalities for inequality(5),
λ1 = H0

λ

≦ λ̃
= λ1 + λ2
M + N − 1 (H − H0 ) min{M, N }
= H0
+
.
2
2
11

2

.

Therefore

1
[(H − H0 ) min{M, N } + H0 (M + N − 1)] .
2

λ≦

Remark 4.1. Under the same assumption of the Main Theorem, suppose
(
xik ykj − aik bkj (k ∈ {1, . . . , H0 })
k
fij =
.
xik ykj
(k ∈ {H0 + 1, . . . , H})
If fijk ≧ 0, i.e. , xik ykj − aik bkj ≧ 0

(k ∈ {1, . . . , H0 }),

1
[(H − H0 ) min{M, N } + H0 (M + N − 1)] .
2

λ=

Proof. Owing to fijk ≧ 0,
H
X

(fijk )2

≦

k=1

H
X

fijk

k=1

!2

≦H

(fijk )2

H
X

∼

k=1

Using the above relation,

(fijk )2 .

k=1

Thus,
H
X

H
X

fijk

k=1

!2

.

kXY − ABk2
=

H0
M X
N
X
X
i=1 j=1

=

∼

i=1 j=1

=

H0
X

fijk

!2

(fijk )2

k=1

H0
N
M X
X
X

2

(xik ykj − aik bkj ) +

Therefore
λ=

H
X

2
x2ik ykj

k=H0 +1

xk ykT − ak bTk

!2

!

k=1

k=1

xik ykj

k=H0 +1

k=1

H
M X
N
X
X
i=1 j=1

=

(xik ykj − aik bkj ) +

k=1

H
M X
N
X
X
i=1 j=1

H
X



x1(H0 +1)
2

..
+ 
.
xM(H0 +1)

!


y(H0 +1)1
. . . x1H
..
..  
..
.
.
. 
yH1
. . . xMH

1
[(H − H0 ) min{M, N } + H0 (M + N − 1)] .
2

12


. . . y(H0 +1)N

..
..

.
.
...
yHN

2

.

5

Discussion

In this section, we discuss the result of this paper from four points of view.

5.1

Minus Log Bayesian Marginal Likelihood

First, let us discuss the asymptotic expansion of the minus log Bayesian marginal
likelihood, which is defined by
F = − log

Z

ϕ(X, Y )

n
Y

p(Wi |X, Y )dXdY,

i=1

and is sometimes used in model selection and hyperparameter optimization. It
is well known [9] that F has an asymptotic expansion,
F =−

n
X

log q(Wi ) + λ log n + op (log n),

i=1

where λ is an RLCT. Therefore, using the results of this paper, we can also
derive the upper bound of the minus log Bayesian marginal likelihood.
There is an information criterion called ”singular BIC” or ”sBIC”. Singular
BIC uses the exact value of RLCT or its upper bound in order to approximate the minus log marginal likelihood of singular model [14]. NMF is singular
statistical model thus our theoretical result can be applied to model selection
problems in NMF by using sBIC.

5.2

Relationship with Reduced Rank Regression

Second, the tightness of the upper bound is considered. If the observed {Wi } are
taken from the same probability density function q(W ) and if the coefficients of
parameter (X, Y ) are not limited to non-negative values but can be made plus
and minus values, then the RLCT is equal to that of the reduced rank regression,
since the zeta function of NMF is equal to that of reduced rank regression. The
RLCT of reduced rank regression was clarified in [11], which is a lower bound of
NMF. The smallest inner matrix dimension of an NMF is called a nonnegative
rank [15]. Suppose C ∈ M(M, N, K). Let rank+ C be the nonnegative rank of
C. The following inequality is known [15]:
rankC ≦ rank+ C ≦ min{M, N }.
We can also prove that the RLCT of reduced rank regression is lower than the
RLCT of NMF in cases where the equality of the Main Theorem is attained,
i.e. , in case of Lemma 3.1 and Lemma 3.2. Especially in the case of Lemma
3.2, the equal sign holds. We prove it below.
Proof. Let λrRRR and λkN MF be the RLCT of reduced rank regression and the
RLCT of NMF in case of Lemma k where r is the rank of true parameter AB.
13

We only have to consider r = 0 or r = 1, respectively. Because of the Main
Theorem and Lemma 3.1 and Lemma 3.2,
λ3.1
N MF

=

λ3.2
N MF

=

H min{M, N }
,
2
M +N −1
.
2

In contrast, λrRRR is equal to the following [11] :

2(H + r)(M + N ) − (M − N )2 − (H + r)2 /8






 if(N + r ≦ M + H ∧ M + r ≦ N + H ∧ H + r ≦ M + N ∧ M + H + N + r : even)

2
2


 2(H + r)(M + N ) − (M − N ) − (H + r) + 1 /8
r
λRRR =
if(N + r ≦ M + H ∧ M + r ≦ N + H ∧ H + r ≦ M + N ∧ M + H + N + r : odd)


(HM − Hr + N r)/2 if(M + H < N + r)





(HN − Hr + M r)/2 if(N + H < M + r)



M N/2 if(M + N < H + r).
(6)
(Case 1) r = 0 i.e. in the case of Lemma 3.1.
(Case 1-1) N ≦ M + H ∧ M ≦ N + H ∧ H ≦ M + N ∧ M + H + N : even.
We assume N ≦ M , i.e. , λ3.1
N MF = HM/2. Owing to equality (6),

λ0RRR =
2H(M + N ) − (M − N )2 − H 2 /8

=
2(M H + HN + N M ) − M 2 − N 2 − H 2 /8.

λ0RRR − λ3.1
= − M 2 + H 2 + N 2 − 2(M H + HN + N M ) − 4HM /8
N MF

= − M 2 + H 2 + N 2 − 2(−M H + HN + N M ) /8


= − (−M )2 + (−H)2 + N 2 + 2{(−M )(−H) + (−H)N + N (−M )} /8
= −(N − M − H)2 /8 ≦ 0

Therefore

λ0RRR ≦ λ3.1
N MF .

If N > M , that can be derived in the same way as above.
(Case 1-2) N ≦ M + H ∧ M ≦ N + H ∧ H ≦ M + N ∧ M + H + N : odd.
We assume N ≦ M , i.e. , λ3.1
N MF = HM/2. In the same way as Case 1-1,
λ0RRR − λ3.1
N MF

=
=

1/8 − (N − M − H)2 /8
12 /8 − (M + H − N )2 /8

=

−(M + H − N + 1)(M + H − N − 1)/8.

If H = 0, N + 1 ≦ M is attained because of that M + N must be odd and
N ≦ M , hence λ0RRR ≦ λ3.1
N MF . Else, owing to H ≧ 1 and N ≦ M ,
N +1 ≦
≦
14

M +1
M + H.

We also get N − 1 ≦ N + 1 ≦ M + H thus
−(M + H − N + 1)(M + H − N − 1) ≦ 0.
Therefore,

λ0RRR ≦ λ3.1
N MF .

If N > M , that can be derived in the same way as above.
(Case 1-3) N + H < M i.e. N < N + H < M .
0
On account of N < M ,λ3.1
N MF = HM/2 = λRRR .
(Case 1-4) M + H < N i.e. M < M + H < N .
0
In the same way as Case 1-3, λ3.1
N MF = HN/2 = λRRR .
(Case 1-5) M + N < H i.e. N < M + N < H ∧ M < M + N < H.
On account of M < H ∧ N < H, M N < HN ∧ M N < HM i.e. M N <
H min{M, N }. Thus
0
λ3.1
N MF = H min{M, N }/2 > M N/2 = λRRR .

From the above, λ0RRR ≦ λ3.1
N MF .
(Case 2) r = 1 ∧ H = 1 i.e. in case of Lemma 3.2.
As for the claim λ1RRR ≦ λ3.2
N MF , we only have to prove in the case of
M + H < N + r and N + H < M + r since it contradicts in other cases. In the
first and second cases, owing to r = H = 1, it becomes M < N ∧ N < M and
contradicts. In the last(fifth) case, it is 2 ≦ M + N < 2 and conflicts. Thus
we consider the third and fourth cases:M + H < N + r and N + H < M + r,
respectively.
(Case 2-1) M + H < N + r i.e. M < N .
According to the equality (6) and Lemma 3.2,
M +N −1
= λ3.2
N MF .
2
(Case 2-2) N + H < M + r, i.e. , N < M .
That can be derived in the same way as Case 2-1.
From the above, λ1RRR = λ3.2
N MF .
λ1RRR =

That is why the RLCT of the reduced rank regression is considered as a
lower bound of the RLCT of the NMF.
H0 and r are denoted by the non-negative rank and the usual rank respectively. We compare the RLCT of reduced rank regression with the exact value and upper bound and summarize in Table1 below in Table section.
Note that there is no non-negative matrix which has larger rank than nonnegative rank. Moreover, if M ≦ 3 or N ≦ 3, then rankC = rank+ C, where
C ∈ M(M, N, K)[15]. For example, the following 4 × 4 non-negative matrix C
satisfying rankC = 3 and rank+ C = 4 is well-known [15]:


1 1 0 0
 1 0 1 0 

C=
 0 1 0 1 .
0 0 1 1
15

5.3

Robustness on Distribution

Third, we study the generalization of the result to another distribution. In
Theorem2.2, we studied a case where the matrix W is generated by normal distributions whose averages are AB and XY . Then Kullback-Leibler divergence
KL(X, Y ) satisfies KL(X, Y ) ∼ kXY − ABk2 , as is well known [11]. If W is
generated by Poisson distributions, due to the conditions of the parameters of
these distributions, elements of XY are restricted by positive elements. However, the Kullback-Leibler divergence has same RLCT as the square error if
elements of XY are positive.
Proposition 5.1. Let the probability density functions of W ∈ M(M, N, K)
be q(W ) and p(W |X, Y ), which represent a true distribution and a learning
machine respectively defined by
q(W ) ∝
p(W |X, Y ) ∝

Poi(W |AB),
Poi(W |XY ),

where Poi(W |C) is a probability density function of the Poisson distribution
with average C. Also let ϕ(X, Y ) be a probability density function such that it
is positive on a compact subset of M(M, H, K0 ) × M(H, N, K0 ). Then, the
Kullback-Leibler divergence has the same RLCT as the square error.
Proof. Let x ∈ N ∪ {0}, a > 0, b > 0. We put
p(x|a)

:=

K(a, b) :=

e−a ax
,
x!
∞
X
p(x|a)
p(x|a)log
.
p(x|b)
x=0

Using
log

p(x|a)
p(x|b)

and

e−a ax
e−b bx
− log
x!
x!
= −a + xloga + b − xlogb
a
= b − a + xlog
b

= log

∞
X

p(x|a) = 1,

x=0

K(a, b) =
=

∞
a X e−a ax−1
a
b x=1 (x − 1)!
a
b − a + alog .
b

b − a + log

16

(7)

Then, owing to
∂a K(a, b) =
∂b K(a, b) =

loga − logb,
1 − a/b,

and that a log function is monotone increasing,
∂a K(a, b) = ∂b K(a, b) = 0 ⇔ a = b.
Signs of the above partial derivations are
∂a K(a, b) > 0 ∧ ∂b K(a, b) < 0 in case of a > b
∂a K(a, b) < 0 ∧ ∂b K(a, b) > 0 in case of a < b
On account of that and smoothness, the increase or decrease and convexity
of K(a, b) is the same as those of (b − a)2 . Hence ∃c1 , c2 > 0 s.t.
c1 (b − a)2 ≦ K(a, b) ≦ c2 (b − a)2 .

(8)

i.e. K(a, b) ∼ (b − a)2 .
We assumed that elements of matrices are generated by Poisson distributions. Using inequality(8) for each element,
Kullback-Leibler divergence ∼ kXY − ABk2
,where a is an element of AB and b is an element of XY .
We also consider the case which W is generated by exponential distribution.
The elements of XY must be restricted by positive ones as in the case of Poisson.
However, also in this case, the Kullback-Leibler divergence has same RLCT as
the square error if elements of XY are positive.
Proposition 5.2. Let the probability density functions of W ∈ M(M, N, K)
be q(W ) and p(W |X, Y ), which represent a true distribution and a learning
machine respectively defined by
q(W ) ∝

Exp(W |AB),

p(W |X, Y ) ∝

Exp(W |XY ),

where Exp(W |C) is a probability density function of the exponential distribution
with average C. Also let ϕ(X, Y ) be a probability density function such that it
is positive on a compact subset of M(M, H, K0 ) × M(H, N, K0 ). Then, the
Kullback-Leibler divergence has same RLCT as the square error.
Proof. Let x > 0, a > 0, b > 0. We put
p(x|a)

:=

K(a, b) :=

e−x/a
,
Z a

p(x|a)log

17

p(x|a)
dx.
p(x|b)

Using
log

p(x|a)
p(x|b)

and

Z
K(a, b) =
=

b
= log e−x/a+x/b
a
x x
= log b − log a − + ,
a
b
Z
p(x|a)dx = 1

xp(x|a)dx = E[x] = a,


x x
dx
p(x|a) log b − log a − +
a
b
a
log b − log a − 1 + .
b

Z

(9)

We immediately get
∂a K(a, b) =
∂b K(a, b) =
=

1 1
− ,
b a
a
1
−
b b2
b−a
b2

and increase or decrease of K. Thus this proposition can be proved in the same
way as Proposition 5.1 as follows.
Remark 5.1. The right side of equality(7) and (9) are respectively equal to ”Idivergence”[16] and ”Itakura-Saito-divergence”[17] which are used as criterion
of difference between observed matrix and reproduced matrix in NMF[4, 18, 19].
As a result, we can apply the Main Theorem to Bayesian learning if we use a
Kullback-Leibler divergence or square error as a criterion of difference between
XY and AB in cases where elements of matrices are generated by normal,
Poisson or exponential distributions. Thus, the upper bound of the average
of generalization errors can be clarified if the size and inner dimension of the
observed matrix and reproduced matrix.

5.4

Future Study

Lastly, we discuss future study in this subsection. There are two works.
The first is about tightness of our upper bound. The upper bound derived in
this paper is the first theoretical result about Bayesian NMF. Hence there is not
yet clarified its upper bound which is lower than our upper bound. However,
in this paper, we also showed that there are cases when the exact values are
18

obtained as Lemma 3.1 and Lemma 3.2. Based on the exact cases, it is expected
that the upper bound is tight, however, its mathematical proof is the future
study.
The second is about numerical experiments. In NMF, even in non-Bayesian
cases[1], the accurate simulation method is not yet established. The numerical
calculation of them strongly depends on the initial values. It seldom goes to the
minimum point because it has many local minima. Also in Bayesian cases[2], the
effective Markov Chain Monte Carlo method ( MCMC ) is not yet devised for
NMF, where MCMC is a method to construct posterior in numerical Bayesian
estimation. Thus we have no method by which can be used for verifying the
upper bound. Our theorem is the first theoretical base for construction of numerical calculation. It is the future study to make the numerical method for
Bayesian NMF based on our theoretical result.

6

Conclusion

The upper bound of a real log canonical threshold of the non-negative matrix
factorization is derived and an application to Bayesian learning is introduced.
Future work is to clarify the tightness of the upper bound mathematically and
make numerical method to consider numerical verification in Bayesian nonnegative matrix factorization.

Acknowledgments
This research was partially supported by the Ministry of Education, Science,
Sports and Culture in Japan, Grant-in-Aid for Scientific Research 15K00331.

Table
In this section, we show the table mentioned in Discussion 5.2. ”-” means that
there is no matrix satisfying the condition written in the table. ”( value )” is
denoted that this value is not lower bound of the RLCT of NMF ( i.e. our result
) since there is no non-negative matrix satisfying conditions such as H0 = 4 and
M = N = 3.

19

Table 1: Comparison the RLCT of NMF with that of reduced rank regression
model size and model / matrix size
M=N=2 M=N=3 M=N=4 M=N=5
H = M ,H0 = 0
NMF ( exact value )
2
9/2
8
25/2
r=0
reduced rank regression
3/2
7/2
6
19/2
H = H0 = 1
NMF ( exact value )
3/2
5/2
7/2
9/2
r=1
reduced rank regression
3/2
5/2
7/2
9/2
H = H0 = 2
NMF ( bound )
3
5
7
9
r=2
reduced rank regression
2
4
6
8
H = H0 = 3
NMF ( bound )
15/2
21/2
27/2
r=3
reduced rank regression
9/2
15/2
21/2
H = H0 = 4
NMF ( bound )
14
18
r=3
reduced rank regression
(9/2)
8
23/2
r=4
8
12
H = H0 = 5
NMF ( bound )
45/2
r=3
(9/2)
(8)
12
r=4
reduced rank regression
(8)
25/2
r=5
25/2

References
[1] P. Paatero, U. Tapper, Positive matrix factorization: A non-negative factor
model with optimal utilization of error estimates of data values, Environmetrics 5 (2) (1994) 111–126, doi:10.1002/env.3170050203.
[2] A. T. Cemgil, Bayesian inference in non-negative matrix factorisation models, Computational Intelligence and Neuroscience 2009 (4) (2009) 17, article
ID 785152.
[3] W. Xu, X. Liu, Y. Gong, Document clustering based on non-negative matrix factorization, in: Proceedings of the 26th annual international ACM
SIGIR conference on Research and development in information retrieval.
Association for Computing Machinery, 2003, pp. 267–273.
[4] D. D. Lee, H. S. Seung, Learning the parts of objects with nonnegative
matrix factorization, Nature 401 (1999) 788–791.
[5] Y. Xu, S. Deng, X. Li, Y. He, A sparse unmixing model based on NMF
and its application in raman image, Neurocomputing 207 (2016) 120 – 130.
[6] X. Jia, F. Sun, H. Li, Y. Cao, X. Zhang, Image multi-label annotation
based on supervised nonnegative matrix factorization with new matching
measurement, Neurocomputing 219 (2017) 518 – 525.
[7] H.Kim, H. Park, Sparse non-negative matrix factorizations via
alternating non-negativity-constrained least squares for microarray data analysis, Bioinformatics 23 (12) (2007) 1495–1502,
doi:10.1093/bioinformatics/btm134. PMID 17483501.
20

[8] M. Kohjima, T. Matsubayashi, H. Sawada, Probabilistic non-negative
inconsistent-resolution matrices factorization, in: Proceeding of CIKM ’15
Proceedings of the 24th ACM International on Conference on Information
and Knowledge Management, Vol. 1, 2015, pp. 1855–1858.
[9] S. Watanabe, Algebraic analysis for non-regular learning machines, Advances in Neural Information Processing Systems (12) (2000) 356–362, denver, USA.
[10] S. Watanabe, Algebraic geometrical methods for hierarchical learning machines, Neural Networks 13 (4) (2001) 1049–1060.
[11] M. Aoyagi, S. Watanabe, Stochastic complexities of reduced rank regression
in bayesian learning, Neural Networks (18) (2005) 924–933.
[12] K. Yamazaki, S. Watanabe, Singularities in mixture models and upper
bounds of stochastic complexity, International Journal of Neural Networks
16 (7) (2003) 1029–1038.
[13] K. Yamazaki, S. Watanabe, Algebraic geometry and stochastic complexity
of hidden markov models, Neurocomputing 69 (2005) 62–84, issue 1-3.
[14] M. Drton, M. Plummer, A bayesian information criterion for singular models, Journal of the Royal Statistical Society Series B 79 (2017) 323–380, with
discussion.
[15] J. E. Cohen, U. G. Rothblum, Nonnegative ranks, decompositions, and
factorizations of nonnegative matrices, Linear Algebra and Its Applications
190 (1993) 149–168.
[16] G. Tusnady, I. Csiszar, G. Tusnády, Information geometry and alternating
minimization procedures, Statistics and Decisions: Supplement Issues 1
(1984) 205–237.
[17] F. Itakura, S. Saito, Analysis synthesis telephony based on the maximum
likelihood method, in: In Proc. 6th of the International Congress on Acoustics, 1968.
[18] C. Févotte, N. Bertin, J.-L. Durrieu, Nonnegative matrix factorization with
the itakura-saito divergence: With application to music analysis, Neural
computation 21 (3) (2009) 793–830.
[19] L. Finesso, P. Spreij, Nonnegative matrix factorization and i-divergence
alternating minimization, Linear Algebra and its Applications 416 (2-3)
(2006) 270–287.

21

