# Kernel-Based Training of Generative Networks

Generative adversarial networks (GANs) are designed with the help of min-max optimization problems that are solved with stochastic gradient-type algorithms which are known to be non-robust. In this work we revisit a non-adversarial method based on kernels which relies on a pure minimization problem and propose a simple stochastic gradient algorithm for the computation of its solution. Using simplified tools from Stochastic Approximation theory we demonstrate that batch versions of the algorithm or smoothing of the gradient do not improve convergence. These observations allow for the development of a training algorithm that enjoys reduced computational complexity and increased robustness while exhibiting similar synthesis characteristics as classical GANs.

## Authors

• 6 publications
• 12 publications
• 1 publication
• ### SGD Learns One-Layer Networks in WGANs

Generative adversarial networks (GANs) are a widely used framework for l...
10/15/2019 ∙ by Qi Lei, et al. ∙ 44

• ### Training generative networks using random discriminators

In recent years, Generative Adversarial Networks (GANs) have drawn a lot...
04/22/2019 ∙ by Babak Barazandeh, et al. ∙ 0

• ### Designing GANs: A Likelihood Ratio Approach

We are interested in the design of generative adversarial networks. The ...
02/03/2020 ∙ by Kalliopi Basioti, et al. ∙ 44

12/26/2019 ∙ by Mingrui Liu, et al. ∙ 0

• ### Reducing Noise in GAN Training with Variance Reduced Extragradient

