arXiv:1412.7091v3 [] 14 Jul 2015

Efficient Exact Gradient Update for training Deep
Networks with Very Large Sparse Targets
[ Technical report ]
Pascal Vincent
Département d’Informatique et de Recherche Opérationnelle
Université de Montréal
Montréal, Québec, CANADA
and CIFAR
vincentp@iro.umontreal.ca
Alexandre de Brébisson
Département d’Informatique et de Recherche Opérationnelle
Université de Montréal
Montréal, Québec, CANADA
alexandre.de.brebisson@umontreal.ca
Xavier Bouthillier
Département d’Informatique et de Recherche Opérationnelle
Université de Montréal
Montréal, Québec, CANADA
xavier.bouthillier@iumontreal.ca

Abstract

An important class of problems involves training deep neural networks with sparse
prediction targets of very high dimension D. These occur naturally in e.g. neural
language models or the learning of word-embeddings, often posed as predicting
the probability of next words among a vocabulary of size D (e.g. 200 000). Computing the equally large, but typically non-sparse D-dimensional output vector
from a last hidden layer of reasonable dimension d (e.g. 500) incurs a prohibitive
O(Dd) computational cost for each example, as does updating the D × d output
weight matrix and computing the gradient needed for backpropagation to previous
layers. While efficient handling of large sparse network inputs is trivial, the case
of large sparse targets is not, and has thus so far been sidestepped with approximate alternatives such as hierarchical softmax or sampling-based approximations
during training. In this work we develop an original algorithmic approach which,
for a family of loss functions that includes squared error and spherical softmax,
can compute the exact loss, gradient update for the output weights, and gradient for backpropagation, all in O(d2 ) per example instead of O(Dd), remarkably
without ever computing the D-dimensional output. The proposed algorithm yields
D
a speedup of 4d
, i.e. two orders of magnitude for typical sizes, for that critical part
of the computations that often dominates the training time in this kind of network
architecture.

1

1

Introduction

Many modern applications of neural networks have to deal with data represented, or representable,
as very large sparse vectors. Such representations arise in natural language related tasks, where the
dimension D of that vector is typically (a multiple of) the size of the vocabulary, but also in the sparse
user-item matrices of collaborative-filtering applications. It is trivial to handle very large sparse
inputs to a neural network in a computationally efficient manner: the forward propagation and update
to the input weight matrix after backpropagation are correspondingly sparse. By contrast, training
with very large sparse prediction targets is problematic: even if the target is sparse, the computation
of the equally large network output and the corresponding gradient update to the huge output weight
matrix are not sparse and thus computationally prohibitive. This has been a practical problem ever
since Bengio et al. [1] first proposed using a neural network for learning a language model, in which
case the computed output vector represents the probability of the next word and is the size of the
considered vocabulary, which is becoming increasingly large in modern applications [2]. Several
approaches have been proposed to attempt to address this difficulty essentially by sidestepping it.
They fall in two categories:
• Sampling or selection based approximations consider and compute only a tiny fraction of
the output’s dimensions sampled at random or heuristically chosen. The reconstruction
sampling of Dauphin et al. [3], the efficient use of biased importance sampling in Jean
et al. [4], the use of Noise Contrastive Estimation [5] in Mnih and Kavukcuoglu [6] and
Mikolov et al. [7] all fall under this category. As does the more recent use of approximate
Maximum Inner Product Search based on Locality Sensitive Hashing techniques[8, 9] to
select a good candidate subset.
• Hierarchical softmax [10, 7] imposes a heuristically defined hierarchical tree structure for
the computation of the normalized probability of the target class.
Compared to the initial problem of considering all D output dimensions, both kinds of approaches
are crude approximations. In the present work, we will instead investigate a way to actually perform
the exact gradient update that corresponds to considering all D outputs, but do so implicitly, in a
computationally efficient manner, without actually computing the D outputs. This approach works
for a relatively restricted class of loss functions, the simplest of which is linear output with squared
error (a natural choice for sparse real-valued regression targets). The most common choice for multiclass classification, the softmax loss is not part of that class, but we may use an alternative spherical
softmax, which will also yield normalized class probabilities. For simplicity, our presentation will
focus on squared error and on an online setting, and we only later briefly mention its extension to
minibatches and to a more general class of loss functions.

2
2.1

The problem
Problem definition and setup

We are concerned with gradient-descent based training of a deep feed-forward neural network with
target vectors of very high dimension D (e.g. D = 200 000) but that are sparse, i.e. a comparatively
small number, at most K  D, of the elements of the target vector are non-zero. Such a Ksparse vector will typically be stored and represented compactly as 2K numbers corresponding
to pairs (index, value). A network to be trained with such targets will naturally have an equally
large output layer of dimension D. We can also optionally allow the input to the network to be a
similarly high dimensional sparse vector of dimension Din . Between the large sparse target, output,
and (optionally large sparse) input, we suppose the network’s intermediate hidden layers to be of
smaller, more typically manageable, dimension d  D (e.g. d = 500)1 .

Mathematical notation: Vectors are denoted using lower-case letters, e.g. h, and are considered
column-vectors; corresponding row vectors are denoted with a transpose, e.g. hT . Matrices are
denoted using upper-case letters, e.g. W , with W T the transpose of W . The ith column of W is
T
denoted Wi , and its ith row W:i (both viewed as a column vector). U −T = U −1 denotes the
transpose of the inverse of a square matrix. Id is the d × d identity matrix.

