arXiv:1711.07925v1 [eess.SP] 21 Nov 2017

Kullback-Leibler Principal Component for
Tensors is not NP-hard
Kejun Huang*

Nicholas D. Sidiropoulos†

Abstract
We study the problem of nonnegative rank-one approximation of a nonnegative
tensor, and show that the globally optimal solution that minimizes the generalized
Kullback-Leibler divergence can be efficiently obtained, i.e., it is not NP-hard.
This result works for arbitrary nonnegative tensors with an arbitrary number of
modes (including two, i.e., matrices). We derive a closed-form expression for the
KL principal component, which is easy to compute and has an intuitive probabilistic interpretation. For generalized KL approximation with higher ranks, the
problem is for the first time shown to be equivalent to multinomial latent variable
modeling, and an iterative algorithm is derived that resembles the expectationmaximization algorithm. On the Iris dataset, we showcase how the derived results
help us learn the model in an unsupervised manner, and obtain strikingly close
performance to that from supervised methods.

1 Introduction
Tensors are powerful tools for big data analytics [1], mainly thanks to the ability to
uniquely identify latent factors under mild conditions [2, 3]. On the other hand, most
detection and estimation problems related to tensors are NP-hard [4]. A similar situation is encountered in nonnegative matrix factorization, which is essentially unique
under certain conditions [5, 6] and computationally NP-hard [7]. In a lot of applications, nonnegativity constraints are natural for tensor latent factors as well.
In practice, the latent factors of tensors and matrices are usually obtained by minimizing the mismatch between the data and the factorization model according to certain
loss measures. The most popular loss measures the sum of the element-wise squared
errors, which is conceptually appealing and conducive for algorithm design, thanks to
the success of least-squares-based methods. For example, an effective algorithm for
minimizing the least-squares loss is AO-ADMM [8], and we refer the readers to the
references therein for other least-squares-based methods. From an estimation theoretic
point of view, the least-squares loss admits a maximum-likelihood interpretation under
i.i.d. Gaussian noise. In a lot of applications, however, it remains questionable whether
Gaussian noise is a suitable model for nonnegative data.
* University
† University

of Minnesota, Minneapolis, MN 55455. Email: huang663@umn.edu
of Virginia, Charlottesville, VA 22904. Email: nikos@virginia.edu

1

We study the problem of fitting a nonnegative data matrix/tensor Y with low rank
factors, using the generalized Kullback-Leibler (KL) divergence as the fitting criterion.
Mathematically, given a N -way tensor data Y ∈ RJ1 ×...×JN and a target rank K, we
try to find factor matrices constituting a canonical polyadic decomposition (CPD) that
best approximates the data tensor Y in terms of generalized KL divergence:
minimize
λ,{P (n)}

X

−Yj1 ...jN log

subject to λ ≥ 0, P

(n)

⊤

N
Y

λk

(n)

Pjn k

k=1

n=1

K
X

N
Y

j1 ,...,jN

+

K
X

λk

(n)
Pjn k

n=1

k=1
(n)

!

(1)

= 1⊤, ∀ n ∈ [N ].

≥ 0, 1 P

The conditions imposed in Problem (1) besides nonnegativity are intended for resolving
the trivial scaling ambiguity inherent in matrix factorization and tensor CPD models,
and thus are without loss of generality: the columns of all the factor matrices are normalized to sum up to one, and the scalings are absorbed into the corresponding values
in the vector λ. We adopt the common notation [[λ; P (1), ..., P (N)]] to denote the tensor
synthesized from the CPD model using these factors.
An instant advantage of this formulation is that the loss function of (1) can be
greatly simplified. Thanks to the constraints 1⊤P (n) = 1⊤, it is easy to see that
K
X X

j1 ,...,jN k=1

λk

N
Y

(n)
Pjn k

=

n=1

K
X

λk .

k=1

Therefore, Problem (1) is mathematically equivalent to
minimize −
λ,{P (n)}

X

Yj1 ...jN log

j1 ,...,jN

K
X

k=1

λk

N
Y

(n)

Pjn k + 1⊤λ

n=1

(2)

subject to λ ≥ 0, P (n) ≥ 0, 1⊤P (n) = 1⊤, ∀ n ∈ [N ].

