Infinite Positive Semidefinite Tensor Factorizationfor Source Separation of Mixture Signals

Kazuyoshi Yoshii
k.yoshii@aist.go.jp
National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki 305-8568, Japan
Ryota Tomioka
tomioka@mist.i.u-tokyo.ac.jp
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Daichi Mochihashi
daichi@ism.ac.jp
The Institute of Statistical Mathematics (ISM), 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan
Masataka Goto
m.goto@aist.go.jp
National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki 305-8568, Japan

Abstract
This paper presents a new class of tensor factorization called positive semideﬁnite tensor
factorization (PSDTF) that decomposes a setof positive semideﬁnite (PSD) matrices into
the convex combinations of fewer PSD basis
matrices. PSDTF can be viewed as a natural extension of nonnegative matrix factorization. One of the main problems of PSDTF isthat an appropriate number of bases should
be given in advance. To solve this problem,
we propose a nonparametric Bayesian model
based on a gamma process that can instantiate only a limited number of necessary bases
from the inﬁnitely many bases assumed toexist. We derive a variational Bayesian algorithm for closed-form posterior inference anda multiplicative update rule for maximumlikelihood estimation. We evaluated PSDTF
on both synthetic data and real music recordings to show its superiority.

1. Introduction
Matrix factorization (MF) has recently been an active
research topic in the ﬁeld of machine learning. Givena matrix X ∈ RM ×N as observed data, the objective
is to ﬁnd a low-rank approximation X ≈ AB whereA ∈ RM ×K , B ∈ RK×N , and K  min{M, N }. This
problem often arises in many application ﬁelds. CanProceedings of the 30 th International Conference on Machine Learning, Atlanta, Georgia, USA, 2013. JMLR:
W&CP volume 28. Copyright 2013 by the author(s).

didates of X include a user-item rating matrix in collaborative ﬁltering (Salakhutdinov & Mnih, 2008), a
set of face images in image processing (Lee & Seung,2000), and a time-frequency spectrogram in audio processing (Smaragdis & Brown, 2003). Many variants ofMF have been proposed by using various measures on
the reconstruction error D(X|AB) and imposing constraints on A and B. In terms of probabilistic modeling, a speciﬁc model is deﬁned by a likelihood functionp(X|A, B) and prior distributions p(A) and p(B).
One of the popular classes of MF is nonnegative matrix
factorization (NMF), in which all elements of A and B
must be no less than zero. This constraint reﬂects thefact that some physical quantities, e.g., pixel brightness and signal energy, cannot be negative. A typical
way to impose this constraint is to place element-wise
gamma priors on A and B (Cemgil, 2009). D(X|AB)
has often been deﬁned in an element-wise manner byusing the Bregman divergence (Bregman, 1967), which
includes as special cases the Kullback-Leibler (KL) divergence (Kullback & Leibler, 1951) and the ItakuraSaito (IS) divergence (Itakura & Saito, 1968). The assumption underlying the likelihood p(X|A, B) is thateach element of X is independently Poisson distributed
in KL-NMF or independently exponentially distributedin IS-NMF.
In audio analysis, IS-NMF is more suitable for decomposing a power spectrogram X over M frequency bins
and N frames as the product of sound-source power
spectra (K columns of A) and the corresponding temporal activations (K rows of B) (Févotte et al., 2009).
IS-NMF is theoretically justiﬁed if the frequency binsof source spectra are independent. Note that the au-

Infinite Positive Semidefinite Tensor Factorization

2. Gamma Process Positive
Semidefinite Tensor Factorization
Second mode

All elements
are nonnegative…

First modeFirst mode

Second mode

Diagonal
elements arenonnegative…

: A set of nonnegative vectors:

A set of positive
semidefinite matrices

Figure 1. Comparison between IS-NMF and LD-PSDTF.

dio signals of pitched instruments have clear periodicities, i.e., those signals are highly autocorrelated at certain time lags. However, the short-term Fourier transform (STFT) is unable to perfectly decorrelate the frequency components forming harmonic structures. A
similar problem arises in electroencephalogram (EEG)
analysis (Lee et al., 2006), in which cross-correlations
between multichannel signals recorded at diﬀerent positions of the head are usually ignored. This indicatesthat it is not appropriate to place gamma priors on A
and deﬁne D(X|AB) in an element-wise manner.
To solve this problem, we propose a new class of tensor
factorization called positive semideﬁnite tensor factorization (PSDTF). As NMF decomposes N nonnegative
vectors (a matrix) as the conic sums of K nonnegative
vectors, PSDTF decomposes N PSD matrices (a tensor) as the conic sums of K PSD matrices. As shownin Figure 1, each nonnegative vector is embeded into a
PSD matrix that represents the covariance structure of
the multivariate elements. We thus place matrix-wise
Wishart priors on the basis matrices. In this paper thereconstruction error is deﬁned by using a kind of the
Bregman matrix divergence called the log-determinant(LD) divergence (Kulis et al., 2009), also in a matrixwise manner. This implies that each slice of the observed tensor is assumed to have a Wishart likelihood.
Since the resulting LD-PSDTF is a natural extensionof IS-NMF, an inherited problem is that the numberof bases, K, should be given in advance.
To estimate an appropriate number of basis matrices, we propose a nonparametric Bayesian model of
LD-PSDTF similar to one of IS-NMF (Hoﬀman et al.,2010). Although the Wishart prior-Wishart likelihood
hierarchy does not satisfy the conjugacy condition, we
can derive an elegant variational algorithm for closedform posterior inference. A multiplicative update rulecan also be used for maximum-likelihood estimation.
In addition, we reveal that IS-NMF in the frequency
domain approximates LD-PSDTF in the time domainthat can consider periodic covariance structures of theaudio signals. This explains why IS-NMF works wellfor music power-spectrogram decomposition.