1
Our approach does not impose any restriction on the architecture nor size of the hidden layers, as long as
they are amenable to usual gradient backpropagation.

2

‣ Training deep neural networks with very large sparse targets is an important problem
‣ Arises e.g. in Neural Language Models [1] with large vocabulary size (e.g. D = 500 000 one-hot target).
‣ Efficient handling of large sparse inputs is trivial.
‣ But backprop training with large sparse targets is prohibitively expensive.
‣ Focus on output layer: maps last hidden representation h of reasonable dimension d (e.g. 500)
to very large output o of dimension D (e.g. 500 000) with a Dxd parameter matrix W:

Problem: expensive computation
Output o

o = Wh

O(Dd)

Backpropagation

(small d)

O(Kd)

Ex: d = 500

W(1)
(D d)
x

Input x

Altogether: O(3Dd )

We will now see how
only U and V respe
date
versions
Q=
5.3
Eﬃcient
with
Q
= W T Wof gra
hSupposing
in Equations
we ????
hav
The for
gradient
the
squ
updateof
b)
).will
2
d×d matrix
(we

h−y�
∂L
W is ∂W
= ∂�W∂W
ĥ = Qh
haswould
a comp
update
to
W
b
Solution:
representation
ofSeco
y,co
learning
rate. Again,
vector,
so that
a) Unew
=U
−comp
2η(
computational
complex
O(K).
So the all
overal
then
to
update
the
b) Vnew
=
V
+
2ηy
O(Kd
+ d2 ) = To
O((K
N
will
be sparse).
ove
of magnitude
cheape
(a
Proof:

If we define inter
computation of L caW

Note that we can d

Prohibitive!

O(d2) O(d2)

Boo

W(1) "W (1)- ! x !aT
O(Kd) cheap!

...

��

(large D, but K-sparse)

Ex: D = 500 000, K=5

�

Us
com

!

�

...

��

Target y (large D, but K-sparse)

�

Current workarounds are approximations:
Figure 1: ‣The
computational
problem
posed
very
large
sparse
Dealing
sparse inSampling
based approximations
compute
only by
a tiny
fraction
of the
output’stargets.
dimensions
sampled atwith
random.
put efficiently
is trivial,sampling
with both
the
and
backward
propagation
achieved in
Reconstruction
[2] and
theforward
use of Noise
Contrastive
Estimation
[3] in [4, 5]phases
fall undereasily
this category.
O(Kd). However
this is not the case with large sparse targets. They incur a prohibitive compu‣ Hierarchical softmax [6, 4] imposes a heuristically defined hierarchical tree structure for the computation of the
tational costnormalized
of O(Dd)
at the
output
layer as forward propagation, gradient backpropagation and
probability
of the
target class.
weight update each require accessing all D × d elements of the large output weight matrix.
References:

b) W ← W + 2ηyh

Vnew Un

...
W(2)(d d)
x

�

Forward propagation

O(Dd)

!h = W T !o

O(d2)

hidden 1

we suppose K << d << D

Prohibitivley expensive!

O(d2) O(d2)

hidden 2
(small d)

O(D)
!o = 2(o-y)

...

h

(small d)

O(D)

large D, not sparse

O(Dd)
W "W- ! !o hT

last hidden

Loss L = �o − y�2

a) W ← W − 2ηW

! Ne
" ac
SHORT FORMU
LINE CASE: !
Ne
" ac

Qn

a) First update
implicitly by updat

Proof:
[1] Bengio, Y., Ducharme, R., and Vincent, P. (2001). A neural probabilistic language model. NIPS 2000.
[2] Dauphin, Y., Glorot, X., and Bengio, Y. (2011). Large-scale learning of embeddings with reconstruction sampling. ICML 2011.

Network architecture: We consider a standard feed forward neural network architecture as depicted in Figure 1. An input vector x ∈ RDin is linearly transformed into a linear activation
a(1) = W (1)T x + b(1) through a Din × d input weight matrix W (1) (and an optional bias vector
b(1) ∈ Rd ). This is typically followed by a non-linear transformation s to yield the representation of
the first hidden layer h(1) = s(a(1) ). This first hidden layer representation is then similarly transformed through a number of subsequent non-linear layers (that can be of any usual kind amenable to
backpropagation) e.g. h(k) = s(a(k) ) with a(k) = W (k)T h(k−1) + b(k) until we obtain last hidden
layer representation h = h(m) . We then obtain the final D-dimensional network output as o = W h
where W is a D × d output weight matrix, which will be our main focus in this work. Finally, the
network’s D-dimensional output o is compared to the D-dimensional target vector y associated with
input x using squared error, yielding loss L = ko − yk2 .
Training procedure: This architecture is a typical (possibly deep) multi-layer feed forward neural
network architecture with a linear output layer and squared error loss. Its parameters (weight matrices and bias vectors) will be trained by gradient descent, using gradient backpropagation Rumelhart
et al. [11], LeCun [12, 13] to efficiently compute the gradients. The procedure is shown in Figure 1.
Given an example from the training set as an (input,target) pair (x, y), a pass of forward propagation
proceeds as outlined above, computing the hidden representation of each hidden layer in turn based
on the previous one, and finally the network’s predicted output o and associated loss L. A pass of
gradient backpropagation then works in the opposite direction, starting from ∇o = ∂L
∂o = 2(o − y)
∂L
and propagating back the gradients ∇h(k) = ∂h∂L
(k) and ∇a(k) = ∂a(k) upstream through the network. The corresponding gradient contributions on parameters (weights and biases), collected along
the way, are straightforward once we have the associated ∇a(k) . Specifically they are ∇b(k) = ∇a(k)
and ∇W (k) = h(k−1) (∇a(k) )T . Similarly for the input layer ∇W (1) = x(∇a(1) )T , and for the
output layer ∇W = (o − y)hT . Parameters are then updated through a gradient descent step
W (k) ← W (k) − η∇W (k) and b(k) ← b(k) − η∇b(k) , where η is a positive learning-rate. Similarly
for the output layer which will be our main focus here: W ← W − η∇W .
3

