{'fraction_non_alphanumeric': 0.03205712945346504, 'fraction_numerical': 0.010125361620057859, 'mean_word_length': 4.321336939185909, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parallel training of DNNs with Natural Gradient and Parameter Averaging

We describe the neural-network training framework used in the Kaldi speech recognition toolkit, which is geared towards training DNNs with large amounts of training data using multiple GPU-equipped or multi-core machines. In order to be as hardware-agnostic as possible, we needed a way to use multiple machines without generating excessive network traffic. Our method is to average the neural network parameters periodically (typically every minute or two), and redistribute the averaged parameters to the machines for further training. Each machine sees different data. By itself, this method does not work very well. However, we have another method, an approximate and efficient implementation of Natural Gradient for Stochastic Gradient Descent (NG-SGD), which seems to allow our periodic-averaging method to work well, as well as substantially improving the convergence of SGD on a single machine.


Introduction
In this paper we introduce two ideas: a data-parallel training method based on periodically averaging the parameters of separate SGD runs; and an efficient and practical implementation of Natural Gradient for Stochastic Gradient Descent (NG-SGD) for deep neural networks.
In Section 2 we describe our problem setting, which is Deep Neural Networks (DNNs) applied to speech recognition-although our ideas are more general than this. In Section 3 we introduce the parallel training method. In Section 4 we describe the general ideas behind our natural gradient method, although most of the technical details have been relegated to appendices. In this paper we don't give any proofs, but we do discuss in Section 5 what we think we can and can't be proven about our methods. Section 6 has experiments on the convergence of SGD with and without natural gradient and parallelism. We conclude in Section 7.
There are two versions of our NG-SGD method: a "simple" version and an "online" one. Technical details for these are in Appendices A and B respectively. Appendix C has background information on our DNN implementation.

Problem setting
When training DNNs for speech recognition, the immediate problem is that of classifying vectors x ∈ R D as corresponding to discrete labels y ∈ Y. The dimension D is typically several hundred, with x being derived from short-time spectral properties of the acoustic signal; and Y corresponds to clustered HMM states of context-dependent phones, |Y| ≃ 5000 being typical. Each (x, y) pair corresponds to a single frame of speech data; frames are typically extracted at the rate of 100 per second, with a duration of 25 milliseconds, and x contains spectral information from several adjacent frames spliced together [1]. We are ultimately not just interested in the top y on each frame, but in the log-probabilities log p(y|x) for all y, since we will use them as costs in a Viterbi search for a path corresponding to the most likely word sequence. The objective function for training is the sum, over all frames of training data, of the log-probability of y given x: i log p(y i |x i ). The supervision labels y are derived from a Viterbi alignment of a Hidden Markov Model (HMM) derived from the reference word sequence of each training utterance.

Parameter-averaging overview
The parameter-averaging aspect of our training is quite simple. We have N machines (e.g. N = 4) each doing SGD separately with different randomized subsets of the training data, and we allow their parameters to gradually diverge. After each machine has processed a fixed number of samples K (typically K = 400 000), we average the parameters across all the jobs and re-distribute the result to the machines. (In practice we do this by spawning new jobs in Gri-dEngine or in whatever job management system we are using). This is repeated until we have processed all the data for a specified number of epochs, e.g. 10.
We define an outer iteration of training as the time it takes for each job to process K training examples. The number of outer iterations per epoch depends on K and the quantity of training data.
We find it useful to define the effective learning rate of the parallel SGD procedure as the learning rate η t being used by the individual jobs, divided by the number of jobs N . As we increase the number of jobs N , in order to get a linear speed up we need to increase the learning rate proportional to N so that the effective learning rate stays the same. The idea is that when we do the parameter averaging, the parameter update from any individual SGD job gets diluted N -fold.

CPU versus GPU-based SGD
Each machine in our parallel computation implements SGD. We have two versions of this, one for GPU and one for CPU. The GPU-based computation is standard minibatch-based SGD, typically with 512 examples per minibatch.
In the CPU-based computation, each job uses a specified number of threads (typically 16) in order to take advantage of multi-core processors. Our parallelization method uses shared parameters but no parameter locking; this is referred to in [2] as asynchronous. The same idea is used in Hogwild! [3]. We also implemented a lock-free parallel SGD. (We avoid the term Asynchronous SGD, which seems to lack a commonly accepted definition.) In order to prevent divergence, each thread processes relatively small minibatchestypically, of size 128.
We should mention at this point that in our formulation, we sum the gradients over the elements of the minibatch, rather than averaging: this ensures that we make the same amount of progress per sample, regardless of minibatch size, and so gives more consistent results when changing the minibatch size. The need to limit the minibatch size in the multithreaded case can be understood as follows: think of the effective minibatch size as being the minibatch size times the number of threads. The product of the learning rate η with the effective minibatch size is relevant for stability of the SGD update: if it becomes too large, there is increased danger of divergence.
We typically use the GPU-based method, because in our experience a GPU can process data many times faster than even a 16-threaded process running on CPUs. Aside from speed, the two methods give very similar results.

Data organization and sequential data access
On spinning hard disks, sequential data access can be orders of magnitude more efficient than random data access or access to small files. In the Kaldi toolkit [4], we try very hard to ensure that any high-volume data access takes the form of sequential reads or writes on large files.
For neural network training, we keep data access sequential by dumping pre-randomized "training examples" to disk. Each training example corresponds to a class label together with the corresponding input features, including left and right temporal context as needed by the network. The randomization is done just once for the entire data, and the data is accessed in the same order on each epoch. This is probably not ideal from the point of view of the convergence of SGD, but our expectation is that for large amounts of data the same-order access will not affect the results noticeably.
We break up the training data into N by M rougly equal-sized blocks, where N is the number of parallel jobs, specified by the user (typically 4 ≤ N ≤ 8), and M ≥ 1 is the number of "outer iterations per epoch", which is chosen to ensure that the number of samples processed per iteration is close to a userspecified value K (e.g. K = 400 000). The process of randomly distributing the data into N by M blocks, and ensuring that the order is randomized within each block, is done in parallel; we won't give further details here, because the problem is straightforward and there is nothing particularly special about our method.

Other aspects of our SGD implementation
At this point we provide some more details of other relevant features of our SGD implementation, namely the learning rate schedule and the way we enforce a maximum parameter change to prevent divergence.

Learning rate schedule
It was found in [5] that when training DNNs for speech recognition, an exponentially decreasing learning rate works well, and we independently found the same thing. We generally use a learning rate that decreases by a factor of 10 during training, on an exponential schedule. Unless mentioned otherwise, for experiments reported here the learning rate starts at 0.01 and ends at 0.0001. We specify the number of epochs in advance; it is typically a number in the range 4 to 20 (if we have more data, we train for fewer epochs).

Maximum parameter change
A common pathology when doing SGD for deep learning is that during training, the parameters will suddenly start getting very large and the objective function will go to negative infinity. This is known as parameter divergence. The normal solution is to decrease the learning rate and start the training again, but this is a very inconvenient. To avoid this pathology, we modified the SGD procedure to enforce a maximum parameter change per minibatch. This limit tends to be active only early in training, particularly for layers closer to the output. We have provided further details on this in Appendix C.1.

Natural gradient for SGD
In this section we describe our natural-gradient modification to SGD, in which we scale the gradients by a symmetric positive definite matrix that is an approximation to the inverse of the Fisher matrix.
Technically speaking, Natural Gradient means taking a step along the gradient of a Riemannian parameter surface, which follows a curving path in conventional parameter space and which is extremely hard to compute. However, previous work [6,7] has used the term "Natural Gradient" to describe methods like ours which use an an approximated inverse-Fisher matrix as the learning rate matrix, so we follow their precedent in calling our method "Natural Gradient".

We can replace the scalar learning rate in SGD with a matrix
In SGD, the learning rate is often assumed to be a scalar η t , decreasing with time, and the update equation is something like where g t is the objective-function gradient sampled on time t (e.g. computed from a training sample or a minibatch). However, it is possible to replace this scalar with a symmetric positive definite matrix, and we can write instead: with E t the matrix component of learning-rate; it is more convenient for proofs to keep η t separate rather than absorbing it into E t . It acceptable for E t to be random: if we can bound the eigenvalues of E t above and below, by positive constants known in advance, and E t and g t are independently sampled given the parameter θ, then we can prove convergence under the same kinds of conditions as if we were using a scalar learning rate [8,Sec. 4.2.2] In general, the learning-rate matrix should not be a function of the data sample which we are currently processing, or it may prevent convergence to a local optimum. As an example of this, a matrix that was systematically smaller for a particular type of training data would clearly bias the learning by downweighting that data.

The inverse Fisher matrix is a suitable learning-rate matrix
There are reasons from statistical learning theory, related to the Natural Gradient idea [9], why we may want to set the learning-rate matrix E t to the inverse of the Fisher information matrix. See for example, [10] and [7]. The Fisher matrix is most directly defined for situations where we are learning a distribution, as opposed to classification problems such as the current one. Suppose x, which may be discrete or continuous, is the variable whose distribution we are modeling, and f (x; θ) is the probability or likelihood of x given parameters θ, then the Fisher information matrix I(θ) is defined as the expected outer product (second moment) of the derivative of the log-probability w.r.t. the parameters, i.e. of This derivative is called the "score" in information theory. Part of the justification for this use of the Fisher matrix is that, under certain conditions, the Fisher matrix is identical to the Hessian; and it is obvious why the inverse Hessian would be a good gradient descent direction. These conditions are quite stringent, and include that the model is correct and θ is at the value corresponding to the true data distribution; but even if these conditions do not apply, the Fisher information matrix is in some sense "dimensionally" the same as the Hessian-that is, it transforms the same way under changes of parameterization-so its inverse may still be a good choice of learning-rate matrix.
It is quite easy to generalize the notion of the Fisher matrix to a prediction task p(y; x, θ). We write p(y, x; θ) = q(x)p(y; x, θ) for a data distribution q(x) that we assume to be independently known (and not a function of θ). It is not hard to see that the score equals just ∂ ∂θ log f (x; y, θ); since q(x) does not depend on θ, there is no additional term involving q(x). The expectation that we take when computing the Fisher matrix is taken over the joint distribution of x and y.
Still more generally, we can compute a quantity that is analogous to Fisher matrix for any objective function, even one that does not represent a log-probability or log-likelihood; and we will still have a matrix that transforms in the same way as the Hessian under changes of variables -i.e. its inverse may still be a reasonable choice for a learning-rate matrix.

We need to approximate the Fisher matrix in practice
For large-scale problems, such as DNNs for speech recognition with millions of parameters, even one inversion of the Fisher matrix is impractical because it would take time O(n 3 ) in the parameter dimension. However, it may be practical to deal with factored forms of it. There has been previous literature on this. In [7], the Fisher matrix was divided into diagonal blocks and each block was approximated by a low-rank matrix. The idea of diagonal blocks was also explored in [11], with one block per weight matrix; our approach uses the same idea. In the unpublished manuscript [12] (some of the material in which was later published in [6]), the authors attempted to show analytically that under certain quite strong assumptions, the Fisher matrix for a single-hidden-layer neural network has the form of a Kronecker product.
Although we are skeptical of this, and anyway we are interested in more general networks than they considered, the Kronecker product does also appear in our factorization of the Fisher matrix.
We should note that there are ways to use Natural Gradient without factorizing the Fisher information ma-trix, if one is willing to accept a significantly increased time per iteration. See for example [13], which uses a truncated Newton method to approximate multiplication by the inverse of the Fisher matrix.

Our factorization of the Fisher matrix
Our factored form of the Fisher matrix is as follows: given a neural network with I weight matrices, we divide the Fisher matrix into I diagonal blocks, one for each weight matrix. Consider the i'th diagonal block of the Fisher matrix, corresponding to the parameters of a weight matrix W i , and assume that there is no separate bias term (we can append a 1 to the inputs and include the bias term in the weight matrix). The i'th block of the Fisher matrix is a Kronecker product of two symmetric positive definite matrices: A i , whose dimension is the input (row) dimension of W i , and B i , whose dimension is the output (column) dimension of W i . We further factorize the matrices A i and B i as a low-rank symmetric matrix plus a multiple of the identity matrix. We write the approximated Fisher matrix in the form where A i and B i are each factorized in the form λI + XX T . The order in which A i and B i appear in the Kronecker product depends on the way in which we vectorize the weight matrices-row-wise or columnwise. In practice we don't ever deal explicitly with these Kronecker products or vectorized weight matrices in the algorithm, so this choice doesn't matter. It is not hard to show that if the Fisher matrix can be factored this way, then its inverse can be factored the same way.

How we estimate the Fisher matrix
We have two different methods for estimating the factorized Fisher matrix: • Simple method: We estimate the Fisher matrix from the other samples in the minibatch we are currently training on, holding out the current sample to avoid bias. This can be done surprisingly efficiently. Details are in Appendix A.
• Online method: We estimate the Fisher matrix from all previous minibatches, using a forgetting factor to downweight minibatches that are distant in time. Details are in Appendix B.
We generally use the online method as it is significantly faster on GPUs and usually seems to lead to faster learning, probably due to the less noisy estimate of the Fisher matrix. We describe the simple method because it is easier to understand and helps to motivate the online method.

Operation on vectors
Although we describe our Fisher matrix as as a Kronecker product, we do not have to explicitly construct this product in our code.
Suppose that we process the training examples one at a time. The SGD update for the i'th weight matrix is: where x ti is the input to the i'th weight matrix computed at the current sample, and y ti is the derivative of the objective function with respect to the output of the i'th weight matrix. These quantities naturally occur in backpropagation.
In our natural gradient method, this is modified as follows: where A ti and B ti are factors of the Fisher matrix. It is easy to show that this is equivalent to multiplying the parameter step by the inverse of the Fisher matrix formed from the A and B quantities as in Equation (2).

Operation on minibatches
Rather than processing training examples one at a time, we process them in minibatches (e.g. 512 at a time). Instead of vector-valued derivatives x ti and y ti , we now have matrices X ti and Y ti , each row of which corresponds to one of the x or y quantities (t is now the index for the minibatch). The update is now as follows: and note that unlike some authors, we don't divide the gradients by the minibatch size-this makes it easier to tune the minibatch size and learning rate independently. The update now has the form with the bar indicating modified X and Y quantities.
In the online version of our natural gradient method, we can write these as: but in the simple method, because the A and B matrices are estimated from the other elements in the minibatch, we can't write it this way-it is a separate matrix multiplication for each row of X and Y-but it can still be computed efficiently; see Appendix A.
In programming terms, we can describe the interface of the core natural-gradient code as follows: • Simple method: Given a minibatch of vectors X ti with each row being one element of the minibatch, estimate the Fisher-matrix factors by holding out each sample, do the multiplication by their inverses, and return the modified vectorsX ti .
• Online method: Given a minibatch of vectors X ti and a previous Fisher-matrix factor A t−1,i , com-puteX ti = X ti A −1 t−1,i and the updated Fishermatrix factor A ti .
The interface of the natural gradient code works the same with the Y and B quantities, as with X and A. We call the interface above 2I times for each minibatch: twice for each weight matrix in the network.

Scaling the factors
In both natural gradient methods, we want to prevent the Fisher-matrix multiplication from affecting the overall magnitude of the update very much, compared with the step-sizes in standard SGD. There are several reasons for this: • Early in training, the x and y quantities may be very small or zero, leading to huge or infinite inverse-Fisher matrices.
• The conventional convergence-proof techniques require that the matrix component of the learning rate matrix should have eigenvalues bounded above and below by constants known in advance, which we cannot guarantee if we use an unmodified Fisher matrix.
• Empirically, we have found that it is hard to prevent parameter divergence if we use the real, unscaled Fisher matrix.
Our method is to scale theX ti andȲ ti quantities so that they have the same Frobenius norm as the corresponding inputs X ti and Y ti . We will introduce notation for this in the Appendices.
This scaling introduces a slight problem for convergence proofs. The issue is that each sample can now affect the value of its own learning-rate matrix (via the scalar factor that we use to rescale the matrices). As we mentioned before, it is not permissible in general to use a per-sample learning rate that is a function of the sample itself. However, we don't view this as a practical problem because we never use a minibatch size less than 100, so the resulting bias is tiny.

Smoothing the Fisher matrix factors with the identity
In both versions of NG-SGD, we smooth our estimates of the factors of the Fisher matrix by adding a multiple of the identity matrix before inverting them. In the simple method this is necessary because in general the Fisher matrix estimated from the minibatch will not be full rank. In the online method it is not strictly necessary because we deal with a factorization of the Fisher matrix that already contains a multiple of the unit matrix, but we found that by adding an additional multiple of the unit matrix, as for the simple method, we can improve the convergence of the SGD training.
In both cases the smoothing is of the following form. If S ∈ R D×D is a Fisher matrix factor estimated directly from data as the uncentered covariance of the x or y quantities, then instead of using S as the Fisher-matrix factor A or B, we use instead S + βI, where where ǫ = 10 −20 is used to stop the smoothed S from ever being exactly zero. That is, we smooth the Fisher with the identity matrix scaled by α times the average diagonal element of S. We found in tuning experiments that the relatively large value α = 4 is suitable under a wide range of circumstances, for both the simple and online methods, and even for settings where the noise in S should not be a big problem-e.g. for large minibatch sizes. Our interpretation is that when α is fairly large, we are using a smaller than normal learning rate only in a few directions where the x or y quantities have quite high covariance, and a relatively constant learning rate in all the remaining directions.

What we think we can, and can't, prove
Although this is not a theoretical paper, we would like to say what we think is, and is not, possible to prove about our methods.

Our factorization of the Fisher matrix
If we assume that the distribution of the x and y quantities is Gaussian and independent (between x and y for a single layer, and between layers), then it should not be hard to show that the Fisher matrix has the form of (2), where the A i and B i quantities correspond to the uncentered covariances of the x and y quantities, and that the inverse-Fisher has the same We don't believe that the Fisher matrix would have this exact form in practice, but we do believe that it's a reasonable factorization. One could try to show this experimentally as follows, for a small task. One could make a linear change of variables to make our approximated Fisher matrix equal the unit matrix, and then try to measure the eigenvalue distribution of the full Fisher matrix in the new co-ordinates. We believe that the eigenvalue distribution of the transformed Fisher matrix would probably be much more closerly centered around 1 than before the change of variables. Our motivation for this work is a practical one, so we have not allocated effort towards this type of experiment.

The convergence of our preconditioned SGD procedure
Regarding the convergence of SGD using our factored-Fisher learning rate matrices, the most we think is provable is that a slightly modified form of this method would converge under similar conditions to unmodified SGD. The smoothing with constant α > 0 can give us a bound on the ratio of the largest to smallest eigenvalues of the A and B factors; using this together with the rescaling of Section 4.8, we can bound from above and below the eigenvalues of the rescaled A and B factors. By multiplying these together, we can get lower and upper bounds on the eigenvalues of the overall inverse-Fisher matrix that we use as the learning-rate matrix E t .
It is necessary for the Fisher matrix to be randomly chosen independent of the identity of the current sample. Unfortunately this is not quite true due to the rescaling being done at the minibatch level; we mentioned in Section 4.8 that this would be a problem for proofs. As mentioned, it would be easy to use the rescaling factor from the previous minibatch; this gives us back the independence, but at the cost of no longer having such easy bounds on the upper and lower eigenvalues of the rescaled A and B factors. Alternately, one could keep the algorithm as it is and try to prove instead that the parameter value we converge to will not differ very much in some sense from an optimum of the true objective function, as the minibatch size gets large.

Online update of a low-rank covariance matrix
There might be some interesting things to say about our online natural gradient method, described in Appendix B, in which estimate the uncentered covariance matrices A and B in a factored form as λI + XX T .
Our online estimation of the covariance matrices involves multiplying X by a weighted combination of (a) the observed covariance matrix from the current minibatch, and (b) the previous value of our factored approximation to it; it is like a matrix version of the power method [14].
Probably the analysis would have to be done initially in the steady state (i.e. assuming the parameter vector θ is constant). If in addition we assume infinite minibatch size so that the covariance matrix equals its expected value, we are confident that we could show that the only stable fixed point of our update equations gives us in some suitable sense the closest approximation to the covariance; and, with a little more effort, that our updates will converge with probability 1 to that best approximation.
The analysis for finite minibatch size would have to involve different methods. Because of the noise and the finite forgetting factor, we would never converge to the true value; but it might be possible to define some objective function that measures some kind of goodness of approximation, and then say something about the convergence of the distribution of that objective function.

Experiments
We show experiments on a speech recognition setup called Fisher English 1 , which is English-language conversational telephone speech, sampled at 8kHz, and transcribed in a quick but relatively low-quality way. The total amount of training data is 1600 hours (only including transcribed segments, i.e. not the silent other half of the telephone conversation). We test on a held-out subset of the data, about 3.3 hours long, that we defined ourselves. Error Rates may seem on the high side, but this is mainly due to the difficulty of the data and the quick transcription method used on this data.

System details and Word Error Rate performance
The GMM system is based on MFCC features, spliced across ±3 frames and processed with LDA+MLLT to 40-dimensional features, then adapted with featurespace MLLR (fMLLR) in both training and test time. See [4] for an explanation of these terms and the normal system build steps. All these systems used the same phonetic context decision tree with 7 880 contextdependent states; the GMM system had 300 000 Gaussians in total.
The DNN1 system uses speaker adapted features from the GMM system, so it requires a first pass of GMM decoding and adaptation. The 40-dimensional features from GMM1 are spliced across ±4 frames of context and used as input to the DNN. DNN1 is a p-norm DNN [15] with 5 hidden layers and p-norm (input, output) dimensions of (5000, 500) respectively, i.e. the nonlinearity reduces the dimension tenfold. We use 15 000 "sub-classes" (see Section C.3 for explanation), and the number of parameters is 19.3 million. It is trained for 12 epochs with learning rate varying from 0.08 to 0.008, trained with 8 parallel jobs with online natural gradient SGD (NG-SGD). For both this and the DNN2 system, we trained with K = 400 000 samples per outer iteration for each machine.
The DNN2 system is trained for our online decoding setup (see Appendix C.7), which is geared towards applications where reduced latency is important and audio data must be processed strictly in the order it is received. The input features are equivalent to unadapted, un-normalized 40-dimensional log-mel filterbank features, spliced for ±7 frames, plus a 100dimensional iVector representing speaker characteristics, extracted from only the speaker's audio up to and including the current time. For the results shown here, we include previous utterances of the same speaker in the same conversation when computing the iVector. Because this system is intended for real-time decoding on a single CPU, we limit the number of parameters by using only 4 hidden layers, p-norm (input, output) dimensions of (350, 3500), and 12 000 sub-classes, for a total of 10.4 million parameters. It was trained using online NG-SGD with 6 parallel jobs for 5 epochs, with the learning rate decreasing exponentially from 0.01 to 0.001. All our experiments below are based on this setup.
Our server hardware is fairly typical: the majority of them are Dell PowerEdge R720 servers with two Intel Xeon E5-2680v2 CPUs having with 10 cores each, running at 2.8GHz; and with a single NVidia Tesla K10 GPU card, providing two GPUs-each GPU cor-responds to a single machine in our notation, and it becomes incidental that they are co-located. We also have some similar machines with K20 GPU cards, and when reporting time taken, we report the slightly more optimistic figures obtained from running the same jobs on the faster K20 GPUs.

Results
Our main result is in Figure 1 (best viewed in color), where we plot the objective function versus amount of training data processed, for our parallel training method with and without natural gradient, and with 1, 2, 4, 8 and 16 jobs. In order to keep the effective learning rate (Section 3.1) constant, we make the initial/final learning rates proportional to the number of jobs, with the default learning rates of 0.01 to 0.001 corresponding to the 6-job case.
Our natural gradient method always helps-the NG-SGD curves are all above the plain-SGD curves. Also, when using online natural-gradient, the curves shown in Figure 1 are close to each other up to about 4 jobsi.e. after processing the same amount of data with different numbers of jobs we get about the same objective function; however, the 8-and 16-job runs converge a little slower. Thus, for small N we are getting a linear speed up in the number N of machines, because the time taken per epoch is proportional to 1/N . As N gets larger than around 4 we need more epochs to get the same improvement, so the speedup becomes sublinear. The plot also shows that the simple and online natural gradient converge about the same (only tested with one job). We show the final Word Error Rates in Table 2; with NG-SGD, they are not very sensitive to the number of jobs. Figure 2 shows the same plots as Figure 1 but with time as the x-axis. This is a simulated clock time, obtained by multiplying the time taken for each "outer iteration" of training, by the number of outer iterations; the actual clock time depends on queue load. The time per outer iteration was 88 seconds for plain SGD, 93 seconds for online NG-SGD, and 208 seconds for plain NG-SGD, all measured on a K20 GPU. The circles mark the end of training, after 5 epochs.

Conclusions
We have introduced an efficient Natural Gradient version of SGD training (NG-SGD). We have shown experimentally that not only does the method improve the convergence versus plain SGD, it also makes it possible for us to to use very simple parallelization method where we periodically average parameters from multiple SGD runs. Although we only show results from one setup, we are confident based on past experience that it holds true for other types of neural network and improves our final results (Word Error Rate) as well as convergence speed.
We do not have a very good explanation why our parallel training method only works when using Natural Gradient, except to say that the statements in [13] that NG prevents large parameter steps and is more robust to reorderings of the training set, may be relevant.
The authors were supported by DARPA BOLT contract No HR0011-12-C-0015, and IARPA BABEL contract No W911NF-12-C-0015. We gratefully acknowledge the support of Cisco Systems, inc. (grant #574560) and Google, Inc. (Award 2012 R2 106, "Deep Neural Networks for Speech Recognition"), funds which were used to buy computer equipment and cloud computing time that were used in the development of these methods.
The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of DARPA, IARPA, DoD/ARL or the U.S. Government.
[5] Andrew Senior, Georg Heigold, Marc'Aurelio Ranzato, and Ke Yang, "An empirical study of learning rates in deep neural networks for speech recognition," in Acoustics, Speech and Signal Processing (ICASSP), 2013.
A Further details on simple natural gradient method

A.1 Overview of simple natural gradient method
In this section we describe the natural gradient method that uses the other elements of the minibatch to estimate the factors of the Fisher matrix.
As mentioned in Section 4.7, the interface can be described as follows. Given a matrix X, each row of which represents one element of the minibatch (and with a number of columns corresponding to either the row or column dimension of one of the weight matrices in the network), do the inverse-Fisher multiplication for each row x i of X and return the modified matrix X.
The core of the inverse-Fisher multiplication is this: let We extend this basic idea by adding smoothing of F i with the identity matrix, and by scaling the outputX to have the same Frobenius norm as the input.

A.2 Details of method (not considering efficiency)
In this section we describe what we compute in our "simple" natural gradient method, without considering how to compute it efficiently. As described in Section 4.9, we smooth the Fisher matrix with the identity. Defining where our normal settings are α = 4 and ǫ = 10 −20 , and N and D are the number of rows and columns respectively of X, we define smoothed Fisher matrices as follows, to be applied to each row x i of X: For each row i we will then definê and then the result of our computation will bē where the rescaling factor γ, intended to make sure thatX has the same Frobenius norm as the input X, is defined as If the denominator of the above is zero, we take γ to be one.
We should note that by computing the scalars β and γ without "holding out" the current sample, we are violating the rule that the randomly sampled learning rate matrix should be independent of the current sample. However, since we always use a fairly large minibatch size (at least 100) and these are only scalar quantities, we don't believe the small amount of "contamination" that takes place here will significantly bias the training. In fact, it might not turn out to be very difficult to modify the equations to properly hold out the current sample for these purposes, but because we don't believe it would perceptibly affect the results, we haven't gone to the trouble of doing this.

A.3 Efficient computation in simple method
We now describe how we efficiently compute what we described above. Define the smoothed Fisher matrix which is like the G i quantities but without holding out the current sample. Next, compute where Q only differs fromX by x i not being held out from the corresonding G i . There are two equivalent methods to compute Q: (i) In column space: (ii) In row space: We derived the rather surprising row-space version of the formulation by expanding the inverted expression on the right of the column-space expression using the Morrison-Woodbury formula, and simplifying the resulting expression.
For efficiency, we choose method (i) above if the minibatch size is greater than the dimension (N > D), and method (ii) otherwise. Our formula below forX is derived by expressing each G −1 i as a rank-one correction to G −1 , and computing the corresponding correction by which each rowx i differs from the corresponding row q i of Q. It turns out that the correction is in the same direction as q i itself, sox i just becomes a scalar multiple of q i . Defining, for each row-index i, and defining the scalar factor then we can use the following efficient formula: for each rowx i ofX,x We then get the outputX by scalingX by γ, as described above.
When working on CPUs with small minibatch sizes (e.g. N = 128) and large hidden-layer dimensions (e.g. D = 1000), the computation above is very efficient, and does not comprise more than about 20% of the time of the overall backprop computation. However, when using GPUs with larger minibatch sizes (e.g. N = 512) it can take the majority of the time. Even though it typically takes considerably less than half of the total floating point operations of the overall computation, it contains a matrix inversion, and matrix inversions are not very easy to compute on a GPU. Our "online" method which we will describe below is designed to solve this efficiency problem.

B.1 Overview of online natural gradient method
The interface of the online natural-gradient method is essentially the same as the simple method: the user provides a matrix X, and we return a matrixX that's been multiplied by the inverse-Fisher and then rescaled to have the same Frobenius norm as X. Again, each row of X corresponds to an element of the minibatch and the column dimension corresponds to the row or column dimension of one of the weight matrices. A difference from the simple method is that the online method is "stateful", because we maintain a running estimate of the Fisher matrix. Each time we process a minibatch, we use the Fisher matrix estimated from the previous minibatches; and we then update that estimate using the current minibatch. For a single neural net, the number of separate copies of this "state" that we need to maintain corresponds to twice the number of trainable weight matrices in the neural net: one for each of the A i and B i quantities in Equation (2).
Let the input be X ∈ R N ×D , where N is the minibatch size (e.g. 512) and D is the row or column size of the weight matrix we're updating (e.g. 2000). We introduce a user-specified parameter R < D which is the rank of the non-identity part of the Fisher matrix. Let the subscript t = 0, 1, . . . correspond to the minibatch.
where R t ∈ R R×D , D t ∈ R R×R and ρ t > 0 will be estimated online from data; R t has orthonormal rows and D t is diagonal and nonnegative. We'll estimate these quantities online from the data with the aim being that F t should be a good estimate of the covariance of the rows of the X t quantities.
We define G t to be a kind of "smoothed" version of F t where we add in more of the unit matrix, controlled by the α parameter we've previously discussed (normally α = 4): and then the output will be: where γ t is computed so as to ensure that the Frobenius norm ofX t equals that of X t : or γ t = 1 if the denominator of the above equation is 0.

B.2 Updating our low-rank approximation to the variance
Next we discuss the method we use to estimate our low-rank approximation F t of the uncentered covariance of the rows of the inputs X t . Define as the uncentered covariance of the rows of X t . We introduce a user-specified "forgetting factor" 0 < η < 1 (we describe how this is set in Section B.4), and we define We will try to set F t+1 to be a good low-rank approximation to T t . The obvious way would be to make D t+1 correspond to the top eigenvalues of T t and R t+1 to the corresponding eigenvectors, but this would be too slow. Instead we use a method inspired by the power method for finding the top eigenvalue of a matrix. On each iteration we compute with Y t ∈ R R×D . It is useful to think of Y t as containing each eigenvector scaled by its corresponding eigenvalue in T t (of course, this is true in a precise sense only at convergence). Our update uses symmetric eigenvalue decomposition of Y t Y T t to find these scaling factors (they are actually the square roots of the eigenvalues of Y t Y T t ), puts them on the diagonal of D t+1 , and puts the corresponding eigenvectors in the rows of R t+1 . We then have to work out the correct amount of the unit-matrix to add into our factorization of the covariance matrix (i.e. set ρ t+1 ) and subtract that amount from the diagonals of D t+1 . We will give equations for this below.
Observant readers might have noted that it would seem more straightforward do do a Singular Value Decomposition (SVD) on Y t instead of a symmetric eigenvalue decomposition on Y t Y T t . We do it this way for speed.
The details of our update are as follows: so Z t ∈ R R×R . Then do the symmetric eigenvalue decomposition with U orthogonal and C t diagonal. The diagonal elements of C t are positive; we can prove this using ρ t > 0 (which makes T t positive definite) and using the fact that R t has full row rank. We define R t+1 as: If we expand out R t+1 R T t+1 using (27), it is easy to see that it reduces to the identity, hence R t+1 has orthonormal rows. In order to make sure that F t+1 has the desired covariance in the directions corresponding to the rows of R t+1 , we will set but note that at this point, ρ t+1 is still unknown. When we say the "desired covariance", we are ensuring that for each dimension r corresponding to a row of R t+1 , the value of the inner product r T F t+1 r equals that of r T T t r, but this is only precisely true at convergence.
We choose ρ t+1 in order to ensure that tr (F t+1 ) = tr(T t ). This value can be worked out as: We then let for ǫ = 10 −10 ; this is is to ensure that if we get a sequence of zero inputs, ρ t will not become exactly zero in its machine representation.

B.3 Efficient computation
The previous section described what we are computing in the online natural gradient method; here we describe how to compute it efficiently. The essential idea here is to reduce the multiplication by G −1 to two multiplications by a "fat" matrix (of dimension R×D).
Since typically R is much smaller than D, this is quite efficient. We also address how to efficiently keep these matrices updated, at the level of optimizing the matrix expressions. This section is mostly derivation, and will likely only be of interest to someone who is considering implementing this method. In Section B.5 below, we will summarize the algorithm we derive here.
We can write G t as: where the factor of β t is inserted arbitrarily to simplify the update equations; a scalar factor onX doesn't matter because we will later rescale it to have the same norm as X. The output of this whole process is where, in the expression for γ t , if the denominator is zero we take γ t = 1. Note: γ t is not the same as in (21) because of the arbitrary factor of β t , so consider (21) to be superseded by (37). To efficiently compute (35), we apply the Woodbury matrix identity to (31), giving us with elements e tii = 1 β t /d tii + 1 (40) In order to reduce the number of matrix multiplies, it is useful to break the expression R T t E t R t into two equal parts, so we define and we will never store R t ; instead, we will work with W t and the small diagonal factors D t and E t . We can now write the following, which is where most of our computation will take place: You may recall the symmetric matrix Z t ∈ R R×R defined in (25), which is involved in the update of our factorization. The following expressions are going to be useful when computing it, and the first of them appears as a sub-expression of (42). For convenience we state the dimensions of these quantities below: After we have H t , we can computeX t using a single matrix multiply as: We can expand Y t = R t T t , defined in (24), into quantities that will be computed, as: Using (52) we can expand Z t = Y t Y T t , as: and we can substitute some of the sub-expressions we defined above into this, to give: Our strategy will be to compute the symmetric quantities L t and K t on the GPU, and transfer them to the CPU where we can then compute Z t using the expression above -this can be done in O(R 2 ) -and then do the symmetric eigenvalue decomposition as in (26), on the CPU. We repeat the equation here for convenience: Here, U t will be orthogonal, and mathematically, no element of the diagonal matrix C t can be less than (1−η) 2 ρ 2 t , so we floor its diagonal to that value to prevent problems later if, due to roundoff, any element is smaller than that.
Below, we'll say how we efficiently compute tr (XX T ) and tr (XX T ); for now, just assume those quantities have been computed.
We compute ρ t+1 as follows, expanding S t in (29): . (56) We can now compute D t+1 and ρ t+1 ; we floor both to ǫ to ensure they never go to exactly zero which could cause problems for our algorithm.
for a small constant ǫ = 10 −10 (the first max is taken per element). We can now compute the scalar β t+1 and the diagonal matrix E t+1 (we show the formula for its diagonal elements): We never construct R t+1 in memory, but instead we directly compute W t+1 . We can factor it as follows: and note that while it might seem like a factor of E −0.5 t is missing from the second term in B t , in fact we use the fact that it commutes with (D t + ρ t I) to move it to the left, into A t . If we're using a GPU, A t will be computed in time O(R 2 ) on the CPU and transferred to the GPU; we then compute B t on the GPU efficiently by scaling the rows of W t and adding J t ; then we multiply A t and B t on the GPU.

B.3.1 Maintaining orthogonality
We have noticed that the invariance R t R T t = I can sometimes be lost due to roundoff. A proper analysis of roundoff in our algorithm is not something we have time to do, but we will describe how we detect and fix this problem in practice. For speed, we only do the following operations if the diagonal matrix C t , has condition number greater than 10 6 , or if any elements were floored as mentioned just after (55). Note: all the computations we describe in this paper were done in single precision.
We compute the symmetric matrix where the part in parentheses is computed on the GPU and transferred to the CPU. If no element of O t differs by more than 10 −3 from the corresponding element of the unit matrix, we consider that R t is sufficiently orthogonal and we do nothing more. Otherwise, we do a Cholesky decomposition O t = CC T , compute the reorthogonalizing factor M = E 0.5 t C −1 E −0.5 t on the CPU and copy to the GPU, and do W t+1 ← MW t+1 to reorthogonalize. Re-orthogonalization happens extremely rarely, and usually only if something bad has already happened such as parameter divergence.

B.3.2 Initialization
In our implementation we don't bother dumping the "state' of the computation to disk so each new process reinitializes them for the first minibatch it processes. We initialize them so as to most closely approximate the covariance of the first minibatch of features. This is done by taking and finding the top R eigenvalues and eigenvectors; the rows of R 0 contain the top eigenvectors. Let λ i be the corresponding eigenvalues, for 1 ≤ i ≤ R, and we set for ǫ = 10 −10 , and for 1 ≤ i ≤ R, we let d 0ii ← max(ǫ, λ i − ρ 0 ).

B.3.3 Computing matrix traces
We mentioned above that we have a fast way of computing the quantities tr (XX T ) and tr (XX T ). These are needed to compute γ t using (37), and to compute ρ ′ t+1 using (56). We compute these as a side effect of the fact that we need, for each rowx ti of the output, its squared normx T tix ti . This will be required to enforce the "maximum parameter change" per minibatch, as described in Section 3.4.2. Suppose we've already com-putedX t using (50). We compute the inner products for all rows 1 ≤ i ≤ N ofX t as using a single GPU kernel invocation. If we are updating the parameters of the Fisher-matrix factorization, then we can most efficiently obtain our desired traces as follows: tr (XX T ) = tr (XX T ) − tr (L t E t ) + 2tr (L t ). (73) The expression for tr (XX T ) was obtained by expanding tr (XX T ) using (50), moving tr (XX T ) to the left, and recognizing sub-expressions that we have already computed. In case we are not updating the parameters of the Fisher-matrix factorization, we have no other need for L t so it will be more efficient to compute tr (XX T ) directly; this can of course be done in O(N D) operations and does not require a matrix multiply. Once we have the scaling factor γ t we can scale the p i quantities by its square, and they will equal the quantitiesx T tix ti that we'll need for enforcing the maximum parameter change.

B.3.4 Multithreading and other issues
Most of what we have written above is geared towards operation using a GPU, but we also support operation with CPUs, where our SGD implementation is multithreaded. In this case, we have to consider the interaction with multithreaded code because of the "stateful" nature of the computation. We wanted to avoid a bottleneck where different threads wait to update the parameters sequentially. Our solution is that before doing the part of the computation where we update the parameters, we try to get a lock, and if this fails, we simply apply the fixed preconditioning matrix but don't update the Fisher-matrix parameters R t and so on. Since the model parameters don't move very fast, we don't expect that this will make any noticeable difference to the SGD convergence, and we have seen no evidence that it does.

B.4 Typical configuration
The most important user-specified parameters for our algorithm are the rank R and the constant α that controls smoothing with the unit matrix. The value α = 4 seems to work well over a wide variety of conditions, so we normally leave it at that value. The rank R should generally increase with the dimension of the vectors we are multiplying. Our experiments here are with "p-norm" networks where the nonlinearity is dimension reducing, like maxout [16], typically reducing the dimension from something like 3000 to 300. So a typical parameter matrix will increase the dimension from something like 301 to 3000 (it's 301 instead of 300 because of the bias term). Our normal rule for ranks is to use R = 20 on the input side of each matrix and R = 80 on the output side. Part of the way we originally tuned this is to look at the diagonal matrices E t . These matrices have diagonal values 0 < e tii < 1, sorted on i from greatest to least, and 1 − e tii can be interpreted as the amount by which the input is scaled in a certain direction in the space. A value of e tii close to 1 means we are strongly scaling down the input, and a value close to 0 means we are leaving it unchanged. If the last e tii has a value of, say, 0.1, then reducing R by one will be like taking a scaling factor of 0.9 applied to a gradient, and setting to 1 instead; this seems unlikely to make any difference to the SGD, as it's like changing the learning rate in some direction from 0.9 to 1. Our final e tii values are normally in the range 0.05 to 0.2.
Another configurable constant is the "forgetting factor" 0 < η < 1: the closer η is to 1, the more rapidly we track changes in the Fisher matrix due to changes in parameters, but the more noise we will have in our estimates. Because we don't want to have to tune η when we change the minibatch size, we set it as follows. The user specifies a parameter S (interpreted as an approximate number of samples to include in our estimate of the Fisher matrix), and we set where N is the minibatch size. We normally set S = 2000; we have no reason to believe that this is a very important parameter.
In order to increase the speed of the algorithm, we normally configure it so that we only actually update the parameters of the Fisher matrix every 4 minibatches, except on the first 10 minibatches in a process, when we always update them.

B.5 Summary of the online natural gradient method
Here we summarize the online natural-gradient SGD method-that is, we summarize the core part of the algorithm that takes a matrix X ∈ R N ×D , and outputs a matrixX ∈ R N ×D . To understand how this fits into the bigger picture of back-propagation and SGD, see Section 4.
For this summary we will ignore issues of multithread-ing. Our explanation here is just for one instance of the algorithm, corresponding to the row or column dimension of one of the weight matrices; if there are I weight matrices, there are 2I separate copies of the variables we describe here.
Typical configuration variables are as follows: α = 4, S = 2000 (this will determine η), rank R = 20 (or 80), ǫ = 10 −10 ; and let's define a variable J = 4 that dictates the period with which we update the Fishermatrix factors. Minibatch size N is normally 128 or 512.
On t = 0, before running the steps below we have to initialize the parameters as described in Section B.3.2. Note: while in Section B.3.2 we describe how to set ρ 0 , R 0 and D 0 , the variables which we actually store are ρ 0 , D 0 , and W 0 ; to compute W 0 we need Equations (34), (40) and (41).
We have an input X ∈ R N ×D , and despite the notation, we do not require that N be the same for all tsometimes the last minibatch we process has a smaller than normal size.
If t < 10 or J divides t exactly, then we will be updating the factored Fisher matrix; otherwise we just apply it and don't update. There are two slightly versions of the algorithm, depending whether we will be updating the Fisher matrix.
In either case, we first compute η from N and S using (74), and then compute From here the two cases begin to differ.
Without updating the Fisher matrix. If we won't be updating the Fisher matrix, then it's simpler. The input is X t . We first compute tr (X T t X). Then we computeX overwriting the input X t . Next, for each 1 ≤ i ≤ N we compute the row-products p i using (71), and compute tr (X TX ) as the sum of p i . Now we can compute γ t using (37). Next we scaleX t by γ t to produceX t . We also output for each i the quantity γ 2 t p i =x T tix ti , which is needed to enforce the "maximum parameter change per minibatch" constraint.
With updating the Fisher matrix. If we're updating the Fisher matrix, which we usually do every four steps, there are some more operations to do. First we compute Next we want to compute L t and K t . We actually have two separate strategies for this. If N > D (the minibatch size exceeds the vector dimension), we do: and in our implementation we combine these into one matrix operation by placing L and K, and W and J, next to each other in memory. Otherwise, we compute K t as above but L t using: At this point, if we're using a GPU, we transfer the symmetric matrices K t and L t to the CPU. We now compute some small derived quantities on the CPU: β t using (34) and E t using (40), as well as E 0.5 t and E −0.5 t ; E t is diagonal so this is not hard. At this point we compute the symmetric R×R matrix Z t using (54); the expression looks scary but it can be computed in O(R 2 ) time.
We do the symmetric eigenvalue decomposition as in (55), on the CPU, to get the orthogonal matrix U t and the diagonal matrix C t , and we floor the diagonal elements of C t to (1−η) 2 ρ 2 t . Next we computeX then compute the row-products p i using (71), compute tr (X T tX ) = i p i , and compute tr (X T t X) using (73). We can now obtain the scaling factor γ t using (37), and use it to compute the main outputX t = γ tXt and the per-row inner products of the output which equal γ 2 t p i (although in our implementation, to save time we actually output γ t and let the user do the scaling later on).
We next compute ρ ′ t+1 using (56), D t+1 using (57) and ρ t+1 using (58). W t+1 is computed using a matrix multiply on the GPU as in (64), after working out the factors A t and B t .
At this point, if we had floored any diagonal elements of C t above or if its condition number after flooring exceeds 10 6 , we do the orthogonality check and possible reorthogonalization that we described in Section B.3.1 above.

C Other aspects of our DNN implementation
Here we describe some aspects of our neural net training implementation that are of less direct relevance to our parallel training and natural gradient methods, so were not included in the main text of the paper.
In Section C.1 we explain how we enforce a maximum parameter-change per minibatch; in C.2 we explain our generalized model-averaging procedure; in C.3 we explain how we use "mixture components" (a.k.a. subclasses) for DNNs; in C.4 we introduce our method of input data normalization; in C.5 we give details on how we initialize the DNN parameters; in C.6 we give an overview of how we implemented sequence training for DNNs; and in C.7 we discuss online (real-time) decoding using iVectors for speaker adaptation.
C.1 Enforcing the maximum parameter change per minibatch As mentioned in Section 3.4.2, in order to prevent instability and parameter divergence we enforce a maximum parameter-change per minibatch, which is applied for each layer of the network separately. Here we explain how this is done. We don't claim that this is an exceptionally good method for preventing excessive parameter changes, but we describe it here anyway for the sake of completeness.
Suppose the update for a single weight matrix is formulated as follows (and to keep things simple, we don't include an index for the layer of the network): where ∆ t is the change that standard SGD would give us, equal to the derivative of the objective function for this minibatch multiplied by the learning rate η t . To enforce the maximum parameter chanbge, we scale the change by a scalar α t : where we would like to choose α t ≤ 1 to ensure that ||α t ∆ t || F does not exceed a specified limit, || · || F being the Frobenius norm. However, we don't implement this scheme exactly as described above because it would involve creating a temporary matrix to store the product of matrices ∆ t just in order to compute its norm, and we don't want to incur this penalty.
Instead we enforce it in a way that involves a sum over elements of the minibatch. If ∆ t = ηX T Y, then ∆ t can be written as a sum over an index i that ranges over the rows of X and Y. By properties of norms, the 2-norm of ∆ t cannot exceed the sum of the 2-norms of the terms in this sum: if the rows of X and Y are written as x i and y i , then It does not take excessive time or memory to compute the vector norms ||x i || 2 and ||y i || 2 , so we compute the right hand side of 84 and use it as a stand-in for ||∆ t || F , giving us α t = min 1, max-change-per-minibatch where max-change-per-minibatch is a user-specified maximum parameter-change per minibatch. Empirically we have found that it tends to be necessary to increase max-change-per-minibatch when using a larger minibatch size, so to simplify the configuration process we define max-change-per-minibatch = N max-change-per-sample (86) where N is the minibatch size. We always set max-change-per-sample to 0.075 for experiments reported here. To clarify how this method interacts with the natural gradient methods described in Section 4: the natural gradient is implemented as a modification to the X and Y matrices, so we simply apply this maximum-change logic on top of the modified X and Y quantitities.
What we've found that this maximum-parameterchange limit is active only early in training for layers closer to the output.

C.2 Generalized model averaging
The convergence theory of Stochastic Gradient Descent [17] suggests that, for convex problems, if we take not the last iteration's model parameters but the average over all iterations, it can improve the convergence rate, particularly in 'poorly-conditioned' problems (i.e. where the condition number of the Hessian is very large). This is not applicable in non-convex problems such as ours, but it does suggest a related method. As mentioned above, we define an outer iteration as the length of time it takes for all jobs to process K samples (e.g. K = 400 000), and on each outer iteration each job dumps its final model to disk and we average these to produce a single model. We store the models (averaged over all jobs) for each outer iteration. At the very end of training, instead of choosing the model from the final outer iteration, we take the models from the last P outer iterations (e.g. P = 20), and search for a generalized weighted combination of these models that optimizes the objective function on a subset of training data-we tried using validation data here, but for our task we found it worked best to use training data. By generalized weighted combination, what we mean is that the parameters are a weighted combination of the parameters of the input models, but each layer can have different weighting factors. Thus, if there are P models and L layers, the number of parameters we learn on the data subdset is LP . A few more details: • The optimization method is L-BFGS.
• To improve the convergence speed of L-BFGS, we optimize in a transformed (preconditioned) space where the preconditioner is related to the Fisher matrix.
• The starting point for the optimization is the best of P + 1 choices, corresponding to each of the P final iterations, and the average of all of them.
• We add a very tiny regularizer (like 10 −10 times the square of the vector of weights) to stop the weights going to infinity in cases (like p-norm networks) where the objective function is invariant to the parameter scale.
• We generally aim to optimize over the last P = 20 models (assuming they share the same parameter structure, e.g. we haven't added layers).
• In cases where P iterations would amount to less than one epoch, we optimize over P models where the individual models are simple averages of model parameters for a duration of about 1/P of the entire epoch.
We have generally found that this model combination slightly improves results, but it is not a focus of the current paper so we don't provide experimental results for this here.

C.3 Mixture components (sub-classes)
When using Gaussians for speech recognition, the usual approach is to use a Gaussian mixture model (GMM) rather than a single Gaussian, to model each speech state. We have generalized this idea to neural networks, by allowing the posterior of each speech state to be written as a sum over the posterior of "sub-classes" that are analogous to the Gaussians in a GMM. About halfway through training, we "mix up" the model by increasing the dimension of the softmax layer to a user-specified number that is greater than the number of classes (usually about double the number of classes). After the softmax layer we introduce a "sum-group" layer which sums its input over fixed groups of indexes to produce a posterior for each class that is a sum over the posteriors of the hidden "subclasses". We also tried sharing the sub-classes across classes in groups, but did not find this helpful.
Rather than distributing the "sub-classes" evenly, we allocate more sub-classes to the more common classes. We allocate them proportional to the 1/3 power of the count of that class in the training data; this is based on the rule we use to allocate Gaussians in our GMMs.
When initializing the parameters of the "mixed-up" final weight matrix, we make it correspond quite closely with the original weight matrix. Each row of the new weight matrix corresponds to a row of the old weight matrix, plus a small noise term to allow the values of the rows to diverge; and we modify the bias term to normalize for the fact that some classes have more sub-classes than others.
We have generally found that this slightly improves results, but again, this is not a focus of the current paper and we won't be showing experimental results about this.

C.4 Input data normalization
As mentioned in [18, Section 4.3], when training neural networks it is helpful to normalize the input data so that it is zero mean and so that more important dimensions of input data have a larger variance. We wanted a generic way to achieve this that would be invariant to arbitrary affine transforms of the input. The technique we developed requires as statistics a within-class covariance W and a between-class covariance B, accumulated from the class-labeled data as if in preparation for multi-class Linear Discriminant Analysis (LDA). Assume in what follows that we have already normalized the data so that it is zero-mean. For the technique we are about to describe to make sense, the number of classes should not be much smaller than the feature dimension; fortunately, in our case it is much larger-5000 > 300, to give typical numbers.
Suppose we were to do multi-class LDA but not actually reduce the dimension. We would transform into a space where W was unit and B was diagonalized. Suppose the B in this space has diagonal elements b i . Then the total covariance in each dimension i is b i + 1. This has the desirable property that the data covariance is higher in "more important" directions, but it doesn't drop as fast as we'd like for unimportant directions-it never goes below 1. In our method, we do the LDA-type transform as mentioned above, then scale each row of the transform by (b i + 0.001)/(b i + 1). After this scaling, the total covariance becomes b i + 0.001, where b i is the ratio of between-class to within-class covariance. This seems to work well.
After creating the transform matrix as described above, we do a singular value decomposition on it, floor the singular values (normally to 5), and reconstruct again. The motivation here is to avoid a rarely encountered pathology that occurs when the training data covariance was close to singular, which leads to a transform with very large elements, that might produce very large transformed data values on mis-matched test data or due to roundoff. This step rarely floors more than a handful of singular values so has little effect on the transform.

C.5 Parameter initialization
We decided not to implement generative pre-training as in [19], because while it is well established that it improves results for small datasets, our understanding is that as the amount of training data gets larger, it eventually gives no improvement compared to a suitable random initialization or discriminative layerwise backpropagation as in [20]. We could not find a published reference for this; it is something we have been told verbally. We refer here specifically to speech recognition tasks; this does not apply to tasks like computer vision where much larger networks are used. In fact, the alternative "nnet1" implmentation of DNNs in Kaldi does support pre-training, and for small datasets (say, 50 hours or less), it generally gives slightly better results than the "nnet2" implementation which we speak of here. For larger datasets, the "nnet1" implementation eventually becomes impractical to run because it takes too long, and a detailed comparison is way beyond the scope of this paper.
Instead of pre-training, we use what is described in [20] as layer-wise back-propagation (BP). What this means is, we initialize a network with one hidden layer, train with BP for a short time (two "outer iterations" for our experiments reported here), then remove the final softmax layer and add a new, randomly initialized hidden layer on top of the existing hidden layer; train for a short time again; and repeat the process until we have the desired number of hidden layers. Similar to [21], we use a standard deviation of 1 √ i for the weights, where i is the fan-in to the weight matrix; but we initialize the parameters of softmax layers to zero. Note: we found it essential to discard the parameters of the final softmax layer when adding each new hidden layer, as prescribed in [20].
For smaller datasets we can improve results versus layer-wise BP by initializing all but the last layer of the network from a network trained on another large dataset, possibly from another language. When initializing this way we typically find it best to use a larger network than we otherwise would have used.
Because we noticed that sometimes on an outer iteration immediately following the random initialization of parameters (including the first outer iteration), the parameter averaging can degrade rather than improve the objective function, we modified our parallel training method so that on these iterations, instead of averaging the parameters we choose the one that had the best objective function computed on the subset of data that it was trained on (this choice of data avoids any extra computation).

C.6 Sequence training
Sequence training [22] is a term that has the the same meaning for DNNs that "discriminative training" [23] has in the speech recognition community for GMMs. It is a collective term for various objective functions used for training DNNs for sequence tasks, that only make sense at the whole-sequence level. This contrasts with the cross-entropy objective function which, given a fixed Viterbi alignment of the HMM states, easily decomposes over the frames of training data. In GMM-based speech recognition, the term "discriminative training" contrasts with Maximum Likelihood estimation; in DNN-based speech recognition it contrasts with cross-entropy training. There are two popular classes of sequence/discriminative objective functions: • Maximum Mutual Information (MMI)-like objective functions [24,23], more properly called conditional maximum likelihood: these have the form of sum over all utterances of the log-posterior of the correct word sequence for each utterance, given the model and the data. These include its popular 'boosted' variant [25] which is inspired by margin-based objective functions.
• Minimum Bayes Risk (MBR)-like objective functions: popular variants include Minimum Phone Error (MPE) [26,23] and state-level Minimum Bayes Risk [27,28]. These have the form of an expectation, given the data and the model, of an edit-distance type of error. We can compute its derivative w.r.t. the model parameters, because the posteriors of the different sequences vary with the model parameters.
This paper is mainly about our parallel approach to standard cross-entropy training. However, we also apply the same ideas (model averaging, preconditioning) to sequence training. We generally use state-level Minimum Bayes Risk (sMBR) [27,28] although we have also implemented Minimum Phone Error (MPE) [26] and Boosted MMI [25]. The high-level details of our lattice-based training procedure are similar to [29], but note that in that paper we describe an alternative implementation of deep neural nets (the "nnet1" setup) that exists within Kaldi; this paper is about the alternative "nnet2" setup. Some items in common with the sequence training described in that paper include the following: • We use a low, fixed learning rate (e.g. 0.00002).
• We generally train for about 4 epochs.
Some differences include the following: • We do parallel SGD on multiple machines, with periodic model averaging.
• Rather than randomizing at the utterance level, we split the lattice into as small pieces as possible given the lattice topology, and excise parts of the lattice that would not contribute nonzero derivatives; and we randomize the order of the remaining pieces.
• To ensure that all layers of the network are trained about the same amount, we modify the learning rates in order to ensure that the relative change in parameters on each "outer iteration" is the same for each layer; their geometric average is constrained to equal the user-specified fixed learning rate (e.g. 0.00002) which we mentioned above.
• In our recipe, we generate the lattices only once.
• The minibatches actually consist of several small chunks of lattice (split as described above), from many different utterances, spliced together.
Something that we should note in connection with the learning rates is that for p-norm networks, since the network output is invariant to (nonzero) scaling of the parameters of the p-norm layers 2 , and since the generalized weighted combination of Section C.2 may output arbitrarily scaled weights, it is hard to specify in advance a suitable learning rate. To solve this problem, we first scale the parameters of p-norm layers so so that the expected square of a randomly chosen matrix element is one.
For sequence training, because the frames in a minibatch are not drawn independently from the training data but consist of sequential frames from one or a few utterances, our "simple" preconditioning method is not applicable, and we only apply the online method.

C.7 Online decoding and iVector inputs
In speech recognition applications it is sometimes necessary to process data continuously as it arrives, so that there will be no latency in response. This makes it necessary that the algorithms used should not have any dependencies that are "backwards" in time.
Backwards-in-time dependencies in our conventional neural net recipes, e.g. as reported in [15], include cepstral mean normalization (CMN), in which we subtract the mean of the input features; and fMLLR adaptation, also known as constrained MLLR adaptation [30], in which we use a baseline GMM system to compute a likelihood-maximizing linear transform of the features. Although we use "online" versions of both of these things for online GMM-based decoding, it makes the system very complex and is not ideal for combination with DNNs.
In order to have a system that is easier to turn into an online algorithm, we use iVectors [31] as an additional input to the neural network, in addition to the spliced cepstral features. An iVector is a vector normally of dimension in the range of several hundred, that represents speaker characteristics in a form suitable for speaker identification, and which is extracted in a Maximum Likelihood way in conjunction with a single mixture-of-Gaussians model (their means are regressed on the iVector). The parameters of the factor analysis model that extracts the iVectors are trained without any supervision, just on a large number of audio recordings. In our case we extract iVectors of dimension 100. Once the iVector extractor is trained, we switch to "online" extraction of iVectors for both training and decoding, in which only frames preceding the current frame are taken as inputs to the iVector estimation process. At the beginning of the utterance, the iVector will be zero due to the prior term.
The actual inputs to the DNN in this setup normally consist of the iVector, plus ±7 frames of Mel frequency cepstral coefficients (MFCCs) [32], without cepstral mean normalization. Some other authors [33] use log Mel filterbank energies; the MFCC features we use here are equivalent to log Mel filterbank energies because MFCCs are a linear transform of them (we use the same number of coefficients as filterbanks, 40 for these experiments) and our input data normalization (Section C.4) is invariant to such transforms; we only use MFCCs because they are more easily compressible and our "training example" data structure (Section 3.3) compresses the input features.
In order to train models that are well matched both to per-speaker decoding, where statistics from previous utterances of the same speaker are included in the iVector estimation, and per-utterance decoding, where we make a fresh start each time, we generally train after splitting the speakers into "fake" speakers that each have no more than two utterances.
In experiments on a number of datasets, we have generally found that this method gives us about the same performance as our previous recipe where we trained a DNN on top of ±4 frames of the standard 40-dimensional features consisting of mean-normalized MFCC features processed with LDA and MLLT, and