2 Motivation: Probabilistic Latent Variable Modeling
Most of the existing literature motivates the use of generalized KL-divergence by modeling the nonnegative integer data as generated from Poisson distributions [9]. Specifically, the model states that each entry of the tensor Yj1 ...jN is generated from a Poisson
distribution with parameter Θj1 ...jN , and the underlying tensor Θ admits an exact CPD
model Θ = [[λ; P (1), ..., P (N)]]. While it is a simple and reasonable model, the physical
meaning behind the CPD model for the underlying Poisson parameters is not entirely
clear.
In this paper we give the choice of generalized KL-divergence as the loss function
a more compelling motivation. Consider N categorical random variables X1 , ..., XN ,
each taking J1 , ..., JN possible outcomes, respectively. Suppose these random variables are jointly drawn for a number of times, each time independently, and the outcome counts are recorded in a N -dimensional count tensor Y , which is the data we are
2

given. Denote the joint probability that X1 = j1 , ..., XN = jN as Πj1 ...jN , i.e.
P [X1 = j1 , ..., XN = jN ] = Πj1 ...jN ,
where we have

X

Πj1 ...jN = 1.

j1 ,...,jN

Then the overall probability that, out of M independent draws, the event of X1 =
j1 , ..., XN = jN occurs Yj1 ...jN times is
P [Y ] = Q

where

Y
M!
Y
(Πj1 ...jN ) j1 ...jN ,
(Yj1 ...jN !) j ,...,j
1

X

N

Yj1 ...jN = M.

(3)

j1 ,...,jN

The maximum likelihood estimates of the parameters are simply
Πj1 ...jN =

Yj1 ...jN
.
M

However, this simple estimate of Π may not be of much practical use. First of all, the
number of parameters we are trying to estimate is the same as the number of possible
outcomes, which is not a parsimonious model; in other words, we are not exploiting any possible structure between the random variables X1 , ..., XN , other than the
fact that exists a joint distribution between them. Furthermore, as a result of overparameterization, we will need the number of independent draws M to be very large
before we can have an accurate estimate of Π, which is often not the case in practice.
A simple and widely used assumption we can impose onto the set of variables is
the naive Bayes model: Suppose there is a hidden random variable Ξ, which is also
categorical and can take K possible outcomes, such that X1 , ..., XN are conditionally
independent given Ξ. The corresponding graphical model is given in Fig. 1. Mathematically, this means
P [X1 = j1 , ..., XN = jN |Ξ = k] =

N
Y

P [Xn = jn |Ξ = k] .

n=1

Then we have
P [X1 = j1 , ..., XN = jN ]
=

K
X

P [X1 = j1 , ..., XN = jN |Ξ = k] P[Ξ = k]

k=1

=

K
X

k=1

P[Ξ = k]

N
Y

P [Xn = jn |Ξ = k] .

n=1

3

Ξ

X1

...

X2

XN

Figure 1: A graphical model depicting the probabilistic dependencies between the hidden variable Ξ and the observed variables X1 , ..., XN , in a naive Bayes model.
Denote

(n)

λk = P[Ξ = k] and Pjn k = P [Xn = jn |Ξ = k] ,
then it is easy to see that
Π = [[λ; P (1), ..., P (N)]],

(4)

which means Π admits an exact CPD model. Even though naive Bayes seems to be a
very specific model, it has recently been shown that it is far more general than meets the
eye—no matter how dependent X1 , ..., XN are, there always exists a hidden variable
Ξ such that the depicted naive Bayes model holds, thanks to a link between tensors and
probability established in [10].
Using the CPD parameterization of the multinomial parameter Π, we formulate the
maximum likelihood estimation of P[Ξ = k] and P[Xn = jn |Ξ = k] as the following
optimization problem:
minimize −
λ,{P (n)}

X

Yj1 ...jN log

j1 ,...,jN

K
X

k=1
⊤

subject to λ ≥ 0, 1 λ = 1,

λk

N
Y

(n)

Pjn k

n=1

(5)

P (n) ≥ 0, 1⊤P (n) = 1⊤, ∀ n ∈ [N ].
Problem (5) is different from Problem (2), but the difference is small—in (5), λ is
constrained to sum up to one, whereas in (2), the sum of the elements of λ is penalized
in the loss function. In fact, the two problems are exactly equivalent, despite their
apparent differences.
(n)

(n)

Theorem 1. Let (λ⋆ , {P⋆ }) be an optimal solution for Problem (2), then (λ⋆ , {P⋆ })
is an optimal solution for Problem (5), where
λ⋆ =

