# A Simplified Natural Gradient Learning

Authors: Michael R. Bastian; Jacob H. Gunther; Todd K. Moon
Journal: Advances in Artificial Neural Systems
 (2011)

---

## Abstract

Adaptive natural gradient learning avoids singularities in the parameter
space of multilayer perceptrons. However, it requires a larger number
of additional parameters than ordinary backpropagation in the form of
the Fisher information matrix. This paper describes a new approach to
natural gradient learning that uses a smaller Fisher information matrix.
It also uses a prior distribution on the neural network parameters and an
annealed learning rate. While this new approach is computationally simpler,
its performance is comparable to that of adaptive natural gradient
learning.


---

## Body

## 1. Introduction

Amari et al. developed the adaptive natural gradient learning (ANGL) algorithm for multilayer perceptrons [1–3]. It is similar to Newton's method in which the gradient of the error function is multiplied by the inverse of the Hessian matrix of the error function. However, Amari’s algorithm replaces the inverse of the Hessian matrix with the inverse of the Fisher information matrix. The simplified natural gradient learning (SNGL) algorithm introduced in this paper uses a new formulation of the Fisher information matrix. SNGL is based on the backpropagation algorithm [4]. In addition, the SNGL algorithm also uses regularization [5] to penalize solutions with large connection weights. This regularization also prevents the Fisher information matrix from approaching a singular matrix. The learning rate is also annealed [6] to improve the general performance of SNGL.The remainder of this section will review the natural gradient and discuss the details of Amari's ANGL algorithm in [2, 7, 8]. The next section will introduce the SNGL algorithm and the theory behind it. The third section will describe two improvements used in the SNGL algorithm. These improvements are regularization and annealed learning rates. The following section will show the experimental results and compare the performance of these two algorithms with that of ordinary backpropagation. The final section will conclude the ideas discussed in this paper.


### 1.1. Mathematical Preliminaries and Notation

The natural gradient learning algorithm is “natural” in the sense that it follows the manifold of the underlying parameter space of multilayer perceptrons. This parameter space has Riemannian geometry. Calculating the direction of steepest descent in a Riemannian manifold is more complicated than it is in a Euclidean space.The parameter space of multilayer perceptrons is an affine space in which probability distributions are the points and random variables are vectors [9, 10]. Given a probability density function p of a random variable x, the translation operation is p+x=pexp(x).


#### 1.1.1. Exponential Families

An exponential family is a set of distributions with probability density functions of the form(1)p(x∣θ1,…,θr)=c(x)exp⁡(∑i=1rθiηi(x)-A(θ1,…,θr)).
Thus, the distributions in the exponential family are parameterized by {θ1,…,θr}, and these parameters are the canonical coordinates of the distribution. The number A(θ1,…,θr) is chosen such that ∫p(x∣θ1,…,θr)dx=1. Given an open set U in the event space of the random variable X, the probability that the result of the random variable X lies within U is(2)P[X∈U]=∫Up(x∣θ1,…,θr)dx.The sufficient statistics for the exponential family are ther random variables η1,…,ηr in (1). These statistics can be used to estimate the coordinates θ1,…,θr from the values of samples of x.The log-likelihood function of a probability distribution is the logarithm of its probability density functionlogp. The score functions are the derivatives of the log-likelihood function with respect to the coordinates θ1,…,θr. The score for parameter i is(3)si=∂log⁡p(x∣θ1,…,θr)∂θi=ηi(x)-∂A∂θi.The Fisher information of two score functions is the expected value of their product. For the score functionssi and sj, the Fisher information is(4)gpij=Ep{sisj},
where Ep is the expectation operator for the distribution with probability density function p. The Fisher information matrix is a real symmetric r×r matrix Gp=(gpij) whose elements are the Fisher information for each pair of score functions.


#### 1.1.2. Direction of Steepest Descent

When finding the minimum of a functionf in a Euclidean space, a descent direction [11] is a unit norm tangent vector vp such that(5)〈vp,∇f(p)〉≜∑i=1nvi∂f∂xi|p<0.
The minimum value of (5) is -α∥∇f(p)∥2 where α is an arbitrary scalar and ∥∇f(p)∥2=∑i=1n(∂f/∂xi∣p)2 which occurs when vp=-∇f(p)/∥∇f(p)∥. Thus, in a Euclidean space, the direction of steepest descent is the unit tangent vector -∇f(p)/∥∇f(p)∥.Finding the direction of steepest descent turns out to be harder than it first appears to be. In a Euclidean geometry, it is relatively straightforward. In a Riemannian geometry, the coordinate system is curved, and a gradient descent algorithm should follow the curvature of the manifold. The direction of steepest descent should be compensated for this curvature.In a Riemannian space, a metricgp exists at every point p. If the tangent vectors vp1,…,vpn form a basis for the tangent space Tp at p, then gpij=〈vpi,vpj〉 is the Riemannian metric where 〈·,·〉 is an inner product defined on the tangent space Tp.Given a tangent vectorwp=∑i=1nαivpi and the gradient ∇f(p)=∑i=1nβivpi, their inner product will be(6)〈wp,∇f(p)〉=∑i=1n∑j=1ngpijαiβj.
If wp is the steepest descent direction, then, according to the analysis above, 〈wp,∇f(p)〉=-α∥∇f(p)∥2. This occurs when wp=-Gp-1∇f(p) where Gp=(gpij) is the matrix whose elements are equal to the inner products of the basis vectors. In other words, the coordinate system matters when doing gradient descent, and it may change at each point on the manifold. The natural gradient learning algorithm uses this “natural” descent direction at each step.


#### 1.1.3. A Probability Model for Multilayer Perceptrons

In order for this analysis to be applied to neural networks, there must be a probability density function associated with a neural network. This section describes an exponential family in which each member of the family has a probability density function of network output errors. This exponential family will then be used to determine the direction of steepest descent and the natural gradient to be used in a learning algorithm.Let the vectors{x1,…,xn} and {y1,…,yn}, respectively, be the inputs and outputs of a training set. A convenient model for training a multilayer perceptron is(7)yi=ϕ(xi)+ϵi,
where ϕ(xi) represents the output of the multilayer perceptron for training input i, and each ϵi~𝒩(0,Σ) represents the error on each training pair. A probability distribution can then be defined with the density function(8)p(yi∣xi)=1(2π)mexp⁡(-12‖yi-ϕ(xi)‖2),
where m is the dimension of the vector space to which y1,…,yn belong. In this model we assume that each error vector is independent and identically distributed so that the joint density of all the training patterns is(9)p(y1,…,yn∣x1,…,xn)=∏i=1np(yi∣xi).With the training set and the network architecture (e.g., number of hidden nodes, number of hidden layers, transfer functions) specified, the probability distributions in (8) can be points in a manifold [9]. The network parameters are the coordinates of the points. This is called a neuromanifold [10]. Random variables are the tangent vectors at the point p. If v is a random variable, then the point p+v is (by definition) the distribution with the density function pexp(v). The score functions in (3) of the network parameters form a basis for the tangent space at a point [12]. If the inner product of this tangent space is defined to be 〈v,w〉=Ep{vw}, where Ep is the expectation operator with density function p, then the Riemannian metric is the Fisher information matrix because the Fisher information between two score functions is gij=Ep{sisj} [9, 10, 12]. These definitions are illustrated in the following example.


#### 1.1.4. Natural Gradient with a Single Neuron