[3] G
[4] M

2.2

The easy part: input layer forward propagation and weight update

It is easy and straightforward to efficiently compute the forward propagation, and the backpropagation and weight update part for the input layer when we have a very large Din -dimensional but
K−sparse input vector x with appropriate sparse representation. Specifically we suppose that x is
represented as a pair of vectors u, v of length (at most) K, where u contains integer indexes and v
the associated real values of the elements of x such that xi = 0 if i ∈
/ u, and xuk = vk .

• Forward propagation through the input layer: The sparse representation of x as the
positions of K elements together with their value makes it cheap to compute W (1)T x.
Even though W (1) may be a huge full Din × d matrix, only K of its rows (those corresponding to the non-zero entries of x) need to be visited and summed to compute
W (1)T x. Precisely, with our (u, v) sparse representation of x this operation can be written
PK
(1)
(1)
asW (1)T x =
k=1 vk W:uk where each W:uk is a d-dimensional vector, making this an
O(Kd) operation rather than O(Dd).
• Gradient and update through input layer: Let us for now suppose that we were able to
get gradients (through backpropagation) up to the first hidden layer activations a(1) ∈ Rd
in the form of gradient vector ∇a(1) = ∂a∂L
(1) . The corresponding gradient-based update
(1)
(1)
to input layer weights W
is simply W
← W (1) − ηx(∇a(1) )T . This is a rank-one
(1)
update to W . Here again, we see that only the K rows of W (1) associated to the (at
most) K non-zero entries of x need to be modified. Precisely this operation can be written
(1)
(1)
as:W:uk ← W:uk − ηvk ∇a(1) ∀k ∈ {1, . . . , K} making this again a O(Kd) operation
rather than O(Dd).

2.3

The hard part: output layer propagation and weight update

Given some network input x we suppose we can compute without difficulty through forward propagation the associated last hidden layer representation h ∈ Rd . From then on:

• Computing the final output o = W h incurs a prohibitive computational cost of O(Dd)
since W is a full D×d matrix. Note that there is a-priori no reason for representation h to be
sparse (e.g. with a sigmoid non-linearity) but even if it was, this would not fundamentally
change the problem since it is D that is extremely large, and we supposed d reasonably
sized already. Computing the residual (o − t) and associated squared error loss ko − tk2
incurs an additional O(D) cost.
T
• The gradient on h that we need to backpropagate to lower layers is ∇h = ∂L
∂h = 2W (o −
y) which is another O(Dd) matrix-vector product.
• Finally, when performing the corresponding output weight update W ← W − η(o − y)hT
we see that it is a rank-one update that updates all D ×d elements of W , which again incurs
a prohibitive O(Dd) computational cost.

For very large D, all these three O(Dd) operations are prohibitive, and the fact that y is sparse, seen
from this perspective, doesn’t help, since neither o nor o − y will be sparse.

3

A computationally efficient algorithm for performing the exact online
gradient update

Previously proposed workarounds are approximate or use stochastic sampling. We propose a different approach that results in the exact same, yet efficient gradient update, remarkably without ever
having to compute large output o.
3.1

Computing the squared error loss L and the gradient with respect to h efficiently

Suppose that, we have, for a network input example x, computed the last hidden representation
h ∈ Rd through forward propagation. The network’s D dimensional output o = W h is then in
principle compared to the high dimensional target y ∈ RD . The corresponding squared error loss
2
is L = kW h − yk . As we saw in Section 2.3, computing it in the direct naive way would have
a prohibitive computational complexity of O(Dd + D) = O(Dd) because computing output W h
4

with a full D × d matrix W and a typically non-sparse h is O(Dd). Similarly, to backpropagate
the gradient through the network, we need to compute the gradient of loss L with respect to last
∂kW h−yk2
hidden layer representation h. This is ∇h = ∂L
= 2W T (W h − y). So again, if
∂h =
∂h
we were to compute it directly in this manner, the computational complexity would be a prohibitive
O(Dd). Provided we have maintained an up-to-date matrix Q = W T W , which is of reasonable
size d × d and can be cheaply maintained as we will see in Section 3.3, we can rewrite these two
operations so as to perform them in O(d2 ):
Loss computation:

Gradient on h:

O(Dd)