1
λ⋆ .
1 λ⋆
⊤

Before proving Theorem 1, we first show the following lemma, which is interesting
in its own right.
Lemma 1. If λ⋆ is optimal for (2), then
1⊤λ =

X

j1 ,...,jN

4

Yj1 ...jN .

Proof. We show this by checking the optimality condition of (2) with respect to λ.
Without loss of generality, we can assume that λ⋆ > 0 strictly, because otherwise
the rank can be reduced. Since the inequality constraints with respect to λ are not
attained as equalities, their corresponding dual variables are equal to zero, according to
complementary slackness. The KKT condition for (2) with respect to λ then reduces
to the gradient of the loss function of (2) with respect to λ at λ⋆ being equal to zero.
Specifically, setting the derivative with respect to λk equal to zero yields
X

j1 ,...,jN

Therefore
K
X

λ⋆k

k=1

P

N
Y
Yj1 ...jN
(n)
Pjn k = 1.
Q (ν)
P
λ
ν jν k n=1
κ ⋆κ


N
Y
Y
j
...j
(n)
1
N
λ⋆k
Pjn k 
=
Q (ν)
P
P
λ
⋆κ
j1 ,...,jN
k=1
ν jν k n=1
κ
!
K
N
X
X
Y
Yj1 ...jN
(n)
=
λ⋆k
Pjn k
Q (ν)
P
n=1
j1 ,...,jN
k=1
ν Pjν k
κ λ⋆κ
X
Yj1 ...jN .
=
K
X



X

j1 ,...,jN

We now prove Theorem 1 with the help of Lemma 1.
(n)

Proof of Theorem 1. We show that (λ⋆ , {P⋆ }) is optimal for (5) via contradiction.
(n)
(n)
Suppose (λ⋆ , {P⋆ }) is not optimal for (5), then there exists a feasible point (λ♭ , {P♭ })
such that
−

X

Yj1 ...jN log

j1 ,...,jN

K
X

λ♭k

X

(n)

P♭jn k

n=1

k=1

<−

N
Y

Yj1 ...jN log

j1 ,...,jN

K
X

k=1

λ⋆k

N
Y

(n)

P⋆jn k

n=1

(n)

Define λ♭ = M λ♭ , then (λ♭ , {P♭ }) is clearly feasible for (2). Furthermore, we have
−

X

Yj1 ...jN log

j1 ,...,jN

=−

X

K
X

k=1

Yj1 ...jN log

j1 ,...,jN

<−

X

X

j1 ,...,jN

K
X

N
Y

Yj1 ...jN log

K
X

λ♭k

Yj1 ...jN log

N
Y

(n)

P♭jn k − M log M + M

n=1

λ⋆k

N
Y

(n)

P⋆jn k − M log M + M

n=1

k=1
K
X

(n)

P♭jn k + 1⊤λ♭

n=1

k=1

j1 ,...,jN

=−

λ♭k

λ♭k

N
Y

n=1

k=1

5

(n)

P⋆jn k + 1⊤λ⋆ ,

(n)

where the equalities stems from Lemma 1. This means (λ♭ , {P♭ }) gives a smaller loss
(n)
(n)
value for (2) than that of (λ⋆ , {P⋆ }), and contradicts our assumption that (λ⋆ , {P⋆ })
is optimal for (2).
A similar but less general result for the case when N = 2 is given in [11], in the
context of nonnegative matrix factorization using generalized KL-divergence loss.
The take home point of this section is that we can find the maximum likelihood estimate of the hidden variable in the naive Bayes model by taking the nonnegative CPD
of the data tensor. In a way, our analysis suggests that the generalized KL-divergence
is a more suitable loss function to fit nonnegative data than, for example, the L2 loss.

3 KL Principal Component
Now that we have established how important Problem (2) is in probabilistic latent variable modeling, we focus on a specific case of (2) when K = 1. This case corresponds
to extracting the “principal component” of a nonnegative tensor under the generalized
KL-divergence loss. We will show, in this section, that this specific problem is not only
computationally tractable, but admits an extremely simple closed form solution.
Let us first rewrite Problem (2) with K = 1. In this case, the diagonal loadings λ
become a scalar λ, and the individual factor matrices P (n) become vectors p(n):
minimize −
λ,{P (n)}

X

Yj1 ...jN log λ

(n)

pjn + λ

n=1

j1 ,...,jN

subject to λ ≥ 0, p