We propose a probabilistic model of positive semideﬁnite tensor factorization (PSDTF) and derive its nonparametric Bayesian extension that in theory allows
observed data (e.g., a music signal) to contain an inﬁnite number of latent bases (e.g., sound sources) byusing the gamma process. An eﬀective number of bases
required for representing the observed data can be efﬁciently estimated in a data-driven manner. We thendiscuss how PSDTF is related to matrix factorizationand how it is applied to music signal analysis.
2.1. Problem Speciﬁcation
We will formalize the problem. Suppose we have as observed data a three-mode tensor X = [X 1 , · · · , X N ] ∈RM ×M ×N , where each slice X n ∈ RM ×M is a real
symmetric positive semideﬁnite (PSD) matrix. Although PSDTF can be deﬁned even if X n ∈ CM ×M
is a complex Hermitian PSD matrix, we focus on thecase of X n ∈ RM ×M for explanatory simplicity.
The goal of factorization is to approximate each PSD
matrix X n by a convex combination of PSD matrices{V k }K
k=1 (K bases) as follows:
Xn ≈

K


θk hkn V k ≡ Y n ,

(1)

k=1

where θk ≥ 0 is a global weight shared over all N slices
and hkn ≥ 0 is a local weight speciﬁc to the n-th slice.
Eq. (1) can also be represented as X ≈ K
k=1 θk hk ⊗V k ≡ Y , where ⊗ indicates the Kronecker product.
To evaluate the reconstruction error between PSD matrices X n and Y n , we propose to use a Bregman matrix divergence (Bregman, 1967) deﬁned as follows:
Dφ (X n |Y n )= φ(X n ) − φ(Y n ) − tr ∇φ(Y n )T (X n − Y n ) , (2)where φ is a strictly convex matrix function. In thispaper we focus on the log-determinant (LD) divergence(φ(Z) = − log |Z|) (Kulis et al., 2009) given by
DLD (X n |Y n ) + tr X n Y −1
− M.
= − log X n Y −1
n
n

(3)

This divergence is always nonnegative and is zero if
and only if X n = Y n holds. A well-known special casewhen M = 1 is the Itakura-Saito (IS) divergence overnonnegative numbers (Itakura & Saito, 1968) given by
DIS (x|y) = − log(x/y) + x/y − 1 and is often usedin signal processing. Sivalingam et al. (2010) formalized a similar tensor factorization problem that uses
DLD (Y n |X n ) as a cost function.

Infinite Positive Semidefinite Tensor Factorization

Our goal here is to estimate unknown variables θ =[θ1 , · · · , θK ]T ∈ RK , H = [h1 , · · · , hK ] ∈ RN ×K , andV = [V 1 , · · · , V K ] ∈RM ×M ×K such that the costfunction C(X|Y ) =n DLD (X n |Y n ) is minimized.
Note that our model imposes the nonnegativity constraint on θ and H and the positive semideﬁnitenessconstraint on V . We call this model LD-PSDTF.
2.2. Probabilistic Formulation
We explain Bayesian treatment of LD-PSDTF deﬁned
by Eq. (1) in terms of probabilistic modeling.
2.2.1. Formulating a Probabilistic ModelWe ﬁrst formulate a ﬁnite model based on a ﬁxed number of bases by specifying prior distributions on θ, H,
and V and a likelihood function of X. In this model
we assume θk = 1 because the eﬀect of θk can be compensated by adjusting the scale of hk .
Since hkn is nonnegative (hkn ≥ 0) and V k is PSD(V k ≥ 0), a natural choice is to place gamma andWishart priors on hkn and V k as follows:hkn ∼ G(a0 , b0 ),V k ∼ W(ν0 , V 0 ),

(4)
(5)

where a0 and b0 are the shape and rate parameters of
the gamma distribution and ν0 and V 0 are the degreeof freedom (DOF) and scale matrix of the Wishart
distribution.
We then assume PSD matrices {νX n }N
n=1 to be independently Wishart distributed as follows:

 K


θk hkn V k ,(6)
νX n θ, H, V ∼ W ν,
k=1

where ν is a DOF of the Wishart distribution. Note−1
Y n , where M
that E[X n ] = Y n and M[X n ] = ν−Mνmeans the mode. When ν
M , M[X n ] ≈ Y n holds.
When ν < M , X n is rank deﬁcient. If M = ν = 1,
the distribution reduces to an exponential distribution.
The log-likelihood of X n is given byν −M −1
log p(X n |Y n ) = C(ν) +log |X n |
2


ν
ν
, (7)
− log |Y n | − tr X n Y −1n22
where C(ν) is a constant term depending only on ν andthe second term can also be considered to be constantbecause X n is the observed data. Therefore, the maximization of the likelihood p(X|Y ) = n p(X n |Y n )with respect to Y is equivalentto the minimization of
the cost function C(X|Y ) = n DLD (X n |Y n ) (compare Eq. (7) with Eq. (3)).