z}|{
k W h −yk2

L =

∇h =

∂L
∂h

=

T

(W h − y) (W h − y)

=

=

hT W T W h − y T W h − hT W T y + y T y
hT Qh − 2hT (W T y) + y T y

=
=

hT ( Qh −2 W T y ) + y T y
|{z}
| {z }
|{z}

=

O(d2 )

O(Kd)

=

=

(1)

∂kW h − yk2
∂h
2W T (W h − y)

2 WTWh − WTy
2( Qh − W T y )
|{z} | {z }
O(d2 )

O(K)



O(Kd)

The terms in O(Kd) and O(K) are due to leveraging the K-sparse representation of target vector
y. With K  D and d  D, we get altogether a computational cost of O(d2 ) which can be several
orders of magnitude cheaper than the prohibitive O(Dd) of the direct approach.
3.2

Efficient gradient update of W

The gradient of the squared error loss with respect to output layer weight matrix W is
∂kW h−yk
∂W

2

∂L
∂W

=

= 2(W h − y)h . And the corresponding gradient descent update to W would be
Wnew ← W − 2η(W h − y)hT , where η is a positive learning rate. Again, computed in this manner,
this induces a prohibitive O(Dd) computational complexity, both to compute output and residual
W h − y, and then to update all the Dd elements of W (since generally neither W h − y nor h will
be sparse). All D × d elements of W must be accessed during this update. On the surface this seems
hopeless. But we will now see how we can achieve the exact same update on W in O(d2 ). The trick
is to represent W implicitly as the factorization |{z}
W = |{z}
V |{z}
U and update U and V instead:
T

D×d

a) Unew

b) Vnew

D×d d×d

= U − 2η(U h)hT
=

V +

−T
2ηy(Unew
h)T ù

(3)

(4)

This results in implicitly updating W as we did explicitly in the naive approach as we now prove:
Vnew Unew

=

−T
(V + 2ηy(Unew
h)T ) Unew

−T
= V Unew + 2ηy(Unew
h)T Unew
−1
= V Unew + 2ηyhT Unew
Unew
−1
= V (U − 2η(U h)hT ) + 2ηyhT (Unew
Unew )

= V U − 2ηV U hhT + 2ηyhT
= V U − 2η(V U h − y)hT
= W − 2η(W h − y)T hT
= Wnew

We see that the update of U in Eq. 3 is a simple O(d2 ) operation. Following this simple rank-one
update to U , we can use the Sherman-Morrison formula to derive the corresponding rank-one update
to U −T which will also be O(d2 ):
5

(2)

−T
Unew

=

U −T +

2η

2 (U

1 − 2η khk

−T

h)hT

(5)

−T
It is then easy to compute the Unew
h, an O(d2 ) operation needed in Eq. 4. The ensuing rank-one
update of V in Eq 4, thanks to the K-sparsity of y is only O(Kd): only theK rows V associated to
non-zero elements in y are accessed and updated, sited of all D rows of W we had to modify in the
naive update!

Note that with the factored representation of W as V U , we only have W implicitly, so the W T y
terms that entered in the computation of L and ∇h in the previous paragraph need to be adapted
slightly as ŷ = W T y = U T (V T y), which becomes O(d2 + Kd) rather than O(Kd) in computational complexity. But this doesn’t change the overall O(d2 ) complexity of these computations.
3.3

Bookkeeping: keeping an up-to-date Q and U −T

We have already seen, in Eq. 5, how we can cheaply maintain an up-to-date U −T following our
update of U . Similarly, following our updates to U and V , we need to keep an up-to-date Q =
W T W which is needed to efficiently compute the loss L (Eq. 1) and gradient ∇h (Eq. 2). We have
shown that updates to U and V in equations 3 and 4 are equivalent to implicitly updating W as
Wnew ← W − 2η(W h − y)hT , and this translates into the following update to Q = W T W :
ẑ
Qnew

=

Qh − U T (V T y)


= Q − 2η hẑ T + ẑhT + (4η 2 L)hhT

(6)

The proof is straightforward but due to space constraints we put it in supplementary material. One
can see that this last bookkeeping operation also has a O(d2 ) computational complexity.

3.4

Putting it all together: detailed algorithm and expected benefits

We have seen that we can efficiently compute cost L, gradient with respect to h (to be later backpropagated further) as well as updating U and V and performing the bookkeeping for U −T and
Q. Algorithm 1 describes the detailed algorithmic steps that we put together from the equations
derived above. Having K  d  D we see that the proposed algorithm requires O(d2 ) operations,
whereas the standard approach required O(Dd) operations. If we take K ≈ d , we may state more
precisely that the proposed algorithm, for computing the loss and the gradient updates will require
roughly 12d2 operations whereas the standard approach required roughly 3Dd operations. So overD
all the proposed algorithm change corresponds to a computational speedup by a factor of 4d
. For
D = 200 000 and d = 500 the expected speedup is thus 100. Note that the advantage is not only
in computational complexity, but also in memory access. For each example, the standard approach
needs to access and change all D × d elements of matrix W , whereas the proposed approach only
accesses the much smaller number K × d elements of V as well as the three d × d matrices U , U −T ,
and Q. So overall we have a substantially faster algorithm, which, while doing so implicitly, will
nevertheless perform the exact same gradient update as the standard approach. We want to emphasize here that our approach is completely different from simply chaining 2 linear layers U and V
and performing ordinary gradient descent updates on these: this would result in the same prohibitive
computational complexity as the standard approach, and such ordinary separate gradient updates to
U and V would not be equivalent to the ordinary gradient update to W = V U .
3.5

Controlling numerical stability and extension to the minibatch case

The update of U in Equation 3 may over time lead U to become ill-conditioned. To prevent this,
we regularly (every 100 updates) monitor its conditioning number2 . If either the smallest or largest
singular value moves outside an acceptable range, we bring it back to 1 by doing an appropriate
2

Largest and smallest singular value can be computed with an SVD or using the power iteration method.

6

Algorithm 1 Efficient computation of cost L, gradient h, and update to parameters U and V
Step Operation
Computational Number of
#
complexity
multiply-adds
1:
2:
3:
4:
5:
6:
7:
8:
9:

