2014 IEEE INTERNATIONAL WORKSHOP ON MACHINE LEARNING FOR SIGNAL PROCESSING, SEPT. 21–24, 2014, REIMS, FRANCETRAINING RESTRICTED BOLTZMANN MACHINES WITH
AUXILIARY FUNCTION APPROACH
Hirokazu Kameoka1),2) and Norihiro Takamune1)
1)2)

Graduate School of Information Science and Technology, The University of TokyoNTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation

ABSTRACT
Restricted Boltzmann Machines (RBMs) are neural networkmodels for unsupervised learning, but have recently found awide range of applications as feature extractors for supervised
learning algorithms. They have also received a lot of attention
recently after being proposed as building blocks for deep belief networks. The success of these models raises the issueof how best to train them. At present, the most popular training algorithm for RBMs is the Contrastive Divergence (CD)learning algorithm. The aim of this paper is to seek for a new
optimization algorithm tailored for training RBMs in the hopeof obtaining a faster algorithm than the CD algorithm. We
propose deriving a new training algorithm for RBMs basedon an auxiliary function approach. Through an experiment onparameter training of an RBM, we conﬁrmed that the present
algorithm converged faster and to a better solution than the
CD algorithm.
Index Terms— Restricted Boltzmann machines, deep belief networks, contrastive divergence learning algorithm, auxiliary function, minorization-maximization
1. INTRODUCTION
Restricted Boltzmann machines (RBMs) [1] are probabilistic
generative models that can be interpreted as stochastic neural networks. They have found applications in dimensionality
reduction, classiﬁcation, feature learning, etc., and have recently attracted considerable attention after being proposed asbuilding blocks for deep belief networks (DBNs)[2]. Over thelast few years, DBNs have been applied with notable successto a wide range of applications including image recognitionand speech recognition [3].
A DBN can be constructed by stacking multiple RBMs
so that the hidden layer of one RBM becomes the visible
layer of another RBM. DBNs can be efﬁciently trained using greedy layerwise training, in which the RBMs are trained
one at a time in a bottom-up fashion. Furthermore, a discriminative neural network for classiﬁcation or regression can beobtained by adding an output layer on top of the trained RBM
or DBN where the additional units represent labels (e.g., classlabels) corresponding to the observed data. The initialization
obtained with the layer-wise pretraining of RBMs facilitatesc
978-1-4799-3694-6/14/$31.00 ⃝2014IEEE

the optimization problem of training the discriminative neuralnetwork, a procedure known as the global ﬁne-tuning.
As their name implies, RBMs are a special case of Boltzmann machines. As with Boltzmann machines, RBMs consistof visible units and hidden units. The visible units constitutethe ﬁrst layer, which correspond to the components of an observation. The hidden units constitute the second layer, whichmodel dependencies between the components of observations
and can be interpreted as non-linear feature detectors. In the
RBMs network graph, each unit is connected to all the units
in the other layer. However, there are no connections between
units in the same layer. This restriction allows for more efﬁcient training algorithms than those applicable for the generalclass of Boltzmann machines, in particular the gradient-based
contrastive divergence (CD) algorithm [4]. The aim of thispaper is to seek for a new optimization algorithm tailored fortraining RBMs in the hope of obtaining a faster algorithm thanthe CD algorithm.
For many nonlinear optimization problems, parameter estimation algorithms constructed using an auxiliary function
have proven to be very effective. The general principle forthe parameter estimation scheme using an auxiliary function
is referred to as the “auxiliary function approach” (or alternatively the “minorization-maximiation (MM) approach” [5]).
Note that the auxiliary function approach itself is not an algorithm, but a description of how to construct an optimization
algorithm. When applying the auxiliary function approach toa certain optimization problem, the ﬁrst step is to design an
auxiliary function that upper-bounds or lower-bounds the objective function. A parameter estimation algorithm can thenbe derived using the auxiliary function. An algorithm thatconsists of iteratively minimizing/maximizing the auxiliary
function is guaranteed to converge to a stationary point of
the objective function. It should be noted that this concept
is adopted in many existing algorithms. For example, theexpectation-maximization (EM) algorithm [6] builds a surrogate for a likelihood function of latent variable models byusing Jensen’s inequality. It is also well known for its use indevising an algorithm for non-negative matrix factorization[7, 8]. In general, if we can build a tight upper/lower boundfunction for the objective function of a speciﬁc optimization
problem, we expect to obtain a fast-converging algorithm.
In fact, the authors and colleagues have thus far proposed
deriving parameter estimation algorithms based on the aux-