Consequently, a Bayesian model of LD-PSDTF is deﬁned by Eqs. (4), (5), and (6). Given the data X, ourgoal is to calculate a posterior distribution p(H, V |X)
over unknown variables H and V .
2.2.2. Taking the Infinite Limit
To overcome the limitation that the number of bases
K should be speciﬁed in advance, we leverage Bayesiannonparametrics for taking the inﬁnite limit of Eq. (6)
as K diverges to inﬁnity. Given the data X, an eﬀective number of bases should be estimated in a datadriven manner. We thus aim to learn a sparse inﬁnitedimensional weight vector θ = [θ1 , · · · , θ∞ ]T as proposed by Hoﬀman et al. (2010).
We place a gamma process (GaP) prior on θ in a socalled weak-approximation manner as follows:θk ∼ G(αc/K, α),(8)where α andc are positive numbers, Eprior [θk ] = c/K,
and Eprior [ k θk ] = c. When the truncation level K
goes to inﬁnity, the vector θ approximates an inﬁnitedimensional discrete measure G that is stochasticallydrawn from the GaP over a space Θ as follows:G ∼ GaP(α, G0 ),(9)where α is called a concentration parameter and G0 abase measure. In our model we assumed that G0 is a
uniform measure such that G0 (Θ) = c. The eﬀective
number of elements, K + , such that θk >  for somenumber  > 0 is almost surely ﬁnite. If we set K to be
suﬃciently larger than α, only a few of the K elementsof θ will be substantially greater than zero.
A nonparametric Bayesian model of GaP-LD-PSDTF
is deﬁned by Eqs. (4), (5), (6), and (8) with a large
truncation level K. Given the data X, our goal is tocalculate a posterior distribution p(θ, H, V |X) andestimate the value of K + at the same time.
2.3. “Augmented” Matrix Factorization
We show here that LD-PSDTF naturally emerges fromthe standard problem of matrix factorization. Suppose
we have a set of N samples X̂ = [x̂1 , · · · , x̂N ] ∈ RM ×N
as observed data, where x̂n ∈ RM is a feature vector of
the n-th sample. Although the case of x̂n ∈ CM canbe dealt with, as in Section 2.1 we here discuss thecase of x̂n ∈ RM . In signal processing, for example, alocal signal sn ∈ RM in the n-th short window (calleda frame) is often regarded as x̂n (Figure 2). Alternatively, x̂n can be a complex spectrum cn = F sn ∈ CM
at the n-th frame, where F ∈ CM ×M is the unitary
discrete Fourier transform matrix.

Infinite Positive Semidefinite Tensor Factorization
Time-domain representation

Frequency-domain representation…
…
Local signal

Frequency

Framesin thex̂n |θ̂, Ŵ , Ĥ ∼ N-th frame

Frame

Power spectrum

Figure 2. Diﬀerent representations of a music audio signal.

The goal is to discover a limited number of bases having characteristic structures (e.g., instrument sounds)
from the observed data X̂ (e.g., music audio signal),i.e., to decompose each sample x̂n into a linear sum of
K variable bases {ŵ kn }K
k=1 as follows:K
k=1x̂nk =

K


0,

K



θ̂k2 ĥ2kn V̂ k

.

(14)

k=1

Fourier transformx̂n =likelihood of x̂n as follows:θ̂k ĥkn ŵ kn ,(10)

If we assume that θk = θ̂k2 ≥ 0, hkn = ĥ2kn ≥ 0, V k =V̂ k ≥ 0, and X n = x̂n x̂Tn ≥ 0, Eq. (14) recovers
Eq. (7) when ν = 1 except for constant terms. We canput the same priors as Eqs. (4), (5), and (8).
This is a special case of the general LD-PSDTF modelin which each X n is restricted to a rank-1 PSD matrix(X n = x̂n x̂Tn ). In general, the DOF of X can be larger
than M N because each X n is allowed to take any PSD
matrix. In this section, the DOF of X is M N because
X n is just an augmented representation of x̂n .
2.3.2. Estimating the Latent Components

k=1

where θ̂k is a global coeﬃcient of the k-th basis, ĥkn is alocal coeﬃcient of the k-th basis, and x̂nk = θ̂k ĥkn ŵkn
is the k-th component in x̂n . Those variables are allowed to take any real values. If we assume that ŵ kn isa ﬁxed basis such that ŵ kn is equal to ŵ k for any n anddeﬁne some symbols as θ̂ = [θ̂1 , · · · , θ̂K ]T ∈ RK , Ŵ =[ŵ 1 , · · · , ŵ k ] ∈ RM ×K , and Ĥ = [ĥ1 , · · · , ĥK ] ∈RN ×K , Eq. (10) can be simply written as follows:TX̂ = Ŵ diag(θ̂)Ĥ ,(11)

where diag(z) means a diagonal matrix having a vector z as its diagonal elements. If diag(θ̂) is an identity
matrix (θ̂ = 1), Eq. (11) reduces to the standard probT

lem of matrix factorization given by X̂ = Ŵ Ĥ . The
optimal values of θ̂, Ŵ , and Ĥ depend on what kindsof constraints are placed on those variables.
2.3.1. Formulating a Probabilistic Model
We aim to formulate a Bayesian model of Eq. (10). A
key feature is to consider essential correlations between
M elements of basis ŵ kn . A natural choice is to put amultivariate Gaussian prior on ŵ kn as follows:ŵ kn ∼ N (0, V̂ k ),(12)