ĥ = Qh
ŷ = U T (V T y)
ẑ = ĥ − ŷ
∇h = 2ẑ
L = hT ĥ − 2hT ŷ + y T y
Unew = U − 2η(U h)hT
−T
Unew
=
2η
−T
h)hT
U −T + 1−2ηkhk
2 (U
−T
Vnew = V + 2ηy(Unew
h)T
Qnew =

Q − 2η hẑ T + ẑhT +
(4η 2 L)hhT
Altogether:

O(d2 )
O(Kd + d2 )
O(d)
O(d)
O(2d + K)
O(d2 )
O(d2 )

d2
Kd + d2
d
d
2d + K + 1
2d2 + d
2
2d + 2d + 3

O(d2 + Kd)
O(d2 )

d2 + K + Kd
4 + 2d + 3d2

O(d2 )
provided
K<dD

≈ 12d2
elementary
operations

rank-1 update to V (which costs Dd operations, but is only done rarely). Our algorithm can also
be straightforwardly extended to the minibatch case (the derivations are given in the supplementary material section) and yields the same theoretical speedup factor with respect to the standard
naive approach. But one needs to be careful in order to keep the computation of U −T h reasonably
efficient: depending on the size of the minibatch m, it may be more efficient to solve the corresponding linear equation for each minibatch from scratch rather than updating U −T with the Woodbury
equation (which generalizes the Sheman-Morrison formula for m > 1).
3.6

Generalization to a broader class of loss functions

The approach that we detailed for linear output and squared error can easily be extended to slightly
more exotic loss functions: basically any
P loss function that can be expressed using only the oc
associated to non-zero yc and kok2 = j o2j the squared norm of the whole output vector, which
we can compute cheaply. This family of loss functions does not include the standard log of softmax,
2
c +)
but includes the so-called spherical softmax: log P(o(o
2 (where c is the correct class label). It
+)
j
j
remains to be seen in practice how this approach performs computationally, and whether we lose
something due to using this more limited family of loss functions.

4

Experimental validation

We implemented both a CPU version using blas and a parallel GPU (Cuda) version using cublas
of the proposed algorithm3 . We evaluated the GPU and CPU implementations by training word
embeddings with simple neural language models, in which a probability map of the next word
given its preceding n-gram is learned by a neural network. We used a Nvidia Titan Black GPU
and a i7-4820K @ 3.70GHz CPU and ran experiments on the one billion word dataset[14], which
is composed of 0.8 billions words belonging to a vocabulary of 0.8 millions words. We evaluated
the resulting word embeddings with the recently introduced Simlex-999 score [15], which measures
the similarity between words. We also compared our approach to unfactorised versions and to a
two-layer hierarchical softmax. Figure 2 and 3 (left) illustrate the practical speedup of our approach
for the output layer only. Figure 3(right) shows that the LST (Large Sparse Target) models are much
faster to train than the softmax models and converge to only slightly lower Simlex-999 scores. Table
1 summarizes the speedups for the different output layers we tried, both on CPU and GPU. We also
emprically verified that our proposed factored algorithm learns the exact same model weights (V U )
3

Open source code will be released upon official publication of this research.

7

Table 1: Speedups with respect to the baseline naive model on CPU, for a minibatch of 128 and the
whole vocabulary of D = 793471 words. This is a two hidden layer model with 300 neurons on all
its layers (so d = 300).
Model
output layer only speedup whole model speedup
cpu unfactorised (naive)
1
1
gpu unfactorised (naive)
6.8
4.7
gpu hierarchical softmax
125.2
178.1
cpu factorised
763.3
501
gpu factorised
3257.3
1852.3

101

un-factorised CPU
un-factorised GPU
factorised GPU
factorised CPU
h_softmax GPU

0.008

Timing (sec) of a minibatch of size 128

Timing (sec) of a minibatch of size 128

0.010

0.006
0.004
0.002
0.000
0

2000

4000
6000
Size of the vocabulary D

8000

100

10-1

10-2

10-3 1
10

10000

un-factorised CPU
un-factorised GPU
factorised GPU
factorised CPU
h_softmax GPU

102

104
103
Size of the vocabulary D

105

106

Figure 2: Timing of different algorithms. Time taken by forward and backward propagations in the
output layer, including weight update, on a minibatch of size 128 for different sizes of vocabulary
D on both CPU and GPU. The input size d is fixed to 300. The Timing of a 2 layer hierarchical
softmax efficient GPU implementation (h_softmax) is also provided for comparison. Right plot is
in log-log scale. As expected, the timings of factorized versions are independent of the size of the
vocabulary.
as the corresponding naive unfactored algorithm’s W , as it theoretically should (up to negligible
numerical precision differences), and followed the exact same learning curves (as a function of
number of iterations, not time!).

5

Conclusion and future work

We introduced a new algorithmic approach to efficiently compute the exact gradient updates for
training deep networks with very large sparse targets. Remarkably the complexity of the algorithm
is independent of the target size, which allows tackling very large problems. Our CPU and GPU
implementation yield similar speedups to the theoretical one and can thus be used in practical applications, which could be explored in further work. In particular, neural language models seem
good candidates. But it remains unclear how using a loss function other than log-softmax may affect
the quality of the resulting word embeddingsm and further research should be carried out in this
direction. Extensions of the approach to other possible losses than the simple squared error should
also be empirically investigated in this light, in particular log-spherical-softmax.

Acknowledgements
We would like to thank the developers of Theano [16, 17] and Blocks [18].
This research is supported by NSERC and Ubisoft.
8

1200

0.20
0.15
SimLex-999

1000
Speedup

0.25

cpu_unfact / cpu_fact, experimental
gpu_unfact / gpu_fact, experimental
unfact / fact, theoretical
cpu_unfact / gpu_fact, experimental
cpu_unfact / gpu_unfact, experimental

1400

800
600

0.10
0.05

400

0.00

200

0.05

0
0

100

200 300 400 500 600 700
Size of the vocabulary D (in thousands)

0.10

800

LST CPU
LST GPU
Softmax CPU
Softmax GPU
H-Softmax GPU
10-1

100

101
Training time (hours)

102

103

Figure 3: Left: Practical and theoretical speedups for different sizes of vocabulary D and fixed input
size d=300. The practical unfact / fact speedup is similar to the theoretical one. Right: Evolution
of the Simlex-999 score obtained with different models as a function of training time (CPU softmax
times were extrapolated from fewer iterations). Softmax models are zero hidden-layer models, while
our large sparse target (LST) models have two hidden layers. These were the best architectures
retained in both cases (surprisingly the softmax models with hidden layers performed no better on
this task). The extra non-linear layers in LST may help compensate for the lack of a softmax. LST
models converge to slightly lower scores at similar speed as the hierarchical softmax model but
significantly faster than softmax models.

References
[1] Yoshua Bengio, Réjean Ducharme, and Pascal Vincent. A neural probabilistic language model. In
NIPS’00, pages 932–938. MIT Press, 2001.
[2] R. Collobert, J. Weston, L. Bottou, M. Karlen, K. Kavukcuoglu, and P. Kuksa. Natural language processing (almost) from scratch. Journal of Machine Learning Research, 12:2493–2537, 2011.
[3] Y. Dauphin, X. Glorot, and Y. Bengio. Large-scale learning of embeddings with reconstruction sampling.
In Proceedings of the 28th International Conference on Machine learning, ICML ’11, 2011.
[4] Sébastien Jean, Kyunghyun Cho, Roland Memisevic, and Yoshua Bengio. On using very large target
vocabulary for neural machine translation. In ACL-IJCNLP’2015, 2015. arXiv:1412.2007.
[5] M. Gutmann and A. Hyvarinen. Noise-contrastive estimation: A new estimation principle for unnormalized statistical models. In Proceedings of The Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS’10), 2010.
[6] Andriy Mnih and Koray Kavukcuoglu. Learning word embeddings efficiently with noise-contrastive
estimation. In C.J.C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors,
Advances in Neural Information Processing Systems 26, pages 2265–2273. Curran Associates, Inc., 2013.
[7] T. Mikolov, I. Sutskever, K. Chen, G.S. Corrado, and J. Dean. Distributed representations of words and
phrases and their compositionality. In NIPS’2013, pages 3111–3119. 2013.
[8] Anshumali Shrivastava and Ping Li. Asymmetric LSH (ALSH) for sublinear time maximum inner product
search (MIPS). In Z. Ghahramani, M. Welling, C. Cortes, N.D. Lawrence, and K.Q. Weinberger, editors,
Advances in Neural Information Processing Systems 27, pages 2321–2329. Curran Associates, Inc., 2014.
[9] Sudheendra Vijayanarasimhan, Jonathon Shlens, Rajat Monga, and Jay Yagnik. Deep networks with large
output spaces. arxiv:1412.7479, 2014.
[10] Frederic Morin and Yoshua Bengio. Hierarchical probabilistic neural network language model. In
Robert G. Cowell and Zoubin Ghahramani, editors, Proceedings of the Tenth International Workshop
on Artificial Intelligence and Statistics, pages 246–252. Society for Artificial Intelligence and Statistics,
2005.