N
Y

(n)

⊤ (n)

≥ 0, 1 p

(6)

= 1, ∀ n ∈ [N ].

A salient feature of this case when K = 1 is that there is no summation in the log, just
a product. Therefore, we can equivalently write Problem (6) as
!
N
X
X
(n)
log pjn + λ
Yj1 ...jN log λ +
minimize −
λ,{P (n)}
(7)
n=1
j1 ,...,jN
subject to λ ≥ 0, p(n) ≥ 0, 1⊤p(n) = 1, ∀ n ∈ [N ].

Noticing that − log is a convex function, the exciting observation we see in Problem (7)
is that it is in fact convex! This already means that it can be solved optimally and
efficiently [12], but we will in fact show a lot more: that it admits very simple and
intuitive closed-form solution.
Remark. It may seem obvious from our derivation that once we try rewriting Problem (2) with K = 1, one will immediately see that this problem has hidden convexity
in it. However, to the best of our knowledge, we are the first to point out this fact. In
hindsight, there is a subtle caveat that plays a key role in spotting this hidden convexity:
we fixed the inherent scaling ambiguities by constraining all the columns of the factor
matrices to sum up to one; as a result, the sum of all the entries of the reconstructed
tensor, which appears in the generalized KL-divergence loss, boils down to simply the
sum of the diagonal loadings. If this trick is not applied to the problem formulation,
6

there is still a multi-linear term in the loss of (7), and the hidden convexity is not at all
obvious.
We now derive an optimal solution for Problem (7) by checking the KKT conditions. For convex problems such as (7), KKT condition is not only necessary, but also
sufficient for optimality. Using Lemma 1, we immediately have that
X
Yj1 ...jN = M.
λ=
j1 ,...,jN

For the variable p(n), let ξ (n) denote the dual variable corresponding to the equality
constraint 1⊤p(n) = 1, and q (n) denote the nonnagative dual variable corresponding
to the inequality constraint p(n) ≥ 0, for all n = 1, ..., N . The loss function in (7)
(n)
separates down to the components, and the term corresponding to pjn is
−

X

(n)

(n)

(n)

Yj1 ...jN log pjn = −yjn log pjn ,

j1 ,...,jn−1 ,
jn+1 ,...,jN

where we denote

X

(n)

yjn =

Yj1 ...jN .

(8)

j1 ,...,jn−1 ,
jn+1 ,...,jN
(n)

(n)

We see that pjn cannot be equal to zero if yjn 6= 0, because otherwise it will drive
the loss value up to +∞; according to complementary slackness, this means the corre(n)
(n)
(n)
sponding qjn = 0. If yjn = 0, then pjn does not directly affect the loss value of (7),
even when it equals to zero, using the convention that 0 log 0 = 0. Since we have the
(n)
constraint 1⊤p(n) = 1, such a pjn should be equal to zero at optimality, otherwise the
(n)
other entries in p will be smaller, leading to a larger loss value in (7).
(n)
(n)
Suppose yjn 6= 0. Setting the derivative of the Lagrangian with respect to pjn equal
to zero, we have
(n)
yjn
(n)
− qjn + ξ (n) = 0.
− (n)
pjn
(n)

As we explained, qjn = 0 according to complementary slackness. Therefore,
(n)

(n)

pjn = yjn /ξ (n).
The dual variable ξ (n) should be chosen so that the equality constraint 1⊤p(n) = 1
(n)
(n)
is satisfied. Together with our argument that pjn = 0 if yjn = 0, we come to the
conclusion that
(n)
(n)
yj
yj
(n)
pjn = ⊤ n(n) = n .
M
1 y
The result we derived in this section is summarized in the following theorem.

7

Theorem 2. The KL principal component of a nonnegative tensor Y , i.e., the solution
to Problem (6), is
(n)
yj
(n)
λ = M and pjn = n ,
(9)
M
(n)

where M and yjn are defined in (3) and (8), respectively.
Now let us take a deeper look at the solution we derived for Problem (6). Suppose
the data tensor Y is generated by drawing from the joint distribution P[X1 , ..., XN ]
M times. Our definition of y (n) in (8) essentially summarizes the number of times
each possible outcomes of Xn occurs, regardless of the outcomes of the other random
variables. As a result, the optimal KL principal component factor p(n) is in fact the
maximum likelihood estimate of the marginal distribution P[Xn ]. On hindsight, the
simple solution we provide for KL principal component in Theorem 2 becomes very
natural. The case when K = 1 means the latent variable Ξ can only take one possible
outcome with probability one, which means Ξ is not random. In other words, we are
basically assuming that X1 , ..., XN are independent from each other. As a result, the
joint distribution factors into the product of the marginal distributions
P [X1 , ..., XN ] =