where V̂ k ∈ RM ×M is a full covariance matrix. In
general, the Gaussian mean is set to a zero vector. For
example, an audio signal is recorded as real numbers
distributed on both sides of zero (see Figure 2).
The linear relationship x̂nk = θ̂k ĥkn ŵkn and Eq. (12)lead to a likelihood of x̂nk as follows:x̂nk |θ̂, Ŵ , Ĥ ∼ N (0, θ̂k2 ĥ2kn V̂ k ).

(13)Then, using the reproducing property
K of the Gaussian
and the linear relationship x̂n = k=1 x̂nk , we get a

In many applications such as source separation, latentcomponent x̂nk is of main interest. One might consider
it necessary to calculate x̂nk = θ̂k ĥkn ŵ kn . Instead,marginalizing out ŵ kn gives the posterior of x̂nk as a
Gaussian whose mean and covariance areE[x̂nk |x̂n , θ, H, V ] = Y nk Y −1(15)n x̂n ,
V[x̂nk |x̂n , θ, H, V ] = Y nk − Y nk Y −1
Y
,
(16)
nk
n
K
where Y nk = θk hkn V k and Y n = k=1 Y nk are PSD
matrices. For a Bayesian treatment, we need to calculate E[x̂nk |x̂n ] and V[x̂nk |x̂n ] by marginalizing out θ,H, and V under a posterior over these variables, but
this is analytically intractable. One alternative is tosubstitute maximum-a-posteriori (MAP) estimates of
θ, H, and V into Eqs. (15) and (16).
2.4. Fourier Trick
We here discuss the formulation of LD-PSDTF in thefrequency domain. Using Eq. (14), the complex spectrum F x̂n (linear transformation of x̂n ) is found to becomplex-Gaussian distributed as follows: Kθ̂k2 ĥ2kn F V̂ k F H . (17)F x̂n |θ̂, Ŵ , Ĥ ∼ Nc 0,
k=1

It is known that V̂ k can be diagonalized by using F
if V̂ k is strictly a circulant matrix. A trivial example
is a case that V̂ k is a scaled identity matrix, i.e., ŵkn
is stationary white Gaussian noise. If V̂ k is a periodic
kernel and its size M is much larger than its period,V̂ k can be roughly viewed as a circulant matrix.
These facts justify IS-NMF for power-spectrogram decomposition. Since music audio signals roughly consistof pitched sounds and percussive sounds, it is reasonable to approximate V̂ k as a convex combination of

Infinite Positive Semidefinite Tensor Factorization

periodic kernels (for pitched sounds) and identity matrices (for percussive sounds). In the frequency domain
LD-PSDTF thus reduces to IS-NMF discarding thecovariance between frequency bins, while in the time
domain the full covariance structure is still taken into
account. This approximation dramatically reduces thecomputational cost of LD-PSDTF from O(M 3 N K) toO(M N K) as suggested in (Liutkus et al., 2011).

3. Variational Inference
We explain an inference method for a Bayesian modelof GaP-LD-PSDTF deﬁned by Eqs. (4), (5), (6), and(8). Given the observed data X, our goal is to calculate a posterior p(θ, H, V |X) by using the Bayes rulep(θ, H, V |X) = p(X, θ, H, V )/p(X). Since p(X) is
analytically intractable, we use a variational Bayesian(VB) method for approximating p(θ, H, V |X) by afactorizable distribution q(θ, H, V ) given by
N



K
q(hkn ) q(V k ) . (18)q(θk )q(θ, H, V ) =
k=1

n=1

These factors can be alternately updated to monotonically increase a log-evidence lower bound L given bylog p(X) ≥ E[log p(X|θ, H, V )]+ E[log p(θ)] + E[log p(H)] + E[log p(V )]− E[log q(θ)] − E[log q(H)] − E[log q(V )] ≡ L. (19)Since the ﬁrst term is still intractable, we need to takea further lower bound L such that L ≥ L . Note that
L can be indirectly maximized by maximizing L . The
updating formulas areq(θ) ∝ p(θ) exp(Eq(H,V ) [log q(X|θ, H, V )]),
q(H) ∝ p(H) exp(Eq(θ,V ) [log q(X|θ, H, V )]), (20)q(V ) ∝ p(V ) exp(Eq(θ,H) [log q(X|θ, H, V )]),
where log q(X|θ, H, V ) is a variational lower boundof log p(X|θ, H, V ), which is given by Eq. (23).
3.1. Log-Evidence Lower BoundTo derive the tractable bound L , we focus on the convexity and concavity of matrix-variate functions over
PSD matrices. For example, f (V ) = log |V | is concave and g(V ) = tr(ZV −1 ) is convex for any PSD
matrix Z. Let M be the dimension of V .

We then use the following matrix inequality, proposedby Sawada et al. (2012), regarding g(V ):
 
K
−1
K
, (22)tr Φk ZΦTk V −1≤
tr Zkk=1 V kk=1