For a model of the form (5) with parameters Θ, the mean
log-likelihood given N training examples v (1) , . . . , v (N ) isthus given byL(Θ) =

1 ∑log p(v (n) |Θ).
N n

(6)

Fig. 1. Graph structure of RBM.
Maximum likelihood training of an RBM can be formulated
as the problem of maximizing L(Θ) with respect to Θ.
iliary function approach for various optimization problems,
some of which have been proven to be signiﬁcantly faster thangradient-based methods (e.g., [8–15]). Motivated by these experiences, we propose deriving a new training algorithm for
RBMs based on the auxiliary function approach.
2. RESTRICTED BOLTZMANN MACHINES
2.1. Likelihood function of RBMs
RBMs are undirected generative models that use a layer ofhidden variables to model a distribution over visible variables.
As depicted in Fig. 1, RBMs can be viewed as a variant of
Boltzmann Machines, with the restriction that their neurons
must form a bipartite graph structure. This restriction gives
the RBM its name.
The standard type of RBM has binary-valued hidden andvisible units. A binary RBM with J hidden units is a parametric model of the joint distribution between a layer of hidden
variables h = (h1 , . . . , hJ )T ∈ {0, 1}J and the visible variables v = (v1 , . . . , vI )T ∈ {0, 1}I , that is given byexp(−E(v, h; Θ))p(v, h|Θ) =,Z(Θ)

E(v, h; Θ) = −aT v − bT h − v T W h
(2)
∑
∑
∑
ai vi −
bj hj −
vi wi,j hj , (3)
=−
j

i,j

with parameters Θ = {W , a, b}. W = (wi,j )I×J is aweight matrix whose element wi,j is associated with the connection between hidden unit hj and visible unit vi , and a =(a1 , . . . , aI )T and b = (b1 , . . . , bJ )T are bias weight vectors(offsets) associated with the visible units and the hidden units,
respectively. Z is called the “partition function” and is givenby the sum of exp(−E(v, h; Θ)) over all possible pairs ofvisible and hidden vectors:∑∑
Z(Θ) =e−E(v,h;Θ) .
(4)
v

h

The marginal distribution of v is given by summing over all
possible hidden vectors:p(v|Θ) =∑h

1 ∑ −E(v,h;Θ)p(v, h|Θ) =
e
.
Z(Θ)
h

Since it is not possible to ﬁnd the parameters maximizing themean log-likelihood analytically, one convenient way would
be to apply a gradient ascent approach. Although here we
consider only the update rule for wi,j , the update rule for thebiases ai , bj can be deﬁned analogously.
The partial derivative of the mean log-likelihood with respect to wi,j is given by
∂L(Θ)
1 ∑∑
(n)
=−
p(h|v (n) , Θ)vi hj
∂wi,j
N n
h
∑∑
+
p(v, h|Θ)vi hj .
v

(7)

h

This leads to a simple update rule based on steepest ascent in
the log-likelihood:wi,j ← wi,j + ϵ

∂L(Θ)
,
∂wi,j

(8)

(1)

where E(v, h; Θ) denotes a bilinear energy function given byi2.2. Contrastive divergence learning algorithm [2, 4]

where ϵ is a learning rate. As can be seen from (7), the gradient ∂L(Θ)/∂wi,j is given by the difference between two
expectations: the expectation of the energy gradient under theRBM distribution p(v, h|Θ) and under the conditional distribution of the hidden variables p(h|v, Θ) given a set of the
training samples. To perform gradient ascent, we need tocompute these two terms.
Since there are no direct connections between hidden
units as well as between visible units in an RBM, hidden
units are conditionally independent given the state of thevisible variables and vice versa:∏p(h|v, Θ) =p(hj |v, Θ),(9)jp(v|h, Θ) =∏p(vi |h, Θ).