This is an example of a neuromanifold. Consider a neural network that consists of a single neuron that is described by a weight vectorw of dimension d and a bias b. This example will show how the Fisher information matrix is determined for a neural network. It also shows some of the pitfalls that are encountered when computing its inverse.The log-likelihood of a distributionp in a neuromanifold is(10)lp=-m2log⁡(2π)-12∑i=1n‖yi-ϕ(xi)‖2.
If the network is a single neuron where x is a d-dimensional input vector,(11)ϕ(x)=f(w⋅x+b),
and f is a sigmoid function such as tanh or the logistic function, then the score functions are(12)∂lp∂w=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)]xi,∂lp∂b=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)].
The gradient of ℓp is(13)∇lp=[∂lp∂w∂lp∂b].
The Fisher information matrix is Gp=Ep{∇ℓp(∇ℓp)T}. In this case, the straightforward calculation yields(14)Gp=∑i=1n[f′(w⋅xi+b)]2[xixiTxixiT1],
assuming that Ep{[yi-f(w·xi+b)]2}=1.One property of sigmoid functions is that their derivatives are small outside of the region near the origin. Therefore, as the norm ofw or the bias b increases, the derivative of the output tends toward zero. Since the derivative is squared in (14), the Fisher information matrix Gp could become nearly singular.The matrixGp will also be singular unless there are d linearly independent input vectors. For example, if the input vectors of the training set were three-dimensional vectors, then there must be three linearly independent vectors in order to successfully invert Gp. Care must be taken when computing the inverse of Gp. Calculating Gp becomes much more complicated with the complex network architectures of multilayer perceptrons.


### 1.2. Adaptive Natural Gradient

Multilayer perceptrons can have many parameters. If the dimension of the tangent spaceTp is m, then Gp will be an m×m matrix. Inverting the matrix has computational complexity O(m3). In the calculation of the natural gradient, the inverse of the Fisher information matrix needs to be known at each step of the algorithm. However, the Fisher information matrix itself does not need to be calculated directly for each step.Amari et al. developed an algorithm called adaptive natural gradient learning (ANGL) [2]. It uses the matrix inversion lemma [12] to perform a rank 1 update on the inverse of the Fisher information matrix at each step. The parameters are arranged and indexed as a column vector. The score function sk=∇ℓp then has the same form and is used in the updated formula(15)Gk+1-1=(1+δ)Gk-1-δGk-1skskTGk-11+skTGk-1sk.
The variable k is an integer representing the current iteration number of the learning algorithm. At each value of k, new values of the network parameters are used to compute the training error vectors. These errors are used to calculate the changes to the network parameters. The inverse of the Fisher information matrix is first initialized as G0-1=δ-1I.As mentioned in the previous paragraph, the ANGL algorithm stacks all the network parameters into one tall column vector. This is typically done when using a nonlinear least squares algorithm like Levenberg-Marquardt [11]. For networks with multiple layers, this leads to tall vectors and a large Fisher information matrix. Amari avoids inverting the matrix by using the matrix inversion lemma. However, the matrix is still usually large enough to be an issue with memory.The ANGL algorithm performs well, but it needs sufficient memory to store a large matrix and is sensitive to initial conditions and learning rates.


## 1.1. Mathematical Preliminaries and Notation

The natural gradient learning algorithm is “natural” in the sense that it follows the manifold of the underlying parameter space of multilayer perceptrons. This parameter space has Riemannian geometry. Calculating the direction of steepest descent in a Riemannian manifold is more complicated than it is in a Euclidean space.The parameter space of multilayer perceptrons is an affine space in which probability distributions are the points and random variables are vectors [9, 10]. Given a probability density function p of a random variable x, the translation operation is p+x=pexp(x).


### 1.1.1. Exponential Families

An exponential family is a set of distributions with probability density functions of the form(1)p(x∣θ1,…,θr)=c(x)exp⁡(∑i=1rθiηi(x)-A(θ1,…,θr)).
Thus, the distributions in the exponential family are parameterized by {θ1,…,θr}, and these parameters are the canonical coordinates of the distribution. The number A(θ1,…,θr) is chosen such that ∫p(x∣θ1,…,θr)dx=1. Given an open set U in the event space of the random variable X, the probability that the result of the random variable X lies within U is(2)P[X∈U]=∫Up(x∣θ1,…,θr)dx.The sufficient statistics for the exponential family are ther random variables η1,…,ηr in (1). These statistics can be used to estimate the coordinates θ1,…,θr from the values of samples of x.The log-likelihood function of a probability distribution is the logarithm of its probability density functionlogp. The score functions are the derivatives of the log-likelihood function with respect to the coordinates θ1,…,θr. The score for parameter i is(3)si=∂log⁡p(x∣θ1,…,θr)∂θi=ηi(x)-∂A∂θi.The Fisher information of two score functions is the expected value of their product. For the score functionssi and sj, the Fisher information is(4)gpij=Ep{sisj},
where Ep is the expectation operator for the distribution with probability density function p. The Fisher information matrix is a real symmetric r×r matrix Gp=(gpij) whose elements are the Fisher information for each pair of score functions.


### 1.1.2. Direction of Steepest Descent

When finding the minimum of a functionf in a Euclidean space, a descent direction [11] is a unit norm tangent vector vp such that(5)〈vp,∇f(p)〉≜∑i=1nvi∂f∂xi|p<0.
The minimum value of (5) is -α∥∇f(p)∥2 where α is an arbitrary scalar and ∥∇f(p)∥2=∑i=1n(∂f/∂xi∣p)2 which occurs when vp=-∇f(p)/∥∇f(p)∥. Thus, in a Euclidean space, the direction of steepest descent is the unit tangent vector -∇f(p)/∥∇f(p)∥.Finding the direction of steepest descent turns out to be harder than it first appears to be. In a Euclidean geometry, it is relatively straightforward. In a Riemannian geometry, the coordinate system is curved, and a gradient descent algorithm should follow the curvature of the manifold. The direction of steepest descent should be compensated for this curvature.In a Riemannian space, a metricgp exists at every point p. If the tangent vectors vp1,…,vpn form a basis for the tangent space Tp at p, then gpij=〈vpi,vpj〉 is the Riemannian metric where 〈·,·〉 is an inner product defined on the tangent space Tp.Given a tangent vectorwp=∑i=1nαivpi and the gradient ∇f(p)=∑i=1nβivpi, their inner product will be(6)〈wp,∇f(p)〉=∑i=1n∑j=1ngpijαiβj.
If wp is the steepest descent direction, then, according to the analysis above, 〈wp,∇f(p)〉=-α∥∇f(p)∥2. This occurs when wp=-Gp-1∇f(p) where Gp=(gpij) is the matrix whose elements are equal to the inner products of the basis vectors. In other words, the coordinate system matters when doing gradient descent, and it may change at each point on the manifold. The natural gradient learning algorithm uses this “natural” descent direction at each step.


### 1.1.3. A Probability Model for Multilayer Perceptrons

In order for this analysis to be applied to neural networks, there must be a probability density function associated with a neural network. This section describes an exponential family in which each member of the family has a probability density function of network output errors. This exponential family will then be used to determine the direction of steepest descent and the natural gradient to be used in a learning algorithm.Let the vectors{x1,…,xn} and {y1,…,yn}, respectively, be the inputs and outputs of a training set. A convenient model for training a multilayer perceptron is(7)yi=ϕ(xi)+ϵi,
where ϕ(xi) represents the output of the multilayer perceptron for training input i, and each ϵi~𝒩(0,Σ) represents the error on each training pair. A probability distribution can then be defined with the density function(8)p(yi∣xi)=1(2π)mexp⁡(-12‖yi-ϕ(xi)‖2),
where m is the dimension of the vector space to which y1,…,yn belong. In this model we assume that each error vector is independent and identically distributed so that the joint density of all the training patterns is(9)p(y1,…,yn∣x1,…,xn)=∏i=1np(yi∣xi).With the training set and the network architecture (e.g., number of hidden nodes, number of hidden layers, transfer functions) specified, the probability distributions in (8) can be points in a manifold [9]. The network parameters are the coordinates of the points. This is called a neuromanifold [10]. Random variables are the tangent vectors at the point p. If v is a random variable, then the point p+v is (by definition) the distribution with the density function pexp(v). The score functions in (3) of the network parameters form a basis for the tangent space at a point [12]. If the inner product of this tangent space is defined to be 〈v,w〉=Ep{vw}, where Ep is the expectation operator with density function p, then the Riemannian metric is the Fisher information matrix because the Fisher information between two score functions is gij=Ep{sisj} [9, 10, 12]. These definitions are illustrated in the following example.