9

[11] D.E. Rumelhart, G.E. Hinton, and R.J. Williams. Learning representations by back-propagating errors.
Nature, 323:533–536, 1986.
[12] Yann LeCun. Une procédure d’apprentissage pour Réseau à seuil assymétrique. In Cognitiva 85: A
la Frontière de l’Intelligence Artificielle, des Sciences de la Connaissance et des Neurosciences, pages
599–604, Paris 1985, 1985. CESTA, Paris.
[13] Yann LeCun. Learning processes in an asymmetric threshold network. In E. Bienenstock, F. FogelmanSoulié, and G. Weisbuch, editors, Disordered Systems and Biological Organization, pages 233–240.
Springer-Verlag, Berlin, Les Houches 1985, 1986.
[14] Ciprian Chelba, Tomas Mikolov, Mike Schuster, Qi Ge, Thorsten Brants, Phillipp Koehn, and Tony
Robinson. One billion word benchmark for measuring progress in statistical language modeling. In
INTERSPEECH 2014, 15th Annual Conference of the International Speech Communication Association,
Singapore, September 14-18, 2014, pages 2635–2639, 2014.
[15] Felix Hill, Roi Reichart, and Anna Korhonen. Simlex-999: Evaluating semantic models with (genuine)
similarity estimation. CoRR, abs/1408.3456, 2014.
[16] James Bergstra, Olivier Breuleux, Frédéric Bastien, Pascal Lamblin, Razvan Pascanu, Guillaume Desjardins, Joseph Turian, David Warde-Farley, and Yoshua Bengio. Theano: a CPU and GPU math expression compiler. In Proceedings of the Python for Scientific Computing Conference (SciPy), 2010. Oral
Presentation.
[17] Frédéric Bastien, Pascal Lamblin, Razvan Pascanu, James Bergstra, Ian J. Goodfellow, Arnaud Bergeron,
Nicolas Bouchard, and Yoshua Bengio. Theano: new features and speed improvements. Deep Learning
and Unsupervised Feature Learning NIPS 2012 Workshop, 2012.
[18] B. van Merriënboer, D. Bahdanau, V. Dumoulin, D. Serdyuk, D. Warde-Farley, J. Chorowski, and Y. Bengio. Blocks and Fuel: Frameworks for deep learning. ArXiv e-prints, June 2015.