Using large mini-batches when training generative adversarial networks (...
04/18/2019 ∙ by Tatjana Chavdarova, et al. ∙ 38

• ### Integration of AI and mechanistic modeling in generative adversarial networks for stochastic inverse problems

The problem of finding distributions of input parameters for determinist...
09/17/2020 ∙ by Jaimit Parikh, et al. ∙ 0

• ### Approximation and convergence of GANs training: an SDE approach

Generative adversarial networks (GANs) have enjoyed tremendous empirical...
06/03/2020 ∙ by Haoyang Cao, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Background

Since their first appearance [2, 11], GANs have gained considerable attention and popularity, mainly due to their remarkable capability to produce, after proper training, synthetic data (usually images) that are realistically close to the data contained in their training set. The main challenge in designing GANs comes from the fact that their training algorithms require heavy computations that are primarily implementable on computationally powerful platforms. Such high computational needs arise not only because the size of the problems is usually large but also because the design of GANs requires the solution of min-max optimization problems. Stochastic gradient type algorithms employed for such cases very often exhibit non-robust behavior and slow rate of convergence, thus raising the computational needs considerably [7, 18].

In this work we focus, primarily, on the computational aspects of the training phase. Our intention is to develop a training algorithm which is simple and requires significantly less computations as compared to the current methods proposed in the literature and used in practice. In particular we will demonstrate, theoretically, that certain ideas as batch processing [17] and/or gradient smoothing [15] that are used for the solution of min-max problems have, in fact, absolutely no effect in the proposed formulation and, therefore, can be ignored. These conclusions will help us shape our algorithmic scheme and suggest a simple and efficient form. In order to be able to develop our algorithm we need to recall certain key notions from the theory behind GANs and understand how it differs from the alternative approach we adopt here.

Figure 1

captures the architecture employed during the training phase of GANs. There is a random vector

with unknownprobability density function (pdf) , with playing the role of a “prototype” random vector. The goal is to design a data-synthesis mechanism that generates realizations for the random vector . For this goal we employ a nonlinear transformation , known as the Generator, that transforms a random vector of known pdf (e.g. Gaussian or Uniform) into a random vector . We would like to design the parameters of the transformation so that is distributed according to . Under general assumptions such a transformation always exists [1, 6] and it can be efficiently approximated [8]

by a sufficiently large neural network, with

summarizing the network parameters.

Adversarial approaches in order to make the proper selection of employ a second nonlinear transformation that transforms and into suitable scalar statistics and and then compute a “mismatch” measure (not necessarily a distance) between the two random scalar quantities . The second transformation is also implemented with the help of a neural network, known as the Discriminator. We are interested in the average mismatch between namely which, after substitution, can be written as

 J(θ,ϑ)=EX,Y[J(D(X,ϑ),D(Y,ϑ))]=EX,Z[J(D(X,ϑ),D(G(Z,θ),ϑ))]. (1)

For every selection of the generator parameters we would like to select the discriminator parameters so that the average mismatch between is maximized. In other words we design the discriminator to differentiate between the synthetic random vector and the prototype random vector , as much as possible. This worst-case performance we then attempt to minimize by selecting properly the generator parameters . This leads to the following min-max optimization problem

 infθsupϑJ(θ,ϑ)=infθsupϑEX,Z[J(D(X,ϑ),D(G(Z,θ),ϑ))]. (2)

Common selections for the mismatch function are:

• , , see [11].

• , see [2].

• , see [2].

It is clear that the generator generates realizations of the random vector by transforming the realizations of . But how can we be assured that these realizations have the correct pdf namely ?. To see that this is indeed the case we need to consider the generator and discriminator transformations as being general functions not limited to the specific classes induced by the two neural networks. This immediately implies that by properly selecting we can shape the pdf of into any pdf we desire [1]. In this idealized situation optimizing over amounts to optimizing over and therefore over and, similarly, optimizing over amounts to optimization over . Consequently, we can redefine the min-max optimization problem in (2) under the following idealized frame

 (3)

where any scalar valued nonlinear transformation, and any pdf. The min-max problems corresponding to the three examples of we mentioned before accept analytic solutions. In particular in the first case, for fixed maximization over is attained for and the resulting functional is minimized over when . In the second case, assuming that , maximization over is achieved when and minimization over yields, again, . Similarly for the third case, maximization over is achieved for and minimization over when . As we can see all three min-max problems result in different optimum discriminator functions but agree in the final solution for , namely is shaped to have the same pdf as the prototype random vector .

Since in the original min-max problem (2) we limit the two transformations to be within the two classes induced by the input/output relationship of the corresponding neural network, it is clear that (2) constitutes an approximation to the ideal setup captured by (3). This implies that the output of the generator will follow a pdf which will be an approximation to the desired pdf of the prototype random vector . This approximation not only depends on the richness of the transformation class induced by the generator structure but, also, on the corresponding richness of the discriminator structure. As long as one of the two structures does not approximate sufficiently close the corresponding ideal transformation (when for example the neural network does not have sufficient number of layers), the design will fail in the sense that the realizations of will not follow the desired pdf of the prototype . The min-max optimization problem becomes more challenging because, as we mentioned, the pdf of is unknown and, instead, we are given a collection of independent realizations of (the training set) drawn from .

###### Remark 1.

Even though the goal is to design a generator network, with GANs we simultaneously require the design of an additional neural network, the discriminator. This requirement increases the number of parameters to be estimated considerably and, consequently, the computational complexity.

Furthermore the algorithmic solution of (2) relies on alternating stochastic gradient-type algorithms and the presence of two antagonistic optimization problems translates into an increased number of updates in the implementation which are also known to be non-robust [3, 7, 18].

Let us now see how we can accomplish a similar approximation for the output pdf of the generator without the need of a discriminator. We are going to revisit the idea suggested in [9, 12] which is based on kernel functions. Let be vectors of the same dimension of and consider a scalar function which is symmetric, i.e. and positive definite, namely, for every scalar function it satisfies

 ∬ϕ(U)k(U,V)ϕ(V)dUdV≥0,

with equality to 0 if and only if .

For two pdfs with fixed and to be determined, we define a distance measure as a function of as follows

 J(g)=∬(f(U)−g(U))k(U,V)(f(V)−g(V))dUdV. (4)

An immediate consequence of the positive definiteness property of the kernel is that the solution to the minimization problem is, obviously, . Let us now write the same distance using expectations. If , , are independent random vectors with following and following then, the double integral in (4) can be expressed as

 (5)

where is constant, not related to . Eq. (5) is simply an alternative way to rewrite the metric introduced in (4), consequently its minimization with respect to still results in the desired equality .

The next step consists in abandoning the ideal world expressed by (5). If is the output of the generator, this suggests that correspond to inputs . The two random input vectors must be statistically independent in order for the same property to be inherited by the two outputs . From (5), by substituting and using the symmetry of the kernel, we can define an average distance as a function of the generator parameters as follows111We prefer the symmetric form in (6) instead of the expectation of adopted in [9].

 J(θ)=EX,Z1,Z2[k(G(Z1,θ),G(Z2,θ))−k(G(Z1,θ),X)−k(G(Z2,θ),X)], (6)

where in the proposed measure we left out the constant term since it does not depend on . Performing the minimization generates a neural network whose output will have a pdf that approximates the desired pdf in the sense of the average distance we introduced in (4). This is clearly the equivalent of the min-max problem in (2) and, as we can see, it involves a pure minimization.

## 3 Properties of training algorithms

Because the problem we are concerned with involves only minimization, this allows for the use of classical stochastic gradient algorithms. What is also appealing, is that we have a rich arsenal of theoretical results coming from Stochastic Approximation Theory [4] that can support the training algorithm we intend to propose. In fact our goal is to arrive at an algorithmic scheme that has reduced computational complexity and demonstrate, theoretically and/or with simulations, that there is no significant performance loss in doing so. Actually the relevant properties we are going to use in order to propose our algorithmic scheme will be presented under a more general frame not limited to the specific optimization problem defined in the previous section.

### 3.1 Stochastic approximation

Suppose we are given a function where denotes a random vector for which we have available a sequence of independent realizations . Consider the following optimization problem

 infθEW[h(W,θ)] (7)

It is then well known that the stochastic gradient algorithm

 θt=θt−1−μH(Wt,θt−1),  H(W,θ)=∇θh(W,θ) (8)

where denotes the learning rate of the algorithm, can lead to a (local) minimizer of (7

) without knowing the probability distribution of

. The algorithm in (8) can be characterized [4] by the average trajectory and the corresponding random perturbations as where

 ¯θt=¯θt−1−μEW[H(W,¯θt−1)],

while for we have the steady-state description with the matrix satisfying the following Lyapunov equation

 CQ+QC⊺=EW[H(W,θ∗)H⊺(W,θ∗)]. (9)

Vector is the true (local) minimizer of (7) and is the Hessian of evaluated at . As it is explained in [4], the mean trajectory captures the transient phase while the perturbation part becomes leading during the steady-state of the corresponding algorithmic run. This is also graphically depicted in Figure 2

for the case of a simple regression model of the form , where has length 5, is i.i.d. zero-mean Gaussian with unit covariance matrix and

is i.i.d. additive zero-mean Gaussian noise with variance

and independent from . We are interested in

. In the stochastic gradient descent version

we select the learning rate . The average trajectory and the steady-state performance introduced before can be used as a means to compare algorithms.

###### Remark 2.

When two algorithms have similar average trajectories and exhibit the same steady-state behavior, they are practically equivalent in performance.

This simple rule which, of course, makes sense will allow us to examine whether certain alternative versions of the classical algorithm in (8) can indeed improve its convergence characteristics.

### 3.2 Does batch processing improve convergence?

A widespread impression [9, 17] is that if we use data in batches of length and approximate with instead of , then the corresponding algorithm

 θ′n=θ′n−1−μ′KK−1∑j=0H(WnK−j,θ′n−1) (10)

converges faster than (8). Considering speed in terms of iterations is actually completely unfair since each iteration in (10) involves the usage of vectors and gradient computations instead of a single and a single gradient computation in the classical scheme (8). In order for the comparison to be correct we need to count speed in terms of the number of vectors already used or the number of gradient computations already performed. Following this principle, (10) should be expressed as

 θ′nK=θ′(n−1)K−μ′KK−1∑j=0H(WnK−j,θ′(n−1)K). (11)

Returning to the classical version, from Stochastic Approximation theory we know that the algorithm in (8) has a natural ability for averaging/smoothing. This can become apparent if we subsample (8) every iterations and expand the formula across consecutive updates

 θnK=θ(n−1)K−μK−1∑j=0H(WnK−j,θnK−j−1). (12)

As we can see by selecting , (11) and (12) become very similar. Of course we observe that in the latter the parameter estimates are different in each term of the sum as opposed to the former where these estimates are all the same. We should however note that, since is very small, changes very slowly resulting in minor differences between and . Following Remark 2, we can make a formal claim by computing the average trajectories and the steady-state perturbation covariance matrices of the two versions. The following lemma compares the two algorithms.

###### Lemma 1.

The average trajectories in (11) and (12) are given respectively by

 ¯θ′nK =¯θ′(n−1)K−μ′EW[H(W,¯θ′(n−1)K)] ¯θnK

while the steady-state perturbation covariance matrix in both cases satisfies and is the solution of (9).

Proof: The first equation is a direct consequence of the definition of the average trajectory. For the second we assume sufficient smoothness of the vector function and apply a Taylor expansion around after expressing . Finally, the computation of the perturbation covariance matrices is also straightforward and since it involves the exact minimizer it is the same for both algorithms.

Lemma 1 implies that batching has actually no noticeable effect during initial convergence and during steady-state. This is also confirmed from Figure 3

where we present the relative estimation-difference power (red) of the two algorithms, again for the case of the simple regression model. As we can see this quantity is very small during the transient phase while during steady-state it becomes proportional to . The latter can be verified by performing simulations with different values and observing the corresponding change in the relative power during steady-state.

Actually, Figure 3 allows us to make a claim that is far stronger than Lemma 1

: Not only the two versions exhibit similar first and second order moments over iterations (i.e. similar average trajectories and steady-state perturbation covariance matrices, as stated in Lemma

1), but their actual estimates are very close to each other, provided of course that the two algorithms use the same data, synchronously.

###### Remark 3.

Even though no convergence speed improvement is observed, batch processing can be beneficial since it can exploit existing parallel or vectorized processing capabilities of the computational platform.

Indeed, as we demonstrated, per gradient computation there is no improvement in convergence, however, if there are parallel processing units (or vectorized computational capabilities) we can perform multiple gradient computations simultaneously and reduce the overall physical computational time [3].

### 3.3 Does smoothing improve convergence?

Another popular variation of the classical stochastic gradient algorithm consists in replacing the instantaneous gradient with a smoothed version updated over each iteration. In particular, instead of (8), in [15] it is proposed as alternative

 ~Ht =ρ~Ht−1+(1−ρ)H(Wt,θ′′t−1) (13) θ′′t =θ′′t−1−μ~Ht. (14)

The smoothing in (13) corresponds to an exponential windowing in place of the orthogonal window employed in the batch implementation. This can be seen from the expansion

 ~Ht=(1−ρ)t−1∑j=0ρjH(Wt−j,θ′′t−j−1).

A typical value of is 0.9 which implies that very quickly the contribution of past gradients in the sum, due to the term , becomes negligible and the sum appears as having practically a fixed number of terms. In fact it is commonly considered in signal processing that an exponential window has the same effect as an orthogonal window of length . This makes smoothing similar to batch processing and, therefore, it is expected not to provide any noticeable difference compared to the original algorithm (8).

We could offer a formal proof to our previous argument by finding, as before, the average trajectory and the steady-state covariance matrix and show that they are similar to the original version. Instead, for simplicity we provide a simulation example in the hope that it is equally convincing. In Figure 3 we plot the relative estimation-difference power (blue) for the regression model example where in (13) we used . As we observe, again the two algorithms provide similar estimates with the relative estimation-difference power being of the order of during steady-state and much smaller during the transient phase.

### 3.4 Does gradient computation using past parameter estimates affect convergence?

Next we would like to examine the effect on the algorithmic performance when in (8) the computation of the gradient is performed not by using but . In other words we consider the algorithm

 θ′′′t=θ′′′t−1−μH(Wt,θ′′′t−k). (15)

Again, computing the average trajectory and the steady-state perturbation covariance matrix we can show that the two algorithms in (8) and (15) are described by similar equations. In particular for (15) we need to use that fact that which will inflict an difference as compared to the average trajectory of (8). The steady-state behavior on the other hand will be the same. In Figure 3, as before, we plot (green) for a delay . As we can see, if for the computation of the gradient we use a delayed version of our parameter estimate, this has only a negligible effect on the overall convergence behavior of the algorithm.

###### Remark 4.

The previous properties apply to every algorithm in the form of (8). We should however emphasize that the computational schemes employed in classical GANs for solving (2) do not fall under this frame.

Indeed for min-max problems each update of (generator parameters) is followed by several updates of (discriminator parameters). Consequently batching/smoothing/delaying may affect these algorithms differently. As far as the class of algorithms captured by (8) is concerned, which are the focus of this work, we believe we have provided sufficient evidence that these modifications have no significant effect on the characteristics of the algorithm.

## 4 Proposed algorithmic scheme

Let us now return to the problem of interest, namely the minimization of which is defined in (6). Because of the properties described in Sections 3.2, 3.3 it is clear that we will adopt a simple version without smoothing (which is common in adversarial approaches [15]). The property mentioned in Section 3.4 will be used after we make the presentation of the first version of our algorithm and it will result in a significant computational reduction without any noticeable sacrifice in performance.

Following (7) and (8), at each iteration we need to provide two statistically independent realizations of the input vector and one realization of . As we pointed out can be generated since their pdf is assumed known (e.g. Gaussian or Uniform) while is available from the training data set .

We would like to point out that stochastic gradient type algorithms for the minimization of were previously proposed in [9, 16]. In [16] expectation is replaced by averaging over the whole set of available training data. Consequently each iteration requires a considerable amount of gradient computations. In [9] this problem is reduced since they propose the use of small-sized (mini) batches. Specifically they define one batch with ’s and a second with pairs . Because each from the first batch is combined with every pair in the second batch, the number of gradient evaluations is still elevated.

Consider now a neural network with two layers. In particular if are the input and output respectively, we define

 W=\mathpzcAZ+\mathpzca,  S=d(W);  T=\mathpzcBS+\mathpzcb,  Y=g(T)

where denote scalar nonlinearities with meaning that is applied to each element of the vector . Matrices and vectors constitute the parameters of the two layers with each matrix forming the linear combination and the vector providing the offset. The input vector is usually of much smaller dimension as compared to and and we reach the target dimension progressively. For our algorithm we consider only fully connected neural networks without any special constraint on their coefficients.

To apply the stochastic gradient algorithm, the most crucial part is the computation of the gradient of the kernel with respect to the parameter matrices and vectors. These parameters affect the kernel value through . Interestingly, there is a very simple recursive formula that allows for the computation of the corresponding derivatives. This is presented in the following lemma.

###### Lemma 2.

Denote with the gradients of the kernel with respect to the elements of the matrices that affect and with the partial derivatives arranged into a matrix of the same dimensions. To compute the gradients define

 \mathpzcV=g′(T)⊙∇Yk(Y,U), \mathpzcU=d′(W)⊙(\mathpzcB⊺V)

with denoting the derivatives of and denoting the element-by-element multiplication of two matrices of the same dimensions. Then the two gradients take the simple form

 ∇[\mathpzcA\mathpzca]k(Y,U)=\mathpzcU[S⊺1],   ∇[\mathpzcB\mathpzcb]k(Y,U)=\mathpzcV[Z⊺1].

Proof: The proof of this lemma presents no particular difficulty. It only requires careful housekeeping of the various partial derivatives. It is also worth mentioning that the gradients of the kernel function with respect to the parameters of each layer turn out to form a rank-one matrix.

Let us now limit ourselves to the Gaussian kernel case . We then have

 ∇Yk(Y,U)=−2he−1h∥Y−U∥2(Y−U).

The function that plays the role of in (7) is

 k(Y1,Y2)−k(Y1,X)−k(Y2,X)

for which we must compute the gradient with respect to the parameters of each layer. The first version of the algorithm is depicted in Table 1.

As we see from Table 1, in each iteration we need to compute the layer outputs and the corresponding gradients for two different inputs . It is exactly here that we intend to use the result of Section 3.4. We propose to compute the layer outputs and gradients corresponding to a single input while as second input we use the one generated during the previous iteration along with its outputs and gradients. This poses no problem since and are independent, as required by our analysis. Finally, a last issue we would like to address is normalization. In order for the learning rate to become data independent, it is necessary that the gradient, before being used in the update, to be normalized. For this reason we adopt the scheme proposed in [20] with . Incorporating all these points into our algorithmic structure produces the final result depicted in Table 2. Regarding our notation, if is a matrix then denotes its -th element.

## 5 Experiments

We applied our algorithm to the MNIST dataset. For the generator, which is the only neural network required by our approach, we used two fully connected layers (as mentioned full connectivity is the only structure considered in this work) with dimensions . Parameter was selected and in order to exploit the parallel processing capabilities of our computational platform we used a batch size of 32. Finally we applied a smoothing factor for the power estimation of each component of the gradient matrices, while for the learning rate we selected .

Similar geometry for the generator was adopted for the GAN implementation, namely two layers with dimensions while the discriminator structure was chosen to be . We also used a batch size of 32 with the same and . We should mention that in this case we also employed smoothing of the average gradient as suggested in [15] with a corresponding smoothing factor equal to . We plot the relative performance of the two methods in Figure 4 where we depict their Inception Score [5, 14] as a function of processing time. Our method attains better score values which is translated into more visually meaningful synthetic images as we can see from Figure 5.

We also applied our algorithm to the CelebA dataset (properly cropped). The geometry used for the generator was again two layers with dimensions , , batch size 64, . In Figure 6 we present results generated by the GAN architecture after (a) ; (b) ; (c)  iterations and batch size , while in (d) the results of our method after iterations and the same batch size. We must clarify that the iterations

of our method that led to Figure 6(d) required less physical time than the iterations with GANs which produced Figure 6(a). As we can see, our algorithm provides far superior synthetic images even when compared to the GAN results captured in Figure 6(c) and produced after the same number of iterations. The latter required five times more physical time than our method.

###### Remark 5.

We would like to mention that we also implemented our algorithm in Matlab. When executed on a laptop computer (MacBook Pro) in the case of the MNIST database, a single sweep of its 60,000 images required approximately 3 min. When, however, we exploited the vectorized computation capacity of Matlab and developed a batched version of the algorithm, the execution time was reduced to 4 sec with a batch of size 32.

We can experience truly exceptional improvement in image quality with GANs if we utilize deep convolutional networks [19] instead of the fully connected version previously discussed. The quality improvement can indeed be verified from Figure 7

where in (a) we repeat the results of our method from Figure 6(d) while in (b) we present the results of a deep convolutional GAN with geometry as the one suggested in [13]. We should however note that this amelioration in GANs comes at an extremely elevated computational cost. Indeed the physical time required for 200,000 iterations (with batch size ) was five times more than the 5 million iterations of our method.

It is also possible to measure quality by adopting the Maximum Mean Discrepancy score [12]. According to this index, better quality images produce smaller values. The scores for the three methods are: for the convolutional GANs, for the classical GANs, and finally for our approach. As we can see, these scores are in complete accordance with the visual perception of the corresponding synthetic images.

One might presume that, in our method, if we continue iterating this could lead to further improvement, possibly matching the quality of convolutional GANs. Unfortunately this is not the case due to saturation phenomena that occur when the parameter estimation method reaches its steady-state phase. Indeed the residual perturbations (explained in Section 3.1) prohibit any additional enhancement.

Finally, we would like to add that, currently our group is targeting the extension of our method to convolutional networks. The goal is to be able to avoid adversarial approaches by developing a non-adversarial training method which enjoys significant reduction in computational complexity and robustness in convergence. Basically, we would like the convolutional networks to inherit the same positive characteristics established for their fully connected counterparts.

## 6 Acknowledgment

This work was supported by the US National Science Foundation under Grant CIF 1513373, through Rutgers University.

## Appendix A Appendix

In the Appendix we provide the Matlab code for the algorithm in Table 2. We recall that this algorithm performs training of a fully connected two-layer neural network. If one chooses to use this program, the training dataset must be in a .MAT file in the form of a single matrix called data. The columns of this matrix must contain the training vectors. With this program we exploit Matlab’s vectorized computations and propose a batched version of the algorithm appearing in Table 2. The corresponding speedup in execution, as was mentioned in Remark 5, is significant even for a single processor platform. We believe that with such execution times it is no longer unrealistic to perform training of generative networks on laptops using Matlab.

function [A,a,B,b]=kerntrain(data_file,n,m,h,mu,rounds,batch);

%% Input Parameters
% file_name: is a string containing the data file name
% n: input size
% m: first layer output size.
% h: is the parameter of the Gaussian Kernel. (Typical value 80)
% mu: learning rate. (Typical value 0.0001)
% rounds: is how many times we would like to recircle the elements of the
%         database. (Typical value 50)
% batch: batch size. (Typical value 32 or 64)
%
%% Output Parameters
% A,a: first layer parameters
% B,b: second layer parameters

% Data must be in the file: data_file.mat. The file must contain a
% single matrix. This matrix must be called:     data
%
% Each matrix column corresponds to a different image reshaped into a
% (column) vector.

%% RMSprop (smoothing) parameter
lambda = 0.999;
e = 10^-8; % small number to be used to avoid divisions by 0.

%% Initialization following Clorot and Bengio, 2010
times = size(data,2); % length of database
k = size(data,1); % final output size
A = randn(m,n)/sqrt(m/2);
a = zeros(m,1);
B = randn(k,m)/sqrt((k+m)/4);
b = zeros(k,1);
Z0 = zeros(n,batch);
S0 = zeros(m,batch);
Y0 = zeros(k,batch);
V0 = zeros(k,batch);
U0 = zeros(m,batch);
M  = zeros(k,m+1);
N  = zeros(m,n+1);

%% Main Algorithm following Table 2, (batched version)
for r=1:rounds
ZZ = randn(n,times); % generate all inputs for one epoch
for t = batch:batch:times
Z = ZZ(:,t-batch+1:t); % select batch inputs
X = data(:,t-batch+1:t); % select batch training vectors
% Compute batch outputs of the two layers
W = A*Z + repmat(a,1,batch);
S = max(W,0);   % ReLU
T = B*S + repmat(b,1,batch);
Y = 1./(1 + exp(-T));  % Sigmoid
YX = Y - X;
YY0  = Y - Y0;
KernYX  = exp(-sum(YX.^2,1)/h);
KernYY0 = exp(-sum(YY0.^2,1)/h);
R = repmat(KernYX,k,1).*YX - repmat(KernYY0,k,1).*YY0;
V = (Y - Y.^2).*R;
U = max(sign(W),0).*(B’*V);
G = V*[S’ ones(batch,1)] + V0*[S0’ ones(batch,1)];
D = U*[Z’ ones(batch,1)] + U0*[Z0’ ones(batch,1)];
% Compute average power of gradient elements
if (t==batch)&&(r==1) % if the first time, don’t smooth
M = G.^2;
N = D.^2;
else  % otherwise smooth
M = lambda*M + (1-lambda)*G.^2;
N = lambda*N + (1-lambda)*D.^2;
end
% Update network parameters
B = B - mu*G(:,1:end-1)./sqrt(M(:,1:end-1)+e);
b = b - mu*G(:,end)./sqrt(M(:,end)+e);
A = A - mu*D(:,1:end-1)./sqrt(N(:,1:end-1)+e);
a = a - mu*D(:,end)./sqrt(N(:,end)+e);
% Update variables needed for the next iteration
Z0 = Z;
S0 = S;
Y0 = Y;
V0 = V;
U0 = U;
end
end