The conditional distributions p(hj |v, Θ) and p(vi |h, Θ) canbe interpreted as the ﬁring rate of each neuron, given by()∑
p(hj = 1|v, Θ) = σ bi +vi wi,j ,
(11)
(

(5)

(10)

i

p(vi = 1|h, Θ) = σ ai +

i

∑
j

)
wi,j hj ,(12)where σ denotes the logistic sigmoid: σ(x) = 1/(1 + e−x ).
By using these facts, the ﬁrst term of (7) can be rewritten as
1 ∑(n)p(hi = 1|v (n) , Θ)vj ,
N n

(13)

which is easy to compute. As for the second term, however,directly calculating the sums that run over all values of v andh leads to a computational complexity which is in generalexponential in the number of variables. To avoid this exponential complexity, one convenient way is to approximate theexpectation by samples from the RBM distribution. Thesesamples can be obtained with Gibbs sampling using the conditional distributions (11) and (12). It has been shown that
estimates obtained after running the Gibbs chain for just a
few steps can be sufﬁcient in practice.
This way of updating the parameters belongs to the family of the contrastive divergence (CD) learning algorithm [4],which has become a standard way of training RBMs.

and h, and the sums inside the exponential function that run
over i and j. Since it is preferable to use a function that canbe maximized analytically as an auxiliary function, we would
like to design an auxiliary function such that the parameters
are separated into individual terms.
First, let us focus on the sum over h inside the logarithm
function in L(Θ). We can invoke the following inequality toconstruct a lower bound function of L(Θ).
Lemma 1 (Jensen’s inequality for concave functions with
non-negative arguments). For any concave function f , any
non-negative values x1 , . . . , xK ∈ [0, ∞) and any positive∑ weights λ1 , . . . , λK ∈ (0, 1) that sums to unity (i.e.,k λk = 1), we have
f

(∑
k

xk
.
m xm

λk = ∑

3.1. Auxiliary function approach
Here we introduce the general principle of the auxiliary function approach (a.k.a the MM approach).
Let us use G(θ) to denote an objective function that wewant to maximize with respect to θ. G+ (θ, θ̃) is deﬁned as an
auxiliary function for G(θ) if it satisﬁes(14)

Theorem 1. G(θ) is non-decreasing under the updates, θ ←argmaxθ G+ (θ, θ̃) and θ̃ ← argmaxθ̃ G+ (θ, θ̃).
Proof: Let us set θ at an arbitrary value θℓ and deﬁne θ̃ℓ+1 =argmaxθ̃ G+ (θℓ , θ̃) and θℓ+1 = argmaxθ G+ (θ, θ̃ℓ+1 ). First,
it is obvious that G(θℓ ) = G+ (θℓ , θ̃ℓ+1 ). Next, we can conﬁrm that G+ (θℓ , θ̃ℓ+1 ) ≤ G+ (θℓ+1 , θ̃ℓ+1 ) since θℓ+1 is the
maximizer of G+ (θ, θ̃ℓ+1 ) with respect to θ. By deﬁnition,
it is obvious that G+ (θℓ+1 , θ̃ℓ+1 ) ≤ G(θℓ+1 ) and so we canﬁnally show that G(θℓ ) ≤ G(θℓ+1 ).
We propose applying this principle for the problem ofmaximizing the mean log-likelihood L(Θ) of RBM.
3.2. Designing auxiliary function
When applying the auxiliary function approach to a certain
maximization problem, the ﬁrst step is to design an auxiliary function that lower-bounds the objective function. The
difﬁculty in solving the optimization problem of maximizing
L(Θ) lies in the nonlinear interdependence of the parameters
due to the sums inside the logarithm function that run over v(16)

Since the logarithm function is a concave function, and ofcourse the value of p(v (n) , h|Θ) is always non-negative, we
can use the above inequality:
L(Θ) =∑1 ∑log
p(v (n) , h|Θ)
N n
h

≥

1 ∑∑
N

n

θ̃

We call θ̃ an auxiliary variable. By using G+ (θ, θ̃), G(θ) canbe iteratively increased according to the following theorem:

(15)

k

where equality holds if and only if

3. AUXILIARY FUNCTION APPROACH FOR
TRAINING RBM

G(θ) = max G+ (θ, θ̃).

(x )
) ∑
k
λk f
xk ≥
,
λk