10

Appendix
A

Minibatch version of the algorithm

The algorithm we derived for online gradient is relatively straightforward to extend to the case of
minibatches containing m examples, and will still yield the same theoretical speedup factor with
respect to the standard naive approach. One may want to be careful in order to keep the computation
of U −T h (or ore precisely U −T H in the minibatch case) reasonably efficient. In the minibatch
version presented below, we update U −T based on the Woodbury equation (which generalizes the
Sheman-Morrison formula for m > 1 and involves inverting an m × m matrix). But depending
on the size of the minibatch m, it may become more efficient to solve the corresponding linear
equations for each minibatch from scratch every time, rather than inverting that m × m matrix. In
which case we won’t need to maintain an U −T at all.
Algorithm 2 Minibatch version of the update algorithm
Initialization
• we can initialize D × d matrix V randomly as we would have initialized W so that we
initially have V = W .
Alternatively we can initialize V to 0 (there won’t be symmetry breaking issues with having
W initially be 0 provided the other layers are initialized randomly, since varying inputs and
targets will naturally break symmetry for the output layer)
• initialize Q ← V T V (or more cheaply initialize Q ← 0 if we have initialized V to 0).
• we initialize U to the identity: U ← Id so that, trivially, we initially have V U = W .
• initialize U −T ← Id

Update
We suppose we receive m target vectors in the m columns of sparse matrix Y , and corresponding
m hidden representations in the m columns of matrix H.
Step Operation
Computation
Computational complexity
#
complexity
with the multiplicative factor
left in.
1:
2:

Ĥ = QH
Ŷ = U T (V T Y )

3:
4:
5:

Ẑ = Ĥ − Ŷ
∇H = 2Ẑ
M = H T Ẑ − Ŷ T H + Y T Y or
alternatively
M = H T Ĥ − (Ŷ T H + H T Ŷ ) + Y T Y
L = Tr(M )
Unew = U − 2η(U H)H T
−T
Unew
= U−T −

1
−T
(U H) (H T H − 2η
Im )−1 H T

6:
7:
8:
9:
10:

−T
Vnew = V + 2ηY (Unew
H)T



Qnew = Q − 2η H Ẑ T + ẐH T +
4η 2 (HM )H T

11

O(md2 )
O(mKd +
md2 )
O(md)
O(md)
O(m2 d +
m2 K)

O(md2 )
O(mKd + md2 )

O(m)
O(md2 )
O(m2 d +
m3 + md2 )

O(m)
O(2md2 )
O(2m2 d + m3 + 2md2 )

O(md)
O(md)
O(2m2 d + m2 K)

O(md2 +
mKd)

O(md2 + mKd)

O(md2 +
dm2 )

O(3md2 + m2 d)

B

Detailed proof for computation of update of Q

Update to Q corresponds to Wnew ← W − 2η(W H − Y )H T
We will use the following precomputed quantities: Q = W T W , Ĥ = QH and Ŷ = W T Y =
U T (V T Y ) and Ẑ = Ĥ − Ŷ .
Qnew

=
=
=
=
=
=

=
=

T
Wnew
Wnew

W − 2η(W H − Y )H T

T

W − 2η(W H − Y )H T



W T W − 2ηH(W H − Y )T W − 2ηW T (W H − Y )H T

+4η 2 H(W H − Y )T (W H − Y )H T


Q − 2η HH T W T W − HY T W − 2η W T W HH T − W T Y H T
+4η 2 H(H T W T W H − H T W T Y − Y T W H + Y T Y )H T


Q − 2η HH T Q − H(W T Y )T − 2η QHH T − (W T Y )H T
+4η 2 H(H T QH − H T (W T Y ) − (W T Y )T H + Y T Y )H T


Q − 2η H Ĥ T − H Ŷ T + ĤH T − Ŷ H T

+4η 2 H(H T Ĥ − H T Ŷ − Ŷ T H + Y T Y )H T


Q − 2η H(Ĥ − Ŷ )T + (Ĥ − Ŷ )H T + 4η 2 H(H T (Ĥ − Ŷ ) − Ŷ T H + Y T Y )H T




Q − 2η H Ẑ T + ẐH T + 4η 2 H H T Ẑ − Ŷ T H + Y T Y H T
{z
}
|
M

This is what is listed as step 10 of the above minibatch algorithm.
In the online case, this becomes:
Qnew



Q − 2η hẑ T + ẑhT + 4η 2 hT ẑ − ŷ T h + y T y hhT



= Q − 2η hẑ T + ẑhT + 4η 2 hT ĥ − hT ŷ − ŷ T h + y T y hhT



= Q − 2η hẑ T + ẑhT + 4η 2 hT ĥ − 2hT ŷ + y T y hhT

= Q − 2η hẑ T + ẑhT + (4η 2 L)hhT

=