N
Y

P[Xn ],

n=1

and the marginal distributions can be simply estimated by “marginalizing” the observations, collected in Y , and then normalizing to sum up to one. This is elementary in
probability. However, in the context of finding the principal component of a nonnegative tensor using generalized KL-divergence, it is not at all obvious. Furthermore, the
argument based on categorical random variables only applies to nonnegative integer
data, whereas our derivation of the KL principal component of a tensor is not restricted
to integers or rational numbers, but works for general real nonnegative numbers as well.

4 KL Approximation with Higher Ranks
If K > 1 in Problem (2), there is more than one term in the logarithm; as a result, the
nice transformation from (6) to (7) cannot be directly applied. There is a way to apply
something similar, and that is through the use of Jensen’s inequality [13] (applied to
the − log function)
− log Ez ≤ −E log z.
We use this inequality to define majorization functions for the design of an iterative
upperbound minimization algorithm [14].
Suppose at the end of iteration t, the obtained updates are λt and {P (n)t }. At the
next iteration, we define
Ψjt1 ...jN k

=

λtk

N
Y

(n)t
Pjn k

n=1

X
K

κ=1

8

λtκ

N
Y

n=1

(n)t

Pjn κ .

(10)

According to this definition, it is easy to see that Ψjt1 ...jN k ≥ 0 and
Assuming Ψjt1 ...jN k > 0, we have
− Yj1 ...jN log

K
X

λk

≤ −Yj1 ...jN

Ψjt1 ...jN k = 1.

(n)

Pjn k

K
X
Ψjt

1 ...jN k

k=1
K
X

k=1

n=1

k=1

= −Yj1 ...jN log

N
Y

PK

Ψjt1 ...jN k

λk

Ψjt1 ...jN k log λk

N
Y

(n)

Pjn k

n=1
N
Y

(n)

Pjn k + const.

(11)

n=1

k=1

Furthermore, equality is attained if λ = λt and P (n) = P (n)t , ∀n ∈ [N ]. This defines a
majorization function for iteration t+1, and the minimizer of (11) is set to be the update
of this iteration. Since (11) and the loss function of (2) are both smooth, the convergence result of the successive upperbound minimization (SUM) algorithm [14] can be
applied to establish that this procedure converges to a stationary point of Problem (2).
We should stress again that this procedure is made easy only after the multi-linear
term in (1) is equivalently replaced by the sum of the diagonal loadings in (2) through
our careful problem formulation. Otherwise, the multi-linear term still remains, which
together with (11) does not end up being a simpler function to optimize.
The majorization function (11) is nice, not only because it is convex, but also since
it is reminiscent of the loss function (7) when K is equal to one—it decouples the variables down to the canonical components, i.e., λk and the k-th columns of P (1), ..., P (N);
each of the sub-problems takes the form of (7), replacing Yj1 ...jN with Yj1 ...jN Ψjt1 ...jN k .
As a result, the update for iteration t + 1 boils down to something similar to what we
have derived in the previous iteration. Specifically, define
X
Yj1 ...jN Ψjt1 ...jN k
Mkt =
j1 ,...,jN

and

(n)t

yjn k =

X

Yj1 ...jN Ψjt1 ...jN k ,

j1 ,...,jn−1 ,
jn+1 ,...,jN

then
λt+1
= Mkt
k

(n)t+1

and Pjn k

.
(n)t
= yjn k Mkt .

(12)

A nice probabilistic interpretation can be made to understand this algorithm. In the
context of probabilistic latent variable modeling, the conditional distributions P[Xn |Ξ]
can be easily estimated if the joint observation X1 = j1 , ..., XN = jN , and Ξ = k
is given, because then we can collect all the observations with Ξ = k and use the
techniques derived in the previous section. Now that we do not observe Ξ, what we can
do is to try to estimate P[Ξ|X1 , ..., XN ] instead. The Ψ defined in (10) does exactly
the job using the current estimate of λ and {P (n)}. Using this estimated posterior
distribution of Ξ, we have a guess of the portion of Yj1 ...jN that jointly occurs with Ξ =
9