where {V k }K
k=1 is a set of arbitrary PSD matrices,matrices that sum to the{Φk }K
k=1 is a set of auxiliaryΦk = I), and the equality isidentity matrix (i.e., k satisﬁed when Φk = V k ( k V k )−1 .

Using Eqs. (21) and (22), we can derive the tractable
lower bound of E[log p(X|θ, H, V )] (the ﬁrst term ofEq. (19)). A term regarding X n is bounded as follows:E[log p(X n |θ, H, V )] (see Eq. (7))
(23)

ν  
ν
−1
+ const.
= − E[log |Y n |] − E tr X n Y n
2
2
 νM
ν
ν  
−1≥ − log |Ωn | −+k E tr θk hkn V k Ωn
2
2
2


ν
T
−1 −1 −1
E tr θk hkn V k Φnk X n Φnk + const.,−
2 k
where Ωn is a PSD matrix and {Φnk }K
k=1 is a set of
auxiliary matrices that sum to an identity matrix. Letting the partial derivatives of Eq. (23) equal to be zero,
we can obtain the optimal values of Ωn and {Φnk }K
k=1
that satisfy the equality as follows:Ωn = k E[Y nk ],(24)
  −1 −1    −1 −1 −1. (25)
Φnk = E Y nkk E Y nk
3.2. Variational Bayesian Update
Here we discuss the functional forms of q(θk ), q(hkn ),and q(V k ). A problem lies in the non-conjugacy of theBayesian model. Eq. (23) involves the expectations
both of the scalar variables and of their reciprocals.
This is, the suﬃcient statistics are x and x−1 , althoughthose of the gamma prior are log(x) and x. This meansthat the functional forms of q(θk ) and q(hkn ) are givenby the generalized inverse Gaussian (GIG) distribution, as shown in Hoﬀman et al. (2010). Note that the
GIG distribution is deﬁned asγ

GIG(x|γ, ρ, τ ) =−1
1(ρ/τ ) 2
xγ−1 e− 2 (ρx+τ x ) , (26)√2Kγ ( ρτ )(21)

where γ, ρ > 0, and τ > 0 are parameters and Kγ isthe modiﬁed Bessel function of the second kind. The
expectations E[x] and E[x−1 ] are given by
√
  √
√
√
τ Kγ+1 ( ρτ )
ρKγ−1 ( ρτ )1√
, E
. (27)
E[x] = √=√√ρKγ ( ρτ )x
τ Kγ ( ρτ )where Ω is an arbitrary PSD matrix (tangent point)
and the equality is satisﬁed when Ω = V .

As to matrix variable V k , we found that the functional
form of q(V k ) is given by the matrix GIG (MGIG) dis-

We ﬁrst calculate a tangent plane of f (V ) by using aﬁrst-order Taylor expansion as follows:
log |V | ≤ log |Ω| + tr(Ω−1 V ) − M,

Infinite Positive Semidefinite Tensor Factorization

tribution (Barndorﬀ-Nielsen et al., 1982). The MGIG
distribution over PSD matrix X is deﬁned as
M +12γM|X|γ− 2
MGIG(X|γ, R, T ) =γ|T | Bγ (RT /4)1 
exp − tr RX + T X −1 , (28)2where γ is a real number, R, T > 0 are PSD matrices,
M is the size of X, and Bγ is the matrix Bessel function of the second kind (Herz, 1955). It includes the
Wishart distribution as a special case (Butler, 1998)
and its suﬃcient statistics are log |X|, X, and X −1 .
To calculate E[X] and E[X −1 ], we use a Monte Carlo
method as described in the supplementary material.
Consequently, we can assume the following forms:q(θk ) = GIG(θk |γkθ , ρθk , τkθ ),hhq(hkn ) = GIG(hkn |γkn
, ρhkn , τkn),
VVVq(V k ) = MGIG(V k |γk , Rk , T k ).

(29)

Pk =

N


hkn Y −1
n , Qk =

n=1

N


−1
hkn Y −1n X n Y n . (33)n=1

Sawada et al. (2012) derived a complicated solution of
Eq. (32), but we can solve it analytically by using theCholesky decomposition Qk = Lk LTk , where Lk is alower triangular matrix. Finally, we get1

V k ← V k Lk (LTk V k P k V k Lk )− 2 LTk V k .

(34)

When all the matrices are diagonal, Eqs. (31) and (34)
reduce to the one in IS-NMF (Nakano et al., 2010).

4. Related WorkWe show that PSDTF has deep connections to nonnegative matrix factorization (NMF), tensor factorization(TF), and principal component analysis (PCA).
4.1. Nonnegative Matrix Factorization

These parameters are iteratively updated as follows: γkθ = αc/K, ρθk = 2α + ν n tr E[hkn V k ] Ω−1
,
n



−1
Φnk X n ΦTnk ,τkθ = ν n tr E h−1
kn V kh,= a0 , ρhkn = 2b0 + νtr E[θk V k ] Ω−1
γkn
n


h
(30)
Φnk X n ΦTnk ,τkn= νtr E θk−1 V −1k−1γkV = ν0 /2, RVk = V −1
0 +νn E[θk hkn ] Ωn ,TT Vk = ν n E θk−1 h−1
kn Φnk X n Φnk .
3.3. Multiplicative Update
The multiplicative update (MU) is a well-known optimization technique often used for maximum-likelihood
estimation of NMF. To show a clear connection of LDPSDTF to IS-NMF, we derive closed-form MU rules
for calculating the point estimates of H and V . Notethat we assume θk = 1 and tr(V k ) = 1 (unit trace)
to remove the scale arbitrariness. If tr(V k ) = s, thescale adjustments V k ← 1s V k and hk ← shk do notchange the LD divergence DLD (X n |Y n ).
We aim to maximize the log-likelihood given by removing the expectation operators from Eq. (23). Letting
the partial derivative with respect to hkn equal to bezero, we get the following update rule: −1
tr Y −1n V kY n Xn


.
(31)
hkn ← hkn
tr Y −1n VkThen, letting the partial derivative with respect to V kequal to be zero, we get the following equation:old
V k P k V k = V old
k Qk V k ,where P k and Qk are PSD matrices given by(32)

PSDTF includes NMF as a special case. If we restrict
PSD matrices X n and V k to diagonal matrices (i.e.,
X n = diag(xn ) and V k = diag(v k ) for some nonnegative vectors xn and v k ) Eq. (1) can be written as
xn ≈

K


θk hkn v k ,

(35)

k=1

where θk and hkn are nonnegative numbers. If θk = 1,
this model reduces to the basic formulation of NMF.
Févotte et al. (2009) showed that the IS divergence istheoretically suitable for evaluating the reconstruction
error of Eq. (35) in the task of audio source separation.
Hoﬀman et al. (2010) proposed an inﬁnite extension of
IS-NMF (GaP-IS-NMF) using a gamma process prioron θ. GaP-LD-PSDTF can therefore be viewed as anatural extension of GaP-IS-NMF.
An interesting ﬁnding is that the positive semideﬁniteness (nonnegative deﬁniteness) constraint on matrices
in PSDTF induces sparse decomposition like the nonnegativity constraint on vectors and scalars as in NMF.
Positive semideﬁniteness can therefore be considered ageneralization of the nonnegativity concept.
4.2. Tensor Factorization
PSDTF is related to a variant of TF called canonicalpolyadic (CP) decomposition (Carroll & Chang, 1970;
Harshman, 1970). If we restrict a PSD matrix V k to arank-1 matrix (i.e., V k = uk uTk for some vector uk ),Eq. (1) can be written as
X≈

K

k=1

θk hk ⊗ uk ⊗ uk .

(36)

Infinite Positive Semidefinite Tensor Factorization

This can be viewed as CP decomposition in which basis
vectors of the second mode, {uk }K
k=1 , are constrained
to be equal to those of the third mode, and θk and hk
should be nonnegative. In addition, PSDTF uses the
LD divergence for evaluating the reconstruction error
while typical TF uses the Euclidean distance.
There are some other related models. Tucker decomposition (Tucker, 1966) is a generalization of CP decomposition and Xu et al. (2012) proposed its inﬁnite
extension (K → ∞) based on the Gaussian or t process. Shashua & Hazan (2005) proposed nonnegative
TF (NTF) that, like NMF, imposes a nonnegativity
constraint on all elements of factors. In this paper thenonnegativity constraint on θk and hk , not on uk , led
to a new class of TF.

LD-PSDTF is related to a major class of matrix factorization (MF) using the Gaussian distribution as a corebuilding block of probabilistic models. Let us recall
the MF model given by Eq. (11):T

Multiplicative update (Maximum-likelihood estimation for LD-PSDTF)

Variational Bayesian update (Bayesian inference for GaP-LD-PSDTF)

Figure 3. Experimental results for synthetic data.
0.24

0.016

-0.065-0.066

-0.46

Front
Back

4.3. Principal Component AnalysisX̂ = Ŵ Ĥ

True positive semidefinite basis matrices (K=6)(37)

where X̂ ∈ RM ×N , Ŵ ∈ RM ×K , and Ĥ ∈ RN ×K (Nis the number of observations).
Several probabilistic models are obtained by putting
Gaussian priors in diﬀerent ways. Placing an isotropic
Gaussian prior on N columns of Ĥ (latent-space coordinates corresponding to the observations) leads toprobabilistic PCA (PPCA) (Bishop, 1999). If we put
an isotropic Gaussian prior on M rows of Ŵ (mapping
functions from the latent space to the observed space),
the resulting model is called dual PPCA. Marginalizing Ŵ out, we can formulate a Gaussian process latent variable model (GPLVM) (Lawrence, 2003). PSDTF, on the other hand, puts a full-covariance Gaussian prior on K columns of Ŵ . If we instead use a GP
prior, PSDTF given by Eq. (14) can be viewed as multiple kernel learning (MKL) (Lanckriet et al., 2004).

5. Evaluation
This section reports experiments to evaluate the performance of LD-PSDTF.
5.1. Synthetic Data
We evaluated the capability of GaP-LD-PSDTF to dis+cover basis PSD matrices {V k }K
k=1 used for generating
an observed tensor X and to estimate the value of K + .
In this experiment we use M = ν = 10, N = 2000, andK + = 6. The synthetic data X was stochastically gen-

Figure 4. Brain activities discovered from EEG data.

erated according to the following process:hkn ∼ G(0.1, 0.1),
V k ∼ W(10, I/10),
νX n ∼ W (ν, k hkn V k ) ,(38)

To learn the GaP-LD-PSDTF model, we used the VB
algorithm with a truncation level K = 100, and hyperparameters α = c = 1, a0 = b0 = 0.1, ν0 = 10, andV 0 = I/ν0 . Since Monte Carlo simulation of E[V k ]and E[V −1k ] was found to be often unreliable, we instead calculated the maximum-a-posteriori estimates.
For maximum-likelihood estimation of the LD-PSDTF
model, we used the MU algorithm with K = 6. The
both methods were initialized randomly.
As shown in Figure 3, the experimental results showedthat the both models successfully discovered the correct basis matrices. The hyperparameters were notsensitive to the results. In GaP-LD-PSDTF, the true
number of bases, K + = 6, was correctly estimated.
5.2. EEG Data
We then tested LD-PSDTF on a popular EEG dataset(Blankertz, 2001). We aimed to predict a left or righthand movement (label -1 or 1) from 500 ms EEG signals recorded at 28 channels of the brain (M = 28)
with a sampling rate of 100 Hz. There were 416 trials(N = 416) of which 100 belong to the test set. For eachtrial we calculated a full-rank covariance matrix over
50 frames. The PSD basis matrices and their activations were estimated by using the MU algorithm with
K = 5 in an unsupervised manner. We used Fisher’s
LDA for binary classiﬁcation of the K-dimensional activations corresponding to the individual trials.

Infinite Positive Semidefinite Tensor Factorization

Note that the best results from the competition wereobtained by combining the ﬁrst-order features (lateralized readiness potential: LRP) and the second-orderones (ERD), whereas our method used only the latter.
Finding discriminative patterns without any supervision in this context is a highly nontrivial task, because
discriminative signals are much weaker than irrelevant
oscillatory activities (e.g., occipital alpha waves).
5.3. Music Data
We evaluated LD-PSDTF for music signal analysis.
As discussed in Section 2.4, LD-PSDTF formulated in
the time domain is equivalent to IS-NMF formulated
in the frequency domain if {V k }K
k=1 are circulant matrices such as periodic kernels, identity matrices, andtheir conic sums. Since this is a reasonable assumption
for music signals, we used the Fourier trick explained in
Section 2.4. We tested an inﬁnite model with a truncation level K = 100 and ﬁnite models (θk = 1) with
diﬀerent K ranging from 1 to 100. For comparison, we
tested an inﬁnite model of KL-NMF (GaP-KL-NMF)
and ﬁnite models with diﬀerent K ranging from 1 to
100 for amplitude-spectrogram decomposition.
We used three songs (No.1, 2, and 3) from the “RWC
Music Database: Popular Music” (Goto et al., 2002).
The CD-quality audio signals were downsampled at 16[kHz] and were analyzed by using short-time Fouriertransform with a window size of 128 [ms] and a shifting
interval of 64 [ms]. The size of X is speciﬁed as M =2048, and N = 3237, 3447, and 3020, respectively. We
used the VB algorithm with hyperparameters α = 1,
c = 1, a0 = b0 = 0.1, ν0 = M , ν = 1, and V 0 = I/ν0 .
The experimental results showed that in each song thevalue of K + chosen by GaP-LD-PSDTF was close to
the best value of K found by ﬁnite models (Figure 5).
This proves that GaP-LD-PSDTF has an ability of automatic model-order selection without expensive gridsearch. Similar results were obtained in KL-NMF, but
GaP-LD-PSDTF achieved much higher log-evidence
lower bounds than GaP-KL-NMF did. This supports
the appropriateness of IS-NMF for music analysis. Additional results of source separation with sound samples are given in the supplementary material.

x10

Log-evidence lower bounds

The most signiﬁcant principal component of each basis matrix is shown in Figure 4, in which each number
indicates the correlation between the ground-truth labels and the estimated activations on the test set. The
accuracies of classiﬁcation were 73% (K = 5) and 79%(K = 10). The leftmost and rightmost bases with highcorrelations are compatible with the well-known physiological process called event-related desynchronization(ERD). We consider these results promising.

7

LD-PSDTF (IS-NMF)
x10
-3.2

K+ = 636

2.9

GaP model

2.8

Finite models

-3.4

50

-3.6
K
100
1

2.7

1

K+ = 61

3.8-3.1

K+ = 40K
100

50
K+ = 45

-3.5

3.6
3.4

KL-NMF

1

50

K
-3.9
100
1
K+ = 66

3.1
3.0

K
100

50

-2.8

K+ = 61

-3.2

2.9
1

50

K
-3.6
100
1

50

K
100

Figure 5. Experimental results for music recordings.

6. Conclusion
This paper presented positive semideﬁnite tensor factorization (PSDTF) as a natural extension of nonnegative matrix factorization (NMF). We used a Bregman
matrix divergence called the log-determinant (LD) divergence as the reconstruction error. This LD-PSDTF
can be viewed as a natural extension of NMF basedon the Itakura-Saito (IS) divergence. We formulated a
nonparametric Bayesian model that allows an observed
tensor to contain an unbounded number of bases andderived a variational Bayesian algorithm and a multiplicative update rule by using matrix inequalities. In
addition, we showed the eﬀectiveness of the Fourier
trick, i.e., the frequency-domain formulation can dramatically reduce the computational cost in some applications such as music signal analysis.
One interesting open question is what kind of PSDTF
can be viewed as an extension of NMF based on the
Kullback-Leibler (KL) divergence. The von-Neumann(vN) divergence (Tsuda et al., 2005) is well known as
another major type of the Bregman matrix divergencethat includes the KL divergence as a special case. Substituting φ(Z) = tr(Z log Z − Z) into Eq. (2), the vN
divergence is given by
DvN (X n |Y n )= tr (X n log X n − X n log Y n − X n + Y n ) . (39)The assumption underlying p(X n |Y n ), however, is not
obvious although the Bregman divergence must correspond one-to-one to an exponential family. We plan toinvestigate vN-PSDTF to formulate a Bayesian model.
Acknowledgment: This study was partially supported by JSPS KAKENHI 23700184, MEXT KAKENHI 25870192, and JST OngaCREST project. We
thank the reviewers for giving insightful comments.

Infinite Positive Semidefinite Tensor Factorization

References
Barndorﬀ-Nielsen, O., Blæsild, P., Jensen, J. L., andJørgensen, B. Exponential transformation models.
Royal Society of London, 379(1776):41–65, 1982.
Bishop, C. M. Variational principal components. In
ICANN, pp. 509–514, 1999.
Blankertz, B., Curio, G., and Müller, K.-L. Classifying
single trial EEG: Towards brain computer interfacing, In: NIPS, 2001.
www.bbci.de/competition/ii/berlin desc.html
Bregman, L. M. The relaxation method of ﬁnding thecommon points of convex sets and its application
to the solution of problems in convex programming.
USSR CMMP, 7(3):200–217, 1967.
Butler, R. W. Generalized inverse Gaussian distributions and their Wishart connections. Scandinavian
Journal of Statistics, 25(1):69–75, 1998.
Butler, R. W. and Wood, A. Laplace approximationfor Bessel functions of matrix argument. J. of Computational and Applied Math., 155(2):359–382, 2003.
Carroll, J. D. and Chang, J. J. Analysis of individual diﬀerences in multidimensional scaling via an
N-way generalization of ‘Eckart-Young’ decomposition. Psychometrika, 35(3):283–319, 1970.
Cemgil, A. T. Bayesian inference for nonnegative
matrix factorisation models. Computational Intelligence and Neuroscience, Article ID 785152, 2009.
Févotte, C., Bertin, N., and Durrieu, J.-L. Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis. Neural Computation, 21(3):793–830, 2009.
Goto, M., Hashiguchi, H., Nishimura, T., and Oka, R.
RWC music database: Popular, classical, and jazzmusic database. In ISMIR, pp. 287–288, 2002.
Harshman, R. A. Foundations of the PARAFAC procedure: Models and conditions for an “explanatory”
multi-modal factor analysis. UCLA Working Papers
in Phonetics, 16(1), 1970.
Herz, C. S. Bessel functions of matrix argument. Annals of Mathematics, 61(3):474–523, 1955.
Hoﬀman, M., Blei, D., and Cook, P. Bayesian nonparametric matrix factorization for recorded music.
In ICML, pp. 439–446, 2010.
Itakura, F. and Saito, S. Analysis synthesis telephonybased on the maximum likelihood method. In ICA,
pp. C17–C20, 1968.
Kulis, B., Sustik, M., and Dhillon, I. Low-rank kernel

learning with Bregman matrix divergences. JMLR,
10:341–376, 2009.
Kullback, S. and Leibler, R. On information and sufﬁciency. Annals of Math. Stat., 22(1):79–86, 1951.
Lanckriet, G., Cristianini, N., Bartlett, P., Ghaoui,
L., and Jordan, M. Learning the kernel matrix withsemideﬁnite programming. JMLR, 5:27–72, 2004.
Lawrence, N. D. Gaussian process latent variable models for visualisation of high dimensional data. In
NIPS, 2003.
Lee, D. and Seung, H. Algorithms for non-negative
matrix factorization. In NIPS, pp. 556–562, 2000.
Lee, H., Cichocki, A., and Choi, S. Nonnegative matrix
factorization for motor imagery EEG classiﬁcation.
In ICANN, pp. 250–259, 2006.
Liutkus, A., Badeau, R., and Richard, G. Gaussian processes for underdetermined source separation. IEEE Trans. on ASLP, 59(7):3155–3167, 2011.
Nakano, M., Kameoka, H., Roux, J. Le, Kitano, Y.,
Ono, N., and Sagayama, S. Convergence-guaranteed
multiplicative algorithms for non-negative matrix
factorization with beta divergence. In MLSP, pp.
283–288, 2010.
Salakhutdinov, R. and Mnih, A. Bayesian probabilistic matrix factorization using Markov chain Monte
Carlo. In ICML, pp. 880–887, 2008.
Sawada, H., Kameoka, H., Araki, S., and Ueda, N.
Eﬃcient algorithms for multichannel extensions of
Itakura-Saito nonnegative matrix factorization. In
ICASSP, pp. 261–264, 2012.
Shashua, A. and Hazan, T. Non-negative tensor factorization with applications to statistics and computer
vision. In ICML, pp. 792–799, 2005.
Sivalingam, R., Boley, D., Morellas, V., and Papanikolopoulos, N. Tensor sparse coding for region
covariances. In ECCV, pp. 722–735, 2010.
Smaragdis, P. and Brown, J. C. Non-negative matrix
factorization for polyphonic music transcription. In
WASPAA, pp. 177–180, 2003.
Tsuda, K., Rätsch, G., and Warmuth, M. K. Matrix
exponentiated gradient updates for on-line learning
and Bregman projection. JMLR, 6:995–1018, 2005.
Tucker, L. R. Some mathematical notes on three-mode
factor analysis. Psychometrika, 31(3):279–311, 1966.
Xu, Z., Yan, F., and Qi, Y. Inﬁnite Tucker decomposition: Nonparametric Bayesian models for multiway
data analysis. In ICML, pp. 1023–1030, 2012.