### 1.1.4. Natural Gradient with a Single Neuron

This is an example of a neuromanifold. Consider a neural network that consists of a single neuron that is described by a weight vectorw of dimension d and a bias b. This example will show how the Fisher information matrix is determined for a neural network. It also shows some of the pitfalls that are encountered when computing its inverse.The log-likelihood of a distributionp in a neuromanifold is(10)lp=-m2log⁡(2π)-12∑i=1n‖yi-ϕ(xi)‖2.
If the network is a single neuron where x is a d-dimensional input vector,(11)ϕ(x)=f(w⋅x+b),
and f is a sigmoid function such as tanh or the logistic function, then the score functions are(12)∂lp∂w=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)]xi,∂lp∂b=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)].
The gradient of ℓp is(13)∇lp=[∂lp∂w∂lp∂b].
The Fisher information matrix is Gp=Ep{∇ℓp(∇ℓp)T}. In this case, the straightforward calculation yields(14)Gp=∑i=1n[f′(w⋅xi+b)]2[xixiTxixiT1],
assuming that Ep{[yi-f(w·xi+b)]2}=1.One property of sigmoid functions is that their derivatives are small outside of the region near the origin. Therefore, as the norm ofw or the bias b increases, the derivative of the output tends toward zero. Since the derivative is squared in (14), the Fisher information matrix Gp could become nearly singular.The matrixGp will also be singular unless there are d linearly independent input vectors. For example, if the input vectors of the training set were three-dimensional vectors, then there must be three linearly independent vectors in order to successfully invert Gp. Care must be taken when computing the inverse of Gp. Calculating Gp becomes much more complicated with the complex network architectures of multilayer perceptrons.


## 1.1.1. Exponential Families

An exponential family is a set of distributions with probability density functions of the form(1)p(x∣θ1,…,θr)=c(x)exp⁡(∑i=1rθiηi(x)-A(θ1,…,θr)).
Thus, the distributions in the exponential family are parameterized by {θ1,…,θr}, and these parameters are the canonical coordinates of the distribution. The number A(θ1,…,θr) is chosen such that ∫p(x∣θ1,…,θr)dx=1. Given an open set U in the event space of the random variable X, the probability that the result of the random variable X lies within U is(2)P[X∈U]=∫Up(x∣θ1,…,θr)dx.The sufficient statistics for the exponential family are ther random variables η1,…,ηr in (1). These statistics can be used to estimate the coordinates θ1,…,θr from the values of samples of x.The log-likelihood function of a probability distribution is the logarithm of its probability density functionlogp. The score functions are the derivatives of the log-likelihood function with respect to the coordinates θ1,…,θr. The score for parameter i is(3)si=∂log⁡p(x∣θ1,…,θr)∂θi=ηi(x)-∂A∂θi.The Fisher information of two score functions is the expected value of their product. For the score functionssi and sj, the Fisher information is(4)gpij=Ep{sisj},
where Ep is the expectation operator for the distribution with probability density function p. The Fisher information matrix is a real symmetric r×r matrix Gp=(gpij) whose elements are the Fisher information for each pair of score functions.


## 1.1.2. Direction of Steepest Descent

When finding the minimum of a functionf in a Euclidean space, a descent direction [11] is a unit norm tangent vector vp such that(5)〈vp,∇f(p)〉≜∑i=1nvi∂f∂xi|p<0.
The minimum value of (5) is -α∥∇f(p)∥2 where α is an arbitrary scalar and ∥∇f(p)∥2=∑i=1n(∂f/∂xi∣p)2 which occurs when vp=-∇f(p)/∥∇f(p)∥. Thus, in a Euclidean space, the direction of steepest descent is the unit tangent vector -∇f(p)/∥∇f(p)∥.Finding the direction of steepest descent turns out to be harder than it first appears to be. In a Euclidean geometry, it is relatively straightforward. In a Riemannian geometry, the coordinate system is curved, and a gradient descent algorithm should follow the curvature of the manifold. The direction of steepest descent should be compensated for this curvature.In a Riemannian space, a metricgp exists at every point p. If the tangent vectors vp1,…,vpn form a basis for the tangent space Tp at p, then gpij=〈vpi,vpj〉 is the Riemannian metric where 〈·,·〉 is an inner product defined on the tangent space Tp.Given a tangent vectorwp=∑i=1nαivpi and the gradient ∇f(p)=∑i=1nβivpi, their inner product will be(6)〈wp,∇f(p)〉=∑i=1n∑j=1ngpijαiβj.
If wp is the steepest descent direction, then, according to the analysis above, 〈wp,∇f(p)〉=-α∥∇f(p)∥2. This occurs when wp=-Gp-1∇f(p) where Gp=(gpij) is the matrix whose elements are equal to the inner products of the basis vectors. In other words, the coordinate system matters when doing gradient descent, and it may change at each point on the manifold. The natural gradient learning algorithm uses this “natural” descent direction at each step.


## 1.1.3. A Probability Model for Multilayer Perceptrons

In order for this analysis to be applied to neural networks, there must be a probability density function associated with a neural network. This section describes an exponential family in which each member of the family has a probability density function of network output errors. This exponential family will then be used to determine the direction of steepest descent and the natural gradient to be used in a learning algorithm.Let the vectors{x1,…,xn} and {y1,…,yn}, respectively, be the inputs and outputs of a training set. A convenient model for training a multilayer perceptron is(7)yi=ϕ(xi)+ϵi,
where ϕ(xi) represents the output of the multilayer perceptron for training input i, and each ϵi~𝒩(0,Σ) represents the error on each training pair. A probability distribution can then be defined with the density function(8)p(yi∣xi)=1(2π)mexp⁡(-12‖yi-ϕ(xi)‖2),
where m is the dimension of the vector space to which y1,…,yn belong. In this model we assume that each error vector is independent and identically distributed so that the joint density of all the training patterns is(9)p(y1,…,yn∣x1,…,xn)=∏i=1np(yi∣xi).With the training set and the network architecture (e.g., number of hidden nodes, number of hidden layers, transfer functions) specified, the probability distributions in (8) can be points in a manifold [9]. The network parameters are the coordinates of the points. This is called a neuromanifold [10]. Random variables are the tangent vectors at the point p. If v is a random variable, then the point p+v is (by definition) the distribution with the density function pexp(v). The score functions in (3) of the network parameters form a basis for the tangent space at a point [12]. If the inner product of this tangent space is defined to be 〈v,w〉=Ep{vw}, where Ep is the expectation operator with density function p, then the Riemannian metric is the Fisher information matrix because the Fisher information between two score functions is gij=Ep{sisj} [9, 10, 12]. These definitions are illustrated in the following example.


## 1.1.4. Natural Gradient with a Single Neuron

This is an example of a neuromanifold. Consider a neural network that consists of a single neuron that is described by a weight vectorw of dimension d and a bias b. This example will show how the Fisher information matrix is determined for a neural network. It also shows some of the pitfalls that are encountered when computing its inverse.The log-likelihood of a distributionp in a neuromanifold is(10)lp=-m2log⁡(2π)-12∑i=1n‖yi-ϕ(xi)‖2.
If the network is a single neuron where x is a d-dimensional input vector,(11)ϕ(x)=f(w⋅x+b),
and f is a sigmoid function such as tanh or the logistic function, then the score functions are(12)∂lp∂w=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)]xi,∂lp∂b=∑i=1nf′(w⋅xi+b)[yi-f(w⋅xi+b)].
The gradient of ℓp is(13)∇lp=[∂lp∂w∂lp∂b].
The Fisher information matrix is Gp=Ep{∇ℓp(∇ℓp)T}. In this case, the straightforward calculation yields(14)Gp=∑i=1n[f′(w⋅xi+b)]2[xixiTxixiT1],
assuming that Ep{[yi-f(w·xi+b)]2}=1.One property of sigmoid functions is that their derivatives are small outside of the region near the origin. Therefore, as the norm ofw or the bias b increases, the derivative of the output tends toward zero. Since the derivative is squared in (14), the Fisher information matrix Gp could become nearly singular.The matrixGp will also be singular unless there are d linearly independent input vectors. For example, if the input vectors of the training set were three-dimensional vectors, then there must be three linearly independent vectors in order to successfully invert Gp. Care must be taken when computing the inverse of Gp. Calculating Gp becomes much more complicated with the complex network architectures of multilayer perceptrons.