k, and use that to obtain a new estimate of P[Ξ] and P[Xn |Ξ]. This is exactly the idea
behind the expectation-maximization (EM) algorithm [15]. Almost the same algorithm
has been derived by Shashanka et al. [16], and the special case when N = 2 is the EM
algorithm for probabilistic latent semantic analysis (pLSA) [17]. However, we should
mention that these algorithms were originally derived in the context of multinomial
latent variable modeling, and without the help of Theorem 1, it was not previously
known that they can be used for generalized KL-divergence fitting as well.
Computationally, although the definition of Ψ in (10) helps us obtain simple expressions (12) and intuitive interpretations, we do not necessarily need to explicitly
form them when implementing the method for memory/computation efficiency considerations. To calculate λt+1 and {P (n)t+1 } as in (12), we first define Ye t such that
Yejt1 ...jN = Yj1 ...jN

X
K

κ=1

λtκ

N
Y

(n)t

Pjn κ .

n=1

This operation requires passing through the data once, and if Y is sparse, Ye t has
exactly the same sparsity structure. Then,
(1)t

(N)t

λt+1
= λtk Ye t ×1 pk ...×N pk ,
k

(13)
(n)t

where ×n denotes the n-mode tensor-vector multiplication [18], and pk
k-th column of P (n)t . As for the factor matrices,


P (n)t+1 = P (n)t ∗ M TTKRP Ye t , {P (ν)t }, n

denotes the

(14)

followed by column normalization to satisfy the sum-to-one constraint, where ∗ denotes matrix Hadamard (element-wise) product, and M TTKRP stands for the n-mode
matricized tensor times Khatri-Rao product of all the factor matrices {P (ν)t } except
the n-th one.
It is interesting to notice that the update rules (13) and (14) are somewhat similar
to the widely used multiplicative-update (MU) rule for NMF [19] and NCP [9]. The
big difference lies in the fact that MU updates the factors alternatingly, whereas EM
updates the factors simultaneously. This makes the EM algorithm extremely easy to
parallelize—for N -way factorizations, we can simply take N cores, each taking care
of the computation for the n-th factor, and we can expect an almost ×N acceleration,
if the sizes of all the modes are similar.

5 Illustrative Example
We give an illustrative example using the Iris data set downloaded from the UCI Machine Learning repository. The data set contains 3 classes of 50 instances each, where
each class refers to a type of iris plant. For each sample, four features are collected:
sepal length/width and petal length/width in cm. The measurements are discretized into
0.1cm intervals, and the four features range between 4.3 ∼ 7.9, 2.0 ∼ 4.4, 1.0 ∼ 6.9,
and 0.1 ∼ 2.5cm, respectively. We therefore form a 4-way tensor with dimension
10

sepal length in cm

0.2

sepal length in cm

0.2

0.1

0.1

0

0
5

6

7

5

sepal width in cm

6

7

sepal width in cm

0.2

0.2

0.1

0.1

0

0
2

2.5

3

3.5

4

2

petal length in cm

2.5

3

3.5

4

petal length in cm

0.2

0.2

0.1

0.1

0

0
2

4

6

2

petal width in cm

4

6

petal width in cm

0.5

0.5

0

0
0.5

1

1.5

2

2.5

0.5

supervised

1

1.5

2

2.5

unsupervised

Figure 2: Learned P[feature|class] in three colors for the three iris classes. Left: class
label is given. Right: class label not given.
37 × 25 × 60 × 25; for a new data sample, the corresponding entry in that tensor is
added with one.
Suppose the four features are conditionally independent given the class label. We
can then collect all the samples from the same class into one tensor, and invoke Theorem 2 to estimate their individual conditional distributions P[feature|class], which are
shown on the left panel of Fig. 2. This learned conditional distribution can then be used
to classify new samples, which is the idea behind the naive Bayes classifier [20].
Now suppose the class labels are not given to us. We then collect all the data samples into a single tensor. We still assume that the features are conditionally independent
given the class label, even though the class label is now latent (unobserved). As per our
discussion in Section 4, we can still try to estimate the conditional distribution using
the EM algorithm (12). We run the EM algorithm from multiple random initializations,

11

and the result that gives the smallest loss is shown on the right of Fig. 2.
The astonishing observation is that the learned conditional distribution is almost
identical to the one learned when the class label is given to us. This suggests that the
conditional independence between features given class labels is actually a reasonable
assumption in this case, contrary to the common belief that naive Bayes is an “oversimplified” model. It is also interesting to notice that this nice result is obtained using
only 150 data samples, which is extremely small considering the size of the tensor.