λn (h) loghp(v (n) , h|Θ),λn (h)(17)

where λn (h) is a positive weight that sums to unity, such that∑h λn (h) = 1. Here, equality holds whenp(v (n) , h|Θ)= p(h|v (n) , Θ).
(n) , h′ |Θ)h′ p(vλn (h) = ∑(18)

Thus, (17) can be rewritten as
1 ∑∑λn (h)E(v (n) , h; Θ)N nh∑− log Z(Θ) −λn (h) log λn (h).

L(Θ) ≥ −(19)h

It is still impossible to obtain a closed-form update rule for Θ
from this lower bound function because of the second term.
In the second term, we notice that the sum is taken over v
and h inside the negative logarithm function. We can use thefollowing inequality to construct a looser lower bound.
Lemma 2. For any differentiable convex function f and any
real numbers x and ζ, we havef (x) ≥ f (ζ) + f ′ (ζ)(x − ζ),
where equality holds if and only if x = ζ.

(20)

αjb (v, h) = bj hj + βjb E(v, h; Θ),

Since the negative logarithm function is a differentiable
convex function, we can apply this inequality to − log Z(Θ):
1− log Z(Θ) ≥ − log ζ − (Z(Θ) − ζ),ζ(21)where ζ is a positive number. Here, equality holds when
ζ = Z(Θ).

(22)

A lower bound function can thus be constructed by replacing − log Z(Θ) with the right-hand side of (21). However, it
is still impossible to obtain a closed-form update rule for Θ
from this lower bound function because of the sums that run
over i and j inside the negative exponential function in Z(Θ).
Since the negative exponential function is a concave function,
it may appear that Lemma 1 can be applied in a similar way
that we obtain the ﬁrst inequality. However, since the argument of the negative exponential function, i.e., −E(v, h; Θ),can possibly take on negative values, the inequality given in
Lemma 1 cannot be applied in this case. Instead, here we introduce the following inequality (which we have previously
employed for a different optimization problem [11]).
Lemma 3 (Jensen’s inequality for concave functions with realarguments). For any concave function f , any real numbers
x1 , . . . , xK ∈ R, any positive
weights β1 , . . . , βK ∈ (0, 1)∑that sums to unity (i.e., k βk = 1), and∑ any real numbers
α1 , . . . , αK ∈ R that sums to zero (i.e., k αk = 0), we have(x − α )(∑ ) ∑
k
k
βk f
xk ≥
f
,
(23)
βk
k

k

where equality holds if and only if∑
αk = xk − βk
xm .

(24)

m

By using this inequality, we obtain
(∑
)
∑
∑
− exp
ai vi +
b j hj +vi wi,j hj
i

≥−

∑

βia exp

i

−

∑

(
(

βjb exp

j

−

∑
i,j

i,j
)
ai vi − αia (v, h)
βia

(
w
βi,j
exp

)

1 ∑∑
λn (h)E(v (n) , h; Θ)
N n
h
1 ∑∑
−
G(v, h; Θ)ζ vh∑− log ζ + 1 −λn (h) log λn (h),L+ (Θ, Θ̃) = −

(28)

h

where
G(v, h; Θ) =∑(βia expi+∑(
βjb expj+∑ai vi − αia (v, h)βia

bj hj − αjb (v, h)

(
w
exp
βi,j

i,j

)
)

βjb

)
w
vi wi,j hj − αi,j(v, h), (29)
w
βi,jwith Θ̃ = {λ, ζ, α}, λ = {λn (h)|n = 1, . . . , N, h ∈ {0, 1}J }w
and α = {αia (v, h), αjb (v, h), αi,j(v, h)|i = 1, . . . , I, j =IJ1, . . . , J, v ∈ {0, 1} , h ∈ {0, 1} }. Then, L+ (Θ, Θ̃) is an
auxiliary function for L(Θ).

Now we can derive the update equations for the model parameters by using the auxiliary function deﬁned above. As statedin Theorem 1, we can iteratively increase the value of L(Θ)
by performing the updates
(30)

Θ̃

Θ ← argmax L+ (Θ, Θ̃).

(31)

Θ

(26)

w
where αia (v, h), αjb (v, h) and αi,j(v, h) are arbitrary real
∑ a∑numbers that sum to zero, i.e., i αi (v, h) + j αjb (v, h) +∑ab
β w are arbitraryi,j αi,j (v, h) = 0, and βi,j , βi,j and∑ i,ja ∑ bpositive numbers that sum to one, i.e., i βi + j βj +∑wi,j βi,j = 1. Here, equality holds when

αia (v, h) = ai vi + βia E(v, h; Θ),Theorem 2 (Auxiliary function for RBMs). Deﬁne L+ asΘ̃ ← argmax L+ (Θ, Θ̃),

βjb

)
w
vi wi,j hj − αi,j(v, h),w
βi,j

We can obtain a slightly looser lower bound function forL(Θ) by further replacing −e−E(v,h;Θ) with the right-handside of (26). It is important to note that this lower bound function is given in such a way that the parameters are separated
in individual terms, thus allowing us to derive closed-form
update equations for the parameters.
To sum up, we obtain the following theorem:

3.3. Update equations(25)

j

bj hj − αjb (v, h)(27)w
αiw (v, h) = vi wi,j hj + βi,jE(v, h; Θ).

(30) is given explicitly as (18), (22) and (27), respectively.
(31) can be obtained by solving ∂L+ /∂ai = 0, ∂L+ /∂bj = 0
and ∂L+ /∂wi,j = 0. For example, ∂L+ /∂wi,j is given by
1 ∑∑∂L+ (Θ)(n)=λn (h)vi hj
(32)
∂wi,j
N n
h
( v w h − α (v, h) )1 ∑∑i i,j ji,j−vi hj exp
.
w
ζ v
βi,j
h

Now, by substituting (30) into (32), ∂L+ /∂wi,j becomesh′where wi,jand Θ′ denote the estimates of wi,j and Θ updated−1.3

at the previous iteration. Since vi hj ∈ {0, 1}, we can write
the second term of (33) as

v

′

(

p(v, h|Θ )vi hj exp

h

= exp

(

′
wi,j − wi,j
w
βi,j

′
vi hj (wi,j − wi,j
)
w
βi,j

)∑∑
v

−1.4
−1.45

Contrastive Divergence
Auxiliary Function 1
Auxiliary Function 1 with gamma
Auxiliary Function 2
Auxiliary Function 2 with gamma

−1.55

)

−1.6
0

p(v, h|Θ′ )vi hj . (34)

100

0.5

h

1 ∑∑
(n)
p(h|v (n) , Θ′ )vi hj
N n
h
w
′
∑∑
log
.
+ βi,jwi,j = wi,jp(v, h|Θ′ )vi hjh(35)
Similarly to the CD approach, the numerator inside the logarithm function of (35) can be easily computed since it can be
written as (13). However, the denominator involves the sums
that run over all values of v and h. To avoid directly calculating these sums, one convenient way is to approximate theexpectation by samples from the RBM distribution, as with
the CD approach. As mentioned in 2.2, these samples can beobtained with Gibbs sampling using the conditional distributions (11) and (12). The update rule for the biases ai , bj canbe derived analogously.
It is interesting to compare the proposed update rule (35)with the update rule (8) of the CD algorithm. While eachstep of the CD algorithm moves the parameters in the direction given by the difference of the two expectations, the proposed update rule moves the parameters in the direction givenby the difference of the logarithms of the two expectations.
w
can be thought of as the learnFrom this perspective, βi,jw
ing rate. Since βi,j must satisfy the sum-to-unity constraint,
the larger the number of the hidden and visible variables bew
comes, the smaller βi,jbecomes in average, thus slowing theconvergence speed of the algorithm. We can thus expect to acw
w γcelerate the algorithm by replacing βi,jwith βi,j, setting γ
at a reasonably small value at the early stage of the algorithm
and moving it towards 1 as the iteration proceeds.
4. EXPERIMENTS
We conducted an experiment to compare the convergence
speeds of the present algorithm and the CD algorithm. In
order to evaluate the values of the likelihoods, We used asmall scale RBM with I = 10 visible units and J = 8 hidden

200
300
loop count

400

500

Contrastive Divergence
Auxiliary Function 1
Auxiliary Function 1 with gamma
Auxiliary Function 2
Auxiliary Function 2 with gamma

0.45This eventually allows us to solve ∂L+ /∂wi,j = 0:

v

−1.35

−1.5

reconstruct error rate

∑∑

x 10

−1.25

log−likelihood∂L+ (Θ)1 ∑∑(n)=p(h|v (n) , Θ′ )vi hj
(33)
∂wi,j
N n
h
(
)
′
∑∑
vi hj (wi,j − wi,j)p(v, h|Θ′ )vi hj exp
−
,
w
βi,j
v

4

−1.2

0.4
0.35
0.3
0.25
0.2
0.15
0.1
0

100

200
300
loop count

400

500

Fig. 2. Evolution of the mean log-likelihood (top) and the reconstruction error (bottom) during training of a small scale RBM using
the CD algorithm (black solid line) and the proposed algorithms (redand blue) with and without the γ scheduling (solid and dashed line),
respectively. “Auxiliary 1” refers to the the present algorithm and“Auxiliary 2” refers to the algorithm proposed in [16].

units. We also compared the present algorithm with anotheralgorithm we proposed in [16] (“auxiliary function 1” and“auxiliary function 2” hereafter).
The experimental conditions were as follows: 2000
randomly-generated binary data were used as the training
examples. All the algorithms were run for T = 500 iterations. The Gibbs sampling was run for 1 iteration at eachw
iterative step. β = {βia , βjb , βi,j|i = 1, . . . , I, j = 1, . . . , J}
was set uniformly, i.e.,w
βia = βjb = βi,j=1,
I + J + IJ

(36)

for all i and j. The learning rate ϵ of the CD algorithm andthe value of γ were scheduled in the following way: At step tof the algorithms, we updated ϵ and γ at(ϵt = ϵinit

ϵend
ϵinit

) Tt−1
−1
,

(37)

(
γt = γinit

γend
γinit

) Tt−1−1,(38)

respectively, where we set ϵinit = 1, ϵend = 0.1, γinit = 0.1
and γend = 1. In addition to the log-likelihood given by (6),
we used the reconstruction error rate E for the measure forcomparisons:1 ∑ ∑ (n)(n)(v − v̂i )2 ,NI n i i

E=

(39)

where
v̂ (n) = argmax p(v|ĥ

(n)

, Θ),

(40)

= argmax p(h|v (n) , Θ).

(41)

v

ĥ

(n)

h

Fig. 2 shows the evolution of the log-likelihoods (6) andthe reconstruction error rate (39) obtained with all the algorithms. With the MATLAB implementation, the execution
time per iteration for the CD algorithm and the proposed algorithms was almost the same. For the present algorithm, thelog-likelihood converged faster when updating γ than whenﬁxing it at 1, as expected. As Fig. 2 shows, the present algorithm increased the log-likelihood and decreased the reconstruction error rate of the RBM faster than the CD algorithm.
It is interesting to note that the present algorithm converged
to a greater value of the log-likelihood than the CD algorithm.
This implies that the present algorithm had a higher ability to
avoid getting trapped into local maxima.
5. CONCLUSIONSThis paper proposed deriving a new training algorithm for
RBMs based on the auxiliary function concept. While eachstep of the CD algorithm moves the parameters in the direction given by the difference of the two expectations, the proposed update rule moves the parameters in the direction givenby the difference of the logarithms of the two expectations.
Through an experiment on parameter training of an RBM, we
conﬁrmed that the present algorithm converged faster and toa better solution than the CD algorithm.
6. ACKNOWLEDGMENT
This work was supported by JSPS KAKENHI Grant Numbers
26730100.
7. REFERENCES[1] P. Smolensky, “Information processing in dynamical systems:
Foundations of harmony theory,” in Parallel Distributed Processing: Explorations in the Microstructure of Cognition, D.E.
Rumelhart and J.L McClelland, Eds. MIT Press, 1986.
[2] G. E. Hinton, S. Osindero, and Y. W. Teh, “A fast learningalgorithm for deep belief nets,” Neural Computation, vol. 18,
no. 7, 2006.

[3] G. E. Hinton, L. Deng, D. Yu, G. E. Dahl, A. Mohamed,
N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. N. Sainath,and B. Kingsbury, “Deep neural networks for acoustic modeling in speech recognition: The shared views of four researchgroups,” IEEE Signal Processing Magazine, vol. 29, no. 6, pp.
82–97, 2012.
[4] G. E. Hinton, “Training products of experts by minimizing
contrastive divergence,” Neural Computation, vol. 14, no. 8,
2002.
[5] J. De Leeuw and W. J. Heiser, “Convergence of correction
matrix algorithms for multidimensional scaling,” in Geometricrepresentations of relational data, J. C. Lingoes, E. E. Roskam,and I. Borg, Eds. Ann Arbor, MI: Mathesis Press, 1977.
[6] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum
likelihood from incomplete data via the EM algorithm,” J. ofRoyal Statistical Society Series B, vol. 39, 1977.
[7] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Adv. Neural Information and Processing
Systems (NIPS). 2001, pp. 556–562, MIT Press.
[8] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama, “Convergence-guaranteed multiplicativealgorithms for non-negative matrix factorization with betadivergence,” in Proc. 2010 IEEE International Workshop on
Machine Learning for Signal Processing (MLSP 2010), 2010,
pp. 283–288.
[9] H. Kameoka, Statistical Approach to Multipitch Analysis,
Ph.D. thesis, The University of Tokyo, 2007.
[10] H. Kameoka, N. Ono, and S. Sagayama, “Auxiliary function
approach to parameter estimation of constrained sinusoidal
model for monaural speech separation,” in Proc. 2008 IEEE
International Conference on Acoustics, Speech and Signal Processing (ICASSP2008), 2008, pp. 29–32.
[11] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A new sparse representation for acoustic signals,”in Proc. 2009 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP2009), 2009, pp. 3437–
3440.
[12] H. Kameoka, T. Nakatani, and T. Yoshioka, “Robust speech
dereverberation based on non-negativity and sparse natureof speech spectrograms,”
in Proc. 2009 IEEE International Conference on Acoustics, Speech and Signal Processing(ICASSP2009), 2009, pp. 45–48.
[13] N. Ono, “Stable and fast update rules for independent vector
analysis based on auxiliary function technique,” in Proc. 2011
IEEE Workshop on Applications of Signal Processing to Audioand Acoustics (WASPAA2011), 2011, pp. 261–264.
[14] N. Yasuraoka, H. Kameoka, T. Yoshioka, and H. G. Okuno, “Idivergence-based dereverberation method with auxiliary function approach,” in Proc. 2011 IEEE International Conferenceon Acoustics, Speech and Signal Processing (ICASSP2011),
2011, pp. 369–372.
[15] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Efﬁcientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” in Proc. 2012 IEEE International Conference on Acoustics, Speech and Signal Processing(ICASSP2012), 2012, pp. 261–264.
[16] N. Takamune and H. Kameoka, “Maximum reconstruction
probability training of restricted Boltzmann machines with
auxiliary function approach,” in Proc. 2014 IEEE International
Workshop on Machine Learning for Signal Processing (MLSP
2014), 2014, submitted.