## 1.2. Adaptive Natural Gradient

Multilayer perceptrons can have many parameters. If the dimension of the tangent spaceTp is m, then Gp will be an m×m matrix. Inverting the matrix has computational complexity O(m3). In the calculation of the natural gradient, the inverse of the Fisher information matrix needs to be known at each step of the algorithm. However, the Fisher information matrix itself does not need to be calculated directly for each step.Amari et al. developed an algorithm called adaptive natural gradient learning (ANGL) [2]. It uses the matrix inversion lemma [12] to perform a rank 1 update on the inverse of the Fisher information matrix at each step. The parameters are arranged and indexed as a column vector. The score function sk=∇ℓp then has the same form and is used in the updated formula(15)Gk+1-1=(1+δ)Gk-1-δGk-1skskTGk-11+skTGk-1sk.
The variable k is an integer representing the current iteration number of the learning algorithm. At each value of k, new values of the network parameters are used to compute the training error vectors. These errors are used to calculate the changes to the network parameters. The inverse of the Fisher information matrix is first initialized as G0-1=δ-1I.As mentioned in the previous paragraph, the ANGL algorithm stacks all the network parameters into one tall column vector. This is typically done when using a nonlinear least squares algorithm like Levenberg-Marquardt [11]. For networks with multiple layers, this leads to tall vectors and a large Fisher information matrix. Amari avoids inverting the matrix by using the matrix inversion lemma. However, the matrix is still usually large enough to be an issue with memory.The ANGL algorithm performs well, but it needs sufficient memory to store a large matrix and is sensitive to initial conditions and learning rates.


## 2. A Simplified Natural Gradient Learning  Algorithm

This section describes the learning rule used in the SNGL algorithm. First, the directional derivative for a single-layer perceptron is determined. Then the Fisher information matrix for this network is described in detail. These two results are then generalized to find the directional derivative and the Fisher information matrix for a multilayer perceptron. Finally, the section will conclude with a discussion of the computational complexity of SNGL compared to that of ANGL.The derivation and analysis of the simplified natural gradient learning algorithm begin with the observation that the neural network is parameterized by(A,b) in the affine transformation Ax+b. There is one affine transformation for each layer of the network. The formula in (11) was for a single neuron that has only one output. In general, the affine transformations in each layer of a multilayer perceptron can have more than one output. The objective of the algorithm is to minimize the sum of squares of the error vector norms. These norms are calculated with the standard L2 norm of a Cartesian coordinate system. However, that does not mean that the matrix parameters need to be column-ordered into tall vectors. While the homomorphism as done in [2] preserves the vector space operations such as addition and multiplication by a scalar, it does not preserve multiplying the matrix by a vector. Therefore, it makes sense that the score function of a matrix parameter should itself be a matrix. The Frobenius norm of a m×n Matrix is tr(ATA)=∑i=1m∑j=1naij2. The inner product associated with this norm is 〈A,B〉=tr(ATB). The following analysis shows how the derivative of the norm of an error vector can be equivalent to an inner product between matrices.


### 2.1. The Directional Derivative

Assume that the network is a single-layer perceptron with the vector-valued sigmoidal transfer functionF(·). Thus, the log-likelihood function ℓ with n training examples that depends on the parameters A and b is(16)lϕ=-12∑i=1n‖yi-F(Axi+b)‖2.
Furthermore, assume that if the affine transformation ψ(x)=Cx+d is added to ϕ(x)=Ax+b as (ϕ+ψ)(x)=(A+C)x+(b+d), then the minimum point, according to (16), is reached. Thus, C and d represent a descent direction.Letγ(t,x)=ϕ(x)+tψ(x)=(A+tC)x+b+td be a curve through the manifold for 0≤t≤1. The affine transformation γ(0,x) is the initial point, and γ(1,x) is the final point. The derivative of this curve with respect to t is γ̇(t,x)=Cx+d and is constant with respect to t. The log-likelihood along this curve is(17)lγ(t)=-12∑i=1n‖yi-F(γ(t,xi))‖2.
The directional derivative of the log-likelihood with respect to t is(18)dlγdt=∑i=1n〈yi-F(γ(t,xi)),F′(γ(t,xi))γ̇(t,xi)〉=∑i=1n〈[F′(γ(t,xi))]T[yi-F(γ(t,xi))],γ̇(t,xi)〉,
where F′(·) is the Jacobian matrix of F(·). Here, we have isolated γ̇(t,x) by using the adjoint of F′(·) in the second argument of the inner product. Given that γ(0,x)=Ax+b and γ̇(0,x)=Cx+d, the directional derivative of the log-likelihood function at t=0 is(19)dlγdt|t=0=∑i=1n〈[F′(Axi+b)]T[yi-F(Axi+b)],Cxi+d〉.Letzi=F′(Axi+b)T[yi-F(Axi+b)]. Then, by using the linearity property of the inner product, (19) can be rewritten as(20)dlγdt|t=0=∑i=1n〈zi,Cxi〉+∑i=1n〈zi,d〉.
The matrix C can be isolated because the matrix inner product 〈A,B〉=tr(ATB) is being used on the right-hand side. Thus, 〈zi,Cxi〉=〈zixiT,C〉.Sinceℓγ=ℓϕ∘γ(t,x), then, by the chain rule, dℓγ/dt=〈∇ℓϕ,γ̇(t,x)〉. Thus, 〈∇ℓϕ,γ̇(t,x)〉=∑i=1n〈∂ℓϕ/∂A,C〉+〈∂ℓϕ/∂b,d〉=∑i=1n〈zixiT,C〉+〈zi,d〉 and the components of the score function are(21)∂lϕ∂A=∑i=1n[F′(Axi+b)]T[yi-F(Axi+b)]xiT,∂lϕ∂b=∑i=1n[F′(Axi+b)]T[yi-F(Axi+b)].
The components in (21) would point in the direction of steepest descent if the space of affine transformations was a Euclidean space. To find the direction of steepest descent, then (21) must be corrected for the curvature of the parameter space using the Riemannian metric G.


### 2.2. The Fisher Information Matrix

The Fisher information matrix is the Riemannian metric in this manifold [1]. For the affine transformation ϕ(x)=Ax+b, the Fisher information matrix is(22)Gϕ=Ep{∇lϕ∇lϕT}.
The score function with the components in (21) can be arranged in an augmented matrix form as(23)∇lϕ=[∂lϕ∂A∂lϕ∂b].Thus, in terms of these components, the Fisher information matrix is(24)Gϕ=Ep{[∂lϕ∂A∂lϕ∂b][(∂lϕ∂A)T(∂lϕ∂b)T]}=Ep{∂lϕ∂A(∂lϕ∂A)T+∂lϕ∂b(∂lϕ∂b)T}.
The partial derivative ∂ℓϕ/∂A in (21) is a sum of n rank 1 matrices. Therefore, the outer product of ∂ℓϕ/∂A with itself is the sum(25)Rij=[F′(Axi+b)]TϵixiTxjϵjTF′(Axj+b),
where ϵi=yi-F(Axi+b). When i≠j then Ep{Rij}=0 because Ep{ϵiϵj}=Ep{ϵi}Ep{ϵj}=0·0=0. When i=j then(26)Ep{Rii}=‖xi‖2[F′(Axi+b)]T𝚺F′(Axi+b).
However, since Σ is not known, it must be estimated. The unbiased estimator is(27)𝚺̂=1n-1∑i=1nϵiϵiT.
The Fisher information of the weight matrix A is estimated as(28)ĜA=1n-1∑i=1n‖xi‖2[F′(Axi+b)]TϵiϵiTF′(Axi+b),
and the Fisher information of the bias vector b is estimated as(29)Ĝb=1n-1∑i=1n[F′(Axi+b)]TϵiϵiTF′(Axi+b).
The Fisher information matrix can be estimated as the sum of these two matrices. Combining like terms gives(30)Ĝϕ=1n-1∑i=1n(‖xi‖2+1)[F′(ϕ(xi))]TϵiϵiTF′(ϕ(xi)).
The direction of steepest descent in (19) is now represented as(31)C=Ĝϕ-1∂lϕ∂A,d=Ĝϕ-1∂lϕ∂b.These equations are for a perceptron with only one layer. Please note that by using the matrix forms, the Fisher information matrix is only ad×d matrix. If we had stacked the weight matrix and bias vector into a single column vector, then the Fisher information matrix would be d(d+1)×d(d+1) matrix.