6 Conclusion
We studied the nonnegative CPD problem with generalized KL-divergence as the loss
function. The most important result of this paper is the discovery that finding the KL
principal component of nonnegative tensor is not NP-hard. To make matters nicer, we
derived a very simple closed-form solution for finding the KL principal component.
This is a surprisingly pleasing result, considering that in the field of tensors “most
problems are NP-hard” [4]. Borrowing the idea for finding the KL principal component, an iterative algorithm for higher rank KL approximation was also derived, which
is guaranteed to converge to a stationary point and is easily and naturally parallelizable.

References
[1] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and
C. Faloutsos, “Tensor decomposition for signal processing and machine learning,”
IEEE Transactions on Signal Processing, vol. 65, no. 13, pp. 3551–3582.
[2] J. B. Kruskal, “Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics,” Linear algebra
and its applications, vol. 18, no. 2, pp. 95–138, 1977.
[3] N. D. Sidiropoulos and R. Bro, “On the uniqueness of multilinear decomposition
of n-way arrays,” Journal of chemometrics, vol. 14, no. 3, pp. 229–239, 2000.
[4] C. J. Hillar and L.-H. Lim, “Most tensor problems are NP-hard,” Journal of the
ACM, vol. 60, no. 6, p. 45, 2013.
[5] K. Huang, N. D. Sidiropoulos, and A. Swami, “Non-negative matrix factorization revisited: Uniqueness and algorithm for symmetric decomposition,” IEEE
Transactions on Signal Processing, vol. 62, no. 1, pp. 211–224, 2014.
[6] X. Fu, K. Huang, and N. D. Sidiropoulos, “On identifiability of nonnegative matrix factorization,” IEEE Signal Processing Letters, 2017 (submitted).
[7] S. A. Vavasis, “On the complexity of nonnegative matrix factorization,” SIAM
Journal on Optimization, vol. 20, no. 3, pp. 1364–1377, 2009.
[8] K. Huang, N. D. Sidiropoulos, and A. P. Liavas, “A flexible and efficient algorithmic framework for constrained matrix and tensor factorization,” IEEE Transactions on Signal Processing, vol. 64, no. 19, pp. 5052–5065.
12

[9] E. C. Chi and T. G. Kolda, “On tensors, sparsity, and nonnegative factorizations,”
SIAM Journal on Matrix Analysis and Applications, vol. 33, no. 4, pp. 1272–
1299, 2012.
[10] N. Kargas and N. D. Sidiropoulos, “Completing a joint PMF from projections:
A low-rank coupled tensor factorization approach,” in Information Theory and
Applications Workshop (ITA), Feb 2017, pp. 1–6.
[11] N.-D. Ho and P. Van Dooren, “Non-negative matrix factorization with fixed row
and column sums,” Linear Algebra and Its Applications, vol. 429, no. 5-6, pp.
1020–1025, 2008.
[12] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge university press,
2004.
[13] J. Jensen, “Sur les fonctions convexes et les inégalités entre les valeurs
moyennes,” Acta mathematica, vol. 30, no. 1, pp. 175–193, 1906.
[14] M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence analysis
of block successive minimization methods for nonsmooth optimization,” SIAM
Journal on Optimization, vol. 23, no. 2, pp. 1126–1153, 2013.
[15] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” Journal of the royal statistical society. Series B
(methodological), pp. 1–38, 1977.
[16] M. Shashanka, B. Raj, and P. Smaragdis, “Probabilistic latent variable models as
nonnegative factorizations,” Computational Intelligence and Neuroscience, 2008.
[17] T. Hofmann, “Probabilistic latent semantic analysis,” in Uncertainty in Artificial
Intelligence (UAI), 1999, pp. 289–296.
[18] B. W. Bader and T. G. Kolda, “Efficient matlab computations with sparse and
factored tensors,” SIAM Journal on Scientific Computing, vol. 30, no. 1, pp. 205–
231, 2007.
[19] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in
Advances in Neural Information Processing Systems (NIPS), 2001, pp. 556–562.
[20] A. Y. Ng and M. I. Jordan, “On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes,” in Advances in Neural Information
Processing Systems (NIPS), 2002, pp. 841–848.

13