which is the update listed as step 9 in the online algorithm (Algorithm 1 on page 7).

12

C

Details regarding controlling numerical stability

The update of U (step 6 of the online algorithm, step 7 in the minibatch version) may over time lead
to U becoming ill-conditioned. Simultaneously, as we update U and U −T (using Sherman-Morrison
or Woodbury) our updated U −T may numerically start to diverge from the true U −T due to numerical
precision. It is thus important to prevent both of these form happening, i.e. make sure U stays well
conditioned, to ensure the numerical stability of the algorithm. We present here progressively refined
strategies for achieving this.

Restoring the system in a pristine stable state
One simple way to ensure numerical stability is to once in a while restore the system in its pristine
state where V = W and U = Id = U −T . This is easily achieved as follows:

V
U
U −T

←
←
←

VU
Id
Id .

This operation doesn’t affects the product V U , so the implicit matrix W remains unchanged, nor
does it affect Q = W T W . And it does restore U to a perfectly well conditioned identity matrix. But
computing V U is an extremely costly O(Dd2 ) operation, so if possible we want to avoid it (except
maybe once at the very end of training, if we want to compute the actual W ). In the next paragraphs
we develop a more efficient strategy.

Stabilizing only problematic singular values
U becoming ill-conditioned is due to its singular values over time becoming too large and/or too
small. Let use define σ1 , . . . , σd as the singular values of U ordered in decreasing order. The
conditioning number of U is defined as σσd1 and it can become overly large when σ1 becomes too
large and/or when σd becomes too small. Restoring the system in its pristine state, as shown in the
previous paragraph, in effect brings back all singular values of U back to 1 (since it brings back
U to being the identity). It is instead possible, and computationally far less costly, to correct when
needed only for the singular values of U that fall outside a safe range. Most often we will only
need to occasionally correct for one singular value (usually the smallest, and only when it becomes
too small). Once we have determined the offending singular value and its corresponding singular
vectors, correcting for that singular value, i.e. effectively bringing it back to 1, will be a O(Dd)
operation. The point is to apply corrective steps only on the problematic singular values and only
when needed, rather than blindly, needlessly and inefficiently correcting for all of them through the
basic O(Dd2 ) full restoration explained in the previous paragraph.
13

Here is the detailed algorithm that achieves this:
Algorithm 3 Numerical stabilization procedure for problematic singular values
• The chosen safe range for singular values is [σlow , σhigh ] (ex: [0.001, 100] )
• The procedures given below act on output layer parameters U , U −T and V .
• For concision, we do not enlist these parameters explicitly in their parameter list.
• Procedure SINGULAR - STABILIZE gets called after every ncheck gradient updates (ex:
ncheck = 100).
procedure SINGULAR - STABILIZE( )
Ū, σ, V̄ = SVD(U ) . Computes singular value decomposition of U as U = Ū diag(σ) V̄T
for all k ∈ {1, . . . , d} do
if σk < σlow OR σk > σhigh then
FIX - SINGULAR - VALUE(σk , Ūk , 1)
end if
end for
end procedure
The following procedure will change singular value σ of U associated to singular vector u to
become target singular value σ ∗ (typically 1). It doesn’t change U ’s singular vectors, only that
one singular value. It also changes V symetrically (with a rank-one update) in such a way that
W = V U remains unchanged.
procedure FIX - SINGULAR - VALUE(σ, u, σ ∗ )
∗
α = σ σ−σ
α
β = − 1+α
U ← U + αu(U T u)T
V ← V + β(V u)uT
U −T ← U −T + βu(U −1 u)T
. Where U −1 is obtained as the transpose of U −T . But
we may instead of this prefer to recompute U −T from scratch by inverting U to ensure it doesn’t
stray too much due to numerical imprecisions.
end procedure
The proof that the FIX - SINGULAR - VALUE procedure achieves what it is supposed to is relatively
straightforward, and left to the reader.
Avoiding the cost of a full singular-value decomposition
Computing the SVD of d × d matrix U as required above, costs roughly 25d3 elementary operations
(use the so-called R - SVD algorithm). But since the offending singular values will typically be only
the smallest or the largest, it is wasteful to compute all d singular values every time. A possibly
cheaper alternative is to use the power iteration method with U to find its largest singular value and
associated singular vector, and similarly with U −1 to obtain the smallest singular value of U (which
corresponds to the inverse of the largest singular value of U −1 ). Each iteration of the power iteration
method requires only O(d2 ) operations, and a few iterations may suffice. In our experiments we
fixed it to 100 power iterations. Also it is probably not critical if the power iteration method is not
run fully to convergence, as correcting along an approximate offending singular vector direction can
be sufficient for the purpose of ensuring numerical stability.
With this refinement, we loop over finding the smallest singular value with the power iteration
method, correcting for it to be 1 by calling FIX - SINGULAR - VALUE if it is too small, and we repeat
this until we find the now smallest singular value to be inside the acceptable range. Similarly for the
largest singular values.
Note that while in principle we may not need to ever invert U from scratch (as we provided update
formulas of U −T with every change we make to U ), it nevertheless proved to be necessary to do so
regularly to ensure U −T doesn’t stray too much from the correct value due to numerical imprecisions. Inverting U using Gaussian-elimination costs roughly d3 operations, so it is very reasonable
14

and won’t affect the computational complexity if we do it no more often than every d training examples (which will typically correspond to less than 10 minibatches of size 128). In practice, we
recompute U −T from scratch every time before we run this check for singular value stabilization.

15