### 2.3. Extension to Multilayer Perceptrons

A multilayer perceptron hasr layers each with its own nonlinear sigmoidal function Fk(·) and affine transformation ϕk(x)=Akx+bk. Each layer k feeds into the one above it with index k+1. Thus, the output of layer k is Fk∘ϕk∘Fk-1∘ϕk-1∘⋯∘F1∘ϕ1(x). The input layer is defined to be the identity transformation ϕ0(x)=x and F0(x)=x. For convenience, let the vector zik be the output of layer k when given input xi. The output zi can be written recursively,(32)zik=Fk(ϕk(zik-1))
with zi0=xi.Two vector spacesV and W can be joined together by means of a direct sum V⊕W in which stacking a vector v∈V on top of a vector w∈W results in the vector(33)[vw]∈V⊕W.
The analog to the direct sum with matrices is to create a block matrix. Given two matrices A and B, their direct sum is(34)A⊕B=[A00B].
The direct sum of the score functions of the layers in a multilayer perceptron is used as the analog for stacking the layers. The Fisher information matrix is then a block-diagonal matrix where each square matrix along the main diagonal is the Fisher information matrix for a specific layer.The score functions for layerk are(35)∂lϕk∂Ak=∑i=1nBik(yi-zir)(zik-1)T,∂lϕk∂bk=∑i=1nBik(yi-zir),
where Bik is the backpropagation matrix that is determined with the recurrence relation(36)Bir=Fr′(ϕr(zir-1))T,Bik=Fk′(ϕk(zik-1))TAk+1TBik+1,
where r is the number of layers in the network. Thus, for the output layer, the layer index is r and the backpropagation matrix is the Jacobian matrix of the output layer's transfer function.Letϵi=yi-zir be the error for pattern i. The estimate of the Fisher information matrix for layer k is(37)Ĝϕk=1n-1∑i=1n(‖zik-1‖2+1)BikϵiϵiT(Bik)T.
Thus, the updates for the parameters in layer k are(38)Ck=Ĝϕk-1∂lϕk∂Ak,dk=Ĝϕk-1∂lϕk∂bk.


### 2.4. Computational Complexity

The new formulation of the Fisher information matrix is significantly smaller than what was used by Amari et al. [2]. Suppose that a two-layer perceptron has m inputs, n outputs, and p hidden nodes. The weight matrices have dimensions (m+1)×p and (p+1)×n. The addition of 1 to both m and p is due to the extra bias vectors of dimensions p and n. If the parameters are column ordered into a single vector, then that vector has dimension (m+1)p+(p+1)n. In the ANGL algorithm's Fisher information matrix, there are ((m+1)p+(p+1)n)2 elements. The SNGL algorithm's block-diagonal Fisher information matrix has two blocks on the main diagonal that together have p2+n2 elements. Clearly,(39)p2+n2<(m+1)2p2+2(m+1)p(p+1)n+(p+1)2n2,
if m, n, and p are positive integers.It is also less computationally intensive to invert the block matrix used in the SNGL algorithm since the matrix inversion algorithms would have complexityO(p3+n3). The complexity to update Amari's recursion in (15) is O((m+1)2p2+2(m+1)p(p+1)n+(p+1)2n2). This takes into account that ANGL exploits the matrix inversion lemma. Table 1 shows the computational complexity in both space and time of the ANGL algorithm and the SNGL algorithm for three problems described by Park et al. [3]. The last two columns labelled “Ratio” in Table 1 are the ratios of the space and time complexities of the SNGL algorithm compared to that of the ANGL written as percentages.Table 1
Complexity comparison between the SNGL and ANGL algorithms on other problems.
AlgorithmSNGLANGLRatioComplexity dimensionSpaceTimeBothSpaceTimeExclusive OR59816%11%Fisher's Iris data25911,2252%7%Mackey-Glass time series1011,0013,7213%27%


### 2.5. SNGL Improvements

There are two more elements of the simplified natural gradient learning algorithm. The first is the regularization of the gradient descent algorithm by adding a prior distribution to the probability density function of the network errors [5]. The second is annealing the learning rate of the algorithm [6]. Neither has any significant impact on the computational complexity. The complexity of the regularization is O(mn+np) for a two-layer network. The learning rate is recalculated once per learning epoch.


#### 2.5.1. Regularization of Network Parameters

The goal of regularization is to penalize larger connection weights. This helps to discourage overtraining. The penalty of the weights is governed by the hyperparameterα. For a neuron with parameters w and b, the probability density function of these parameters in the prior distribution is(40)pα(w,b)=(α2π)q+1exp⁡(-α2(‖w‖2+b2)),
where q is the dimension of the vector space to which w belongs. Thus, the weight vector is distributed normally as w~𝒩(0,(1/α)I) and the bias as b~𝒩(0,1/α). The log-likelihood function is(41)lα(w,b)=q+12[log⁡α-log⁡(2π)]-α2(‖w‖2+b2).
The score functions of w and b with respect to ℓα are, respectively(42)∂lα∂w=-αw,∂lα∂b=-αb.
The probability density function of this prior distribution can be added to the probability density function of the network error distribution. Then the updates for the parameters in a multilayer perceptron are(43)∂lϕk,α∂Ak=∑i=1nBik(yi-zir)(zik-1)T-αAk,∂lϕk,α∂bk=∑i=1nBik(yi-zir)-αbk.When the algorithm is overtraining, the parameters are still increasing. Subtractingα times the current parameter from the update effectively counteracts the update when the norm of that parameter is large enough. The algorithm will eventually reach a steady state where the parameters are longer increasing.SinceEpα{w}=0, Epα{b}=0, Epα{wwT}=(1/α)I, and Epα{b2}=1/α, the Fisher information matrix for the regularized algorithm is(44)Ĝϕk,α=1n-1∑i=1n(‖zik-1‖2+1)BikϵiϵiT(Bik)T+αI.AddingαI to the Fisher information matrix helps the algorithm avoid inverting an almost singular matrix. Thus, if all the errors are zero, then the Fisher information matrix will be αI. The effect will be multiplying all the weight update elements by 1/α. This is the largest value the inverse will reach during the SNGL algorithm's execution.In [13] Park adds regularization to the ANGL algorithm with an Akaike Information Criterion.


#### 2.5.2. Annealed Learning Rates

Annealing the learning rate [6] is a technique in which a learning algorithm will slowly decay its learning rate after a specific epoch during the execution of the algorithm. Let the learning rate of the SNGL algorithm be ηt, where t is the current learning epoch. It is defined as(45)ηt=min⁡⁡(η0,τt),
where η0 is the initial learning rate and τ is a factor that determines which learning epoch begins the annealing process. The epoch in which the decay begins is t0=⌈τ/η0⌉.


## 2.1. The Directional Derivative

Assume that the network is a single-layer perceptron with the vector-valued sigmoidal transfer functionF(·). Thus, the log-likelihood function ℓ with n training examples that depends on the parameters A and b is(16)lϕ=-12∑i=1n‖yi-F(Axi+b)‖2.
Furthermore, assume that if the affine transformation ψ(x)=Cx+d is added to ϕ(x)=Ax+b as (ϕ+ψ)(x)=(A+C)x+(b+d), then the minimum point, according to (16), is reached. Thus, C and d represent a descent direction.Letγ(t,x)=ϕ(x)+tψ(x)=(A+tC)x+b+td be a curve through the manifold for 0≤t≤1. The affine transformation γ(0,x) is the initial point, and γ(1,x) is the final point. The derivative of this curve with respect to t is γ̇(t,x)=Cx+d and is constant with respect to t. The log-likelihood along this curve is(17)lγ(t)=-12∑i=1n‖yi-F(γ(t,xi))‖2.
The directional derivative of the log-likelihood with respect to t is(18)dlγdt=∑i=1n〈yi-F(γ(t,xi)),F′(γ(t,xi))γ̇(t,xi)〉=∑i=1n〈[F′(γ(t,xi))]T[yi-F(γ(t,xi))],γ̇(t,xi)〉,
where F′(·) is the Jacobian matrix of F(·). Here, we have isolated γ̇(t,x) by using the adjoint of F′(·) in the second argument of the inner product. Given that γ(0,x)=Ax+b and γ̇(0,x)=Cx+d, the directional derivative of the log-likelihood function at t=0 is(19)dlγdt|t=0=∑i=1n〈[F′(Axi+b)]T[yi-F(Axi+b)],Cxi+d〉.Letzi=F′(Axi+b)T[yi-F(Axi+b)]. Then, by using the linearity property of the inner product, (19) can be rewritten as(20)dlγdt|t=0=∑i=1n〈zi,Cxi〉+∑i=1n〈zi,d〉.
The matrix C can be isolated because the matrix inner product 〈A,B〉=tr(ATB) is being used on the right-hand side. Thus, 〈zi,Cxi〉=〈zixiT,C〉.Sinceℓγ=ℓϕ∘γ(t,x), then, by the chain rule, dℓγ/dt=〈∇ℓϕ,γ̇(t,x)〉. Thus, 〈∇ℓϕ,γ̇(t,x)〉=∑i=1n〈∂ℓϕ/∂A,C〉+〈∂ℓϕ/∂b,d〉=∑i=1n〈zixiT,C〉+〈zi,d〉 and the components of the score function are(21)∂lϕ∂A=∑i=1n[F′(Axi+b)]T[yi-F(Axi+b)]xiT,∂lϕ∂b=∑i=1n[F′(Axi+b)]T[yi-F(Axi+b)].
The components in (21) would point in the direction of steepest descent if the space of affine transformations was a Euclidean space. To find the direction of steepest descent, then (21) must be corrected for the curvature of the parameter space using the Riemannian metric G.


## 2.2. The Fisher Information Matrix

The Fisher information matrix is the Riemannian metric in this manifold [1]. For the affine transformation ϕ(x)=Ax+b, the Fisher information matrix is(22)Gϕ=Ep{∇lϕ∇lϕT}.
The score function with the components in (21) can be arranged in an augmented matrix form as(23)∇lϕ=[∂lϕ∂A∂lϕ∂b].Thus, in terms of these components, the Fisher information matrix is(24)Gϕ=Ep{[∂lϕ∂A∂lϕ∂b][(∂lϕ∂A)T(∂lϕ∂b)T]}=Ep{∂lϕ∂A(∂lϕ∂A)T+∂lϕ∂b(∂lϕ∂b)T}.
The partial derivative ∂ℓϕ/∂A in (21) is a sum of n rank 1 matrices. Therefore, the outer product of ∂ℓϕ/∂A with itself is the sum(25)Rij=[F′(Axi+b)]TϵixiTxjϵjTF′(Axj+b),
where ϵi=yi-F(Axi+b). When i≠j then Ep{Rij}=0 because Ep{ϵiϵj}=Ep{ϵi}Ep{ϵj}=0·0=0. When i=j then(26)Ep{Rii}=‖xi‖2[F′(Axi+b)]T𝚺F′(Axi+b).
However, since Σ is not known, it must be estimated. The unbiased estimator is(27)𝚺̂=1n-1∑i=1nϵiϵiT.
The Fisher information of the weight matrix A is estimated as(28)ĜA=1n-1∑i=1n‖xi‖2[F′(Axi+b)]TϵiϵiTF′(Axi+b),
and the Fisher information of the bias vector b is estimated as(29)Ĝb=1n-1∑i=1n[F′(Axi+b)]TϵiϵiTF′(Axi+b).
The Fisher information matrix can be estimated as the sum of these two matrices. Combining like terms gives(30)Ĝϕ=1n-1∑i=1n(‖xi‖2+1)[F′(ϕ(xi))]TϵiϵiTF′(ϕ(xi)).
The direction of steepest descent in (19) is now represented as(31)C=Ĝϕ-1∂lϕ∂A,d=Ĝϕ-1∂lϕ∂b.These equations are for a perceptron with only one layer. Please note that by using the matrix forms, the Fisher information matrix is only ad×d matrix. If we had stacked the weight matrix and bias vector into a single column vector, then the Fisher information matrix would be d(d+1)×d(d+1) matrix.


## 2.3. Extension to Multilayer Perceptrons

A multilayer perceptron hasr layers each with its own nonlinear sigmoidal function Fk(·) and affine transformation ϕk(x)=Akx+bk. Each layer k feeds into the one above it with index k+1. Thus, the output of layer k is Fk∘ϕk∘Fk-1∘ϕk-1∘⋯∘F1∘ϕ1(x). The input layer is defined to be the identity transformation ϕ0(x)=x and F0(x)=x. For convenience, let the vector zik be the output of layer k when given input xi. The output zi can be written recursively,(32)zik=Fk(ϕk(zik-1))
with zi0=xi.Two vector spacesV and W can be joined together by means of a direct sum V⊕W in which stacking a vector v∈V on top of a vector w∈W results in the vector(33)[vw]∈V⊕W.
The analog to the direct sum with matrices is to create a block matrix. Given two matrices A and B, their direct sum is(34)A⊕B=[A00B].
The direct sum of the score functions of the layers in a multilayer perceptron is used as the analog for stacking the layers. The Fisher information matrix is then a block-diagonal matrix where each square matrix along the main diagonal is the Fisher information matrix for a specific layer.The score functions for layerk are(35)∂lϕk∂Ak=∑i=1nBik(yi-zir)(zik-1)T,∂lϕk∂bk=∑i=1nBik(yi-zir),
where Bik is the backpropagation matrix that is determined with the recurrence relation(36)Bir=Fr′(ϕr(zir-1))T,Bik=Fk′(ϕk(zik-1))TAk+1TBik+1,
where r is the number of layers in the network. Thus, for the output layer, the layer index is r and the backpropagation matrix is the Jacobian matrix of the output layer's transfer function.Letϵi=yi-zir be the error for pattern i. The estimate of the Fisher information matrix for layer k is(37)Ĝϕk=1n-1∑i=1n(‖zik-1‖2+1)BikϵiϵiT(Bik)T.
Thus, the updates for the parameters in layer k are(38)Ck=Ĝϕk-1∂lϕk∂Ak,dk=Ĝϕk-1∂lϕk∂bk.


## 2.4. Computational Complexity

The new formulation of the Fisher information matrix is significantly smaller than what was used by Amari et al. [2]. Suppose that a two-layer perceptron has m inputs, n outputs, and p hidden nodes. The weight matrices have dimensions (m+1)×p and (p+1)×n. The addition of 1 to both m and p is due to the extra bias vectors of dimensions p and n. If the parameters are column ordered into a single vector, then that vector has dimension (m+1)p+(p+1)n. In the ANGL algorithm's Fisher information matrix, there are ((m+1)p+(p+1)n)2 elements. The SNGL algorithm's block-diagonal Fisher information matrix has two blocks on the main diagonal that together have p2+n2 elements. Clearly,(39)p2+n2<(m+1)2p2+2(m+1)p(p+1)n+(p+1)2n2,
if m, n, and p are positive integers.It is also less computationally intensive to invert the block matrix used in the SNGL algorithm since the matrix inversion algorithms would have complexityO(p3+n3). The complexity to update Amari's recursion in (15) is O((m+1)2p2+2(m+1)p(p+1)n+(p+1)2n2). This takes into account that ANGL exploits the matrix inversion lemma. Table 1 shows the computational complexity in both space and time of the ANGL algorithm and the SNGL algorithm for three problems described by Park et al. [3]. The last two columns labelled “Ratio” in Table 1 are the ratios of the space and time complexities of the SNGL algorithm compared to that of the ANGL written as percentages.Table 1
Complexity comparison between the SNGL and ANGL algorithms on other problems.
AlgorithmSNGLANGLRatioComplexity dimensionSpaceTimeBothSpaceTimeExclusive OR59816%11%Fisher's Iris data25911,2252%7%Mackey-Glass time series1011,0013,7213%27%


## 2.5. SNGL Improvements

There are two more elements of the simplified natural gradient learning algorithm. The first is the regularization of the gradient descent algorithm by adding a prior distribution to the probability density function of the network errors [5]. The second is annealing the learning rate of the algorithm [6]. Neither has any significant impact on the computational complexity. The complexity of the regularization is O(mn+np) for a two-layer network. The learning rate is recalculated once per learning epoch.


### 2.5.1. Regularization of Network Parameters

The goal of regularization is to penalize larger connection weights. This helps to discourage overtraining. The penalty of the weights is governed by the hyperparameterα. For a neuron with parameters w and b, the probability density function of these parameters in the prior distribution is(40)pα(w,b)=(α2π)q+1exp⁡(-α2(‖w‖2+b2)),
where q is the dimension of the vector space to which w belongs. Thus, the weight vector is distributed normally as w~𝒩(0,(1/α)I) and the bias as b~𝒩(0,1/α). The log-likelihood function is(41)lα(w,b)=q+12[log⁡α-log⁡(2π)]-α2(‖w‖2+b2).
The score functions of w and b with respect to ℓα are, respectively(42)∂lα∂w=-αw,∂lα∂b=-αb.
The probability density function of this prior distribution can be added to the probability density function of the network error distribution. Then the updates for the parameters in a multilayer perceptron are(43)∂lϕk,α∂Ak=∑i=1nBik(yi-zir)(zik-1)T-αAk,∂lϕk,α∂bk=∑i=1nBik(yi-zir)-αbk.When the algorithm is overtraining, the parameters are still increasing. Subtractingα times the current parameter from the update effectively counteracts the update when the norm of that parameter is large enough. The algorithm will eventually reach a steady state where the parameters are longer increasing.SinceEpα{w}=0, Epα{b}=0, Epα{wwT}=(1/α)I, and Epα{b2}=1/α, the Fisher information matrix for the regularized algorithm is(44)Ĝϕk,α=1n-1∑i=1n(‖zik-1‖2+1)BikϵiϵiT(Bik)T+αI.AddingαI to the Fisher information matrix helps the algorithm avoid inverting an almost singular matrix. Thus, if all the errors are zero, then the Fisher information matrix will be αI. The effect will be multiplying all the weight update elements by 1/α. This is the largest value the inverse will reach during the SNGL algorithm's execution.In [13] Park adds regularization to the ANGL algorithm with an Akaike Information Criterion.


### 2.5.2. Annealed Learning Rates

Annealing the learning rate [6] is a technique in which a learning algorithm will slowly decay its learning rate after a specific epoch during the execution of the algorithm. Let the learning rate of the SNGL algorithm be ηt, where t is the current learning epoch. It is defined as(45)ηt=min⁡⁡(η0,τt),
where η0 is the initial learning rate and τ is a factor that determines which learning epoch begins the annealing process. The epoch in which the decay begins is t0=⌈τ/η0⌉.


## 2.5.1. Regularization of Network Parameters

The goal of regularization is to penalize larger connection weights. This helps to discourage overtraining. The penalty of the weights is governed by the hyperparameterα. For a neuron with parameters w and b, the probability density function of these parameters in the prior distribution is(40)pα(w,b)=(α2π)q+1exp⁡(-α2(‖w‖2+b2)),
where q is the dimension of the vector space to which w belongs. Thus, the weight vector is distributed normally as w~𝒩(0,(1/α)I) and the bias as b~𝒩(0,1/α). The log-likelihood function is(41)lα(w,b)=q+12[log⁡α-log⁡(2π)]-α2(‖w‖2+b2).
The score functions of w and b with respect to ℓα are, respectively(42)∂lα∂w=-αw,∂lα∂b=-αb.
The probability density function of this prior distribution can be added to the probability density function of the network error distribution. Then the updates for the parameters in a multilayer perceptron are(43)∂lϕk,α∂Ak=∑i=1nBik(yi-zir)(zik-1)T-αAk,∂lϕk,α∂bk=∑i=1nBik(yi-zir)-αbk.When the algorithm is overtraining, the parameters are still increasing. Subtractingα times the current parameter from the update effectively counteracts the update when the norm of that parameter is large enough. The algorithm will eventually reach a steady state where the parameters are longer increasing.SinceEpα{w}=0, Epα{b}=0, Epα{wwT}=(1/α)I, and Epα{b2}=1/α, the Fisher information matrix for the regularized algorithm is(44)Ĝϕk,α=1n-1∑i=1n(‖zik-1‖2+1)BikϵiϵiT(Bik)T+αI.AddingαI to the Fisher information matrix helps the algorithm avoid inverting an almost singular matrix. Thus, if all the errors are zero, then the Fisher information matrix will be αI. The effect will be multiplying all the weight update elements by 1/α. This is the largest value the inverse will reach during the SNGL algorithm's execution.In [13] Park adds regularization to the ANGL algorithm with an Akaike Information Criterion.


## 2.5.2. Annealed Learning Rates

Annealing the learning rate [6] is a technique in which a learning algorithm will slowly decay its learning rate after a specific epoch during the execution of the algorithm. Let the learning rate of the SNGL algorithm be ηt, where t is the current learning epoch. It is defined as(45)ηt=min⁡⁡(η0,τt),
where η0 is the initial learning rate and τ is a factor that determines which learning epoch begins the annealing process. The epoch in which the decay begins is t0=⌈τ/η0⌉.


## 3. Experimental Results

The exclusive OR function was the first Boolean function to be learned by a multilayer perceptron that was not linearly separable [4]. Encoding data with a low density parity check (LDPC) code [14] could be considered the exclusive OR problem's big brother. An LDPC is computed by multiplying (0,1)-vectors with a (0,1)-matrix over GF(2) to get a longer (0,1)-vector. Remember that in GF(2) addition is the XOR operation and multiplication is the AND operation. Formally, a low density parity check code is a code in which extra bits that are parity checks are added to the message. A matrix is used to represent the transformation of the bits over GF(2). The encoding function is the matrix(46)H=[00111100110110111010111001000001000001000001000001].
This code has a rate of 1/2 because the code words are twice as long as the message words.The training set consists of 32 training patterns. The inputs are the 32 possible messages, and the outputs are the 32 10-bit code words corresponding to each message. Thus, the network has 5 input units and 10 output units. The multilayer perceptron will also have 10 hidden units. During the training the bits are changed according to a bit error rate of 0.01. This will help the network deal with bit errors.First, the performance of the three algorithms will be compared. This will be followed by a comparison of their effective learning rates. Finally, an experiment in which the effective learning rate of the SNGL and ANGL algorithms when these learning rates are used in the ordinary gradient learning (OGL) algorithm as described in [4] will be discussed.


### 3.1. Performance Comparison

Figure1 shows the sum-squared error of the three algorithms when learning this LDPC encoding function. As can be seen, the performance of the SNGL algorithm is comparable to that of the OGL algorithm, but it reaches convergence in less than 1% of the time.Figure 1
Sum-squared error of the OGL, SNGL, and ANGL algorithms when learning the low density parity check (LDPC).Table2 lists the parameters used in all three algorithms and some of the performance metrics. The ANGL algorithm achieves a much lower MSSE. Using the prior probability by adding the current parameters multiplied by the hyperparameter α caused the components of the Fisher information matrix to approach infinity. Therefore, it was omitted, and the weights of the network were allowed to grow without constraint. Thus, the ANGL overtrained the network significantly.Table 2
LDPC experiment results.
LDPCOGLSNGLANGLHidden units101010Initial learning rate0.20.20.01Annealed learning rateNoYesNoPrior probability on parametersYesYesNoParameter of prior distribution0.0010.001N/AFisher information matrix sizeN/A10×10×2170×170Adaptation rateN/AN/A1/tNumber of extra parameters020028,900Final learning rate0.20.00972.6015Minimum sum-squared error0.35400.26910.00098Final learning epoch24,852105100,000Table2 may at first seem to be an unfair comparison. However, each of the algorithms were run with and without the regularization and annealed learning rates. The algorithm that performed the best is listed above. This was done to conserve space and simplify the reports. Thus, OGL performs better without an annealed learning rate on this problem. There were several difficulties implementing ANGL with the regularization. It did not perform well with regularization.


### 3.2. The Effective Learning Rate

The effective learning rate is defined to be the product of the baseline learning rate and the smallest eigenvalue of the matrix chosen to represent the inverse of the Fisher information matrix. For OGL the effective learning rate is just the learning rate chosen for the algorithm. For ANGL it is the learning rate multiplied by the smallest eigenvalue of the estimated Fisher information matrix. For SNGL it is the inverse of the largest eigenvalue of the Fisher information matrix multiplied by the learning rate. Figure2 shows the effective learning rates for these three algorithms. It is a simplified number to give the reader a sense of the step size used by the algorithm during each learning epoch.Figure 2
Effective learning rate of the OGL, SNGL, and ANGL algorithms when learning the LDPC encoding function.The effective learning rate for the SNGL algorithm is small for the first 100 learning epochs. Between the hundredth and the first thousandth epoch, the learning rate becomes quite large, for a step size, becoming larger than one for a few epochs. When the algorithm reaches convergence, then the learning rate begins to decrease due to annealing. The LDPC function is more complicated than XOR, and this is reflected in the number of learning epochs required to reach the three phases of learning for SNGL.


### 3.3. Substituting Learning Rates

One might argue that SNGL and ANGL perform better than OGL because they both use higher learning rates. However, as Figure2 shows, both of these algorithms have smaller learning rates at first, and the learning rates do not increase until progress is made towards good parameter values.Figure3 shows the performance of OGL with the effective learning rates of ANGL and SNGL in addition to OGL with a flat learning rate along with ANGL and SNGL. The plot shows that using the effective learning rate of ANGL improves the performance of OGL. However, the performance of OGL does not improve with the SNGL algorithm's effective learning rate. For this problem, the Fisher information matrix used in SNGL is a block-diagonal matrix with two 10-by-10 blocks on the main diagonal. In the ANGL algorithm, it is a 61-by-61 matrix. These matrices are positive definite. The eigenvalues represent how much information is associated with each direction. The eigenvectors associated with these eigenvalues are these directions. Using the maximum eigenvalue to determine an effective learning rate is good for demonstration purposes. However, it will not improve performance of the algorithm because the directional information is missing.Figure 3
Sum-squared error of the OGL, SNGL, and ANGL algorithms compared to the OGL algorithm with the effective LR of the ANGL and SNGL algorithms.


## 3.1. Performance Comparison

Figure1 shows the sum-squared error of the three algorithms when learning this LDPC encoding function. As can be seen, the performance of the SNGL algorithm is comparable to that of the OGL algorithm, but it reaches convergence in less than 1% of the time.Figure 1
Sum-squared error of the OGL, SNGL, and ANGL algorithms when learning the low density parity check (LDPC).Table2 lists the parameters used in all three algorithms and some of the performance metrics. The ANGL algorithm achieves a much lower MSSE. Using the prior probability by adding the current parameters multiplied by the hyperparameter α caused the components of the Fisher information matrix to approach infinity. Therefore, it was omitted, and the weights of the network were allowed to grow without constraint. Thus, the ANGL overtrained the network significantly.Table 2
LDPC experiment results.
LDPCOGLSNGLANGLHidden units101010Initial learning rate0.20.20.01Annealed learning rateNoYesNoPrior probability on parametersYesYesNoParameter of prior distribution0.0010.001N/AFisher information matrix sizeN/A10×10×2170×170Adaptation rateN/AN/A1/tNumber of extra parameters020028,900Final learning rate0.20.00972.6015Minimum sum-squared error0.35400.26910.00098Final learning epoch24,852105100,000Table2 may at first seem to be an unfair comparison. However, each of the algorithms were run with and without the regularization and annealed learning rates. The algorithm that performed the best is listed above. This was done to conserve space and simplify the reports. Thus, OGL performs better without an annealed learning rate on this problem. There were several difficulties implementing ANGL with the regularization. It did not perform well with regularization.


## 3.2. The Effective Learning Rate

The effective learning rate is defined to be the product of the baseline learning rate and the smallest eigenvalue of the matrix chosen to represent the inverse of the Fisher information matrix. For OGL the effective learning rate is just the learning rate chosen for the algorithm. For ANGL it is the learning rate multiplied by the smallest eigenvalue of the estimated Fisher information matrix. For SNGL it is the inverse of the largest eigenvalue of the Fisher information matrix multiplied by the learning rate. Figure2 shows the effective learning rates for these three algorithms. It is a simplified number to give the reader a sense of the step size used by the algorithm during each learning epoch.Figure 2
Effective learning rate of the OGL, SNGL, and ANGL algorithms when learning the LDPC encoding function.The effective learning rate for the SNGL algorithm is small for the first 100 learning epochs. Between the hundredth and the first thousandth epoch, the learning rate becomes quite large, for a step size, becoming larger than one for a few epochs. When the algorithm reaches convergence, then the learning rate begins to decrease due to annealing. The LDPC function is more complicated than XOR, and this is reflected in the number of learning epochs required to reach the three phases of learning for SNGL.


## 3.3. Substituting Learning Rates

One might argue that SNGL and ANGL perform better than OGL because they both use higher learning rates. However, as Figure2 shows, both of these algorithms have smaller learning rates at first, and the learning rates do not increase until progress is made towards good parameter values.Figure3 shows the performance of OGL with the effective learning rates of ANGL and SNGL in addition to OGL with a flat learning rate along with ANGL and SNGL. The plot shows that using the effective learning rate of ANGL improves the performance of OGL. However, the performance of OGL does not improve with the SNGL algorithm's effective learning rate. For this problem, the Fisher information matrix used in SNGL is a block-diagonal matrix with two 10-by-10 blocks on the main diagonal. In the ANGL algorithm, it is a 61-by-61 matrix. These matrices are positive definite. The eigenvalues represent how much information is associated with each direction. The eigenvectors associated with these eigenvalues are these directions. Using the maximum eigenvalue to determine an effective learning rate is good for demonstration purposes. However, it will not improve performance of the algorithm because the directional information is missing.Figure 3
Sum-squared error of the OGL, SNGL, and ANGL algorithms compared to the OGL algorithm with the effective LR of the ANGL and SNGL algorithms.


## 4. Conclusion

This paper developed a simplified natural gradient learning algorithm. Using an algebraic justification to form the gradient of a multilayer perceptron as a block matrix, a smaller block-diagonal Fisher information matrix was discovered. For the three problems studied in the experiments, the Fisher information matrix was much smaller than that in ANGL. The minimum sum-squared error performance of this algorithm is comparable to Amari et al. [2] but uses less memory and computational resources.



---

*Source: 407497-2011-07-24.xml*