Year: 2021

Abstract:
The inclusion of second-order information into Deep Learning optimization has drawn consistent interest as a way forward to improve upon gradient descent methods. Estimating the second-order update is computationally expensive, which drastically limits its usage scope and forces the use of various truncations and approximations. This work demonstrates that it is possible to solve the Newton direction in the stochastic case exactly. We consider feedforward networks as a base model, build a second-order Lagrangian which we call Sifrian, and provide a closed-form formula for the exact stochastic Newton direction under some monotonicity and regularization conditions. We propose a convexity correction to escape saddle points, and we reconsider the intrinsic stochasticity of the online learning process to improve upon the formulas. We finally compare the performance of the developed solution with well-established training methods and show its viability as a training method for Deep Learning. Optimization in Deep Learning is mainly dominated by first-order methods built around the central concept of backpropagation (LeCun et al., 1988). Second-order methods have exceptional theoretical properties in deterministic optimization (Boyd et al., 2004; Nocedal & Wright, 2006), but such properties do not translate well into the stochastic case (LeCun et al., 2012; Bottou et al., 2018). First-order methods such as the Stochastic Gradient Descent (SGD) (Robbins & Monro, 1951) are relatively simple and have been adaptively improved for Deep Learning (Duchi et al., 2011; Tieleman & Hinton, 2012; Kingma & Ba, 2014; Reddi et al., 2019; Yao et al., 2020). Inherently, secondorder methods might seems inadequate for Deep learning due to increased computational cost, poor clock-wall performance, and the non-convex nature of Deep Learning (LeCun et al., 2012). Furthermore, Newton method might even reduce the generalization capabilities of training (Wadia et al., 2021; Amari et al., 2020). Despite these various limitations, substantial effort was deployed to include Hessian information into the optimization process (e.g. Byrd et al., 2011; Sohl-Dickstein et al., 2014; Byrd et al., 2016; Agarwal et al., 2017; Berahas et al., 2019; Anil et al., 2020; Goldfarb et al., 2020; Castera et al., 2021). Several approaches exist such as the Gauss-Newton method (e.g Schraudolph, 2002; Botev et al., 2017), diagonal approximation of the Hessian (e.g. Bordes et al., 2009; Schaul et al., 2013), iterative low-rank updates such as BFGS (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970) (see also, Liu & Nocedal, 1989; Schraudolph et al., 2007; Bollapragada et al., 2018), or Hessian-Free methods which combine the fast Hessian matrix multiplication (Pearlmutter, 1994) and the conjugate gradient algorithm (e.g. Martens, 2010; Martens & Sutskever, 2012; Dauphin et al., 2014). The use of the Fisher matrix to capture curvature information in the space of distributions, instead of the Hessian, is another approach which yields the natural gradient (Amari, 1998), but suffers from the same computational issues as the Newton method. Within this context K-FAC method (Martens & Grosse, 2015; Ba et al., 2016a; George et al., 2018) mitigates some of the computational issues of the natural gradient method. In essence, several drawbacks and flaws limits the adoption of second-order methods as a standard for neural networks training . In this paper, we propose a novel approach to characterize the Newton update which help us derive an exact closed-form solution for the stochastic Newton method. Our method requires a suitable regularization of the neural network and strict monotonicity of the activation functions. We start this paper by introducing useful notations while deriving the well-known backpropagation algorithm for a feedforward network. Then we introduce a second-order Lagrangian which we call Sifrian, that will serve to characterize the Newton direction. We derive four types of equations from the Sifrian, and we provide an exact closed-form solution for the Newton direction in the stochastic case. We

text: The inclusion of second-order information into Deep Learning optimization has drawn consistent interest as a way forward to improve upon gradient descent methods. Estimating the second-order update is computationally expensive, which drastically limits its usage scope and forces the use of various truncations and approximations. This work demonstrates that it is possible to solve the Newton direction in the stochastic case exactly. We consider feedforward networks as a base model, build a second-order Lagrangian which we call Sifrian, and provide a closed-form formula for the exact stochastic Newton direction under some monotonicity and regularization conditions. We propose a convexity correction to escape saddle points, and we reconsider the intrinsic stochasticity of the online learning process to improve upon the formulas. We finally compare the performance of the developed solution with well-established training methods and show its viability as a training method for Deep Learning.
Optimization in Deep Learning is mainly dominated by first-order methods built around the central concept of backpropagation (LeCun et al., 1988). Second-order methods have exceptional theoretical properties in deterministic optimization (Boyd et al., 2004; Nocedal & Wright, 2006), but such properties do not translate well into the stochastic case (LeCun et al., 2012; Bottou et al., 2018). First-order methods such as the Stochastic Gradient Descent (SGD) (Robbins & Monro, 1951) are relatively simple and have been adaptively improved for Deep Learning (Duchi et al., 2011; Tieleman & Hinton, 2012; Kingma & Ba, 2014; Reddi et al., 2019; Yao et al., 2020). Inherently, secondorder methods might seems inadequate for Deep learning due to increased computational cost, poor clock-wall performance, and the non-convex nature of Deep Learning (LeCun et al., 2012). Furthermore, Newton method might even reduce the generalization capabilities of training (Wadia et al., 2021; Amari et al., 2020). Despite these various limitations, substantial effort was deployed to include Hessian information into the optimization process (e.g. Byrd et al., 2011; Sohl-Dickstein et al., 2014; Byrd et al., 2016; Agarwal et al., 2017; Berahas et al., 2019; Anil et al., 2020; Goldfarb et al., 2020; Castera et al., 2021). Several approaches exist such as the Gauss-Newton method (e.g Schraudolph, 2002; Botev et al., 2017), diagonal approximation of the Hessian (e.g. Bordes et al., 2009; Schaul et al., 2013), iterative low-rank updates such as BFGS (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970) (see also, Liu & Nocedal, 1989; Schraudolph et al., 2007; Bollapragada et al., 2018), or Hessian-Free methods which combine the fast Hessian matrix multiplication (Pearlmutter, 1994) and the conjugate gradient algorithm (e.g. Martens, 2010; Martens & Sutskever, 2012; Dauphin et al., 2014). The use of the Fisher matrix to capture curvature information in the space of distributions, instead of the Hessian, is another approach which yields the natural gradient (Amari, 1998), but suffers from the same computational issues as the Newton method. Within this context K-FAC method (Martens & Grosse, 2015; Ba et al., 2016a; George et al., 2018) mitigates some of the computational issues of the natural gradient method. In essence, several drawbacks and flaws limits the adoption of second-order methods as a standard for neural networks training .
In this paper, we propose a novel approach to characterize the Newton update which help us derive an exact closed-form solution for the stochastic Newton method. Our method requires a suitable regularization of the neural network and strict monotonicity of the activation functions. We start this paper by introducing useful notations while deriving the well-known backpropagation algorithm for a feedforward network. Then we introduce a second-order Lagrangian which we call Sifrian, that will serve to characterize the Newton direction. We derive four types of equations from the Sifrian, and we provide an exact closed-form solution for the Newton direction in the stochastic case. We
further propose a saddle-free version of our method, and add a randomization process to enhance our solution. In the last part of this paper, we show the applicability of our method for Deep Learning through diverse classification tasks using feedforward architectures.

heading: 1 PRELIMINARIES
text: Deep Learning and neural networks training could be seen as an optimization problem of an expected loss function ` over a distribution D of labeled samples (x,d). In general, weights and biases (W,β) are the sought-after parameters of the network. The labeled samples distribution D is often unknown, but a large number of samples (database D) allows the approximation of the expected loss with an empirical risk:
min W ,β
( E(x,d)∼D [` (W ,β,x,d)] ) → 1 |D| ∑ p∈D ` ( W ,β,x(p),d(p) ) , ( x(p),d(p) ) p∈D i.i.d∼ D. (1)
The loss or cost function ` is typically a cross entropy function or an l2 norm of the mismatches between the network outputs and the labels: ` ( W ,β, ( x(p),d(p) )) = 12 ∥∥∥x(p)n − d(p)∥∥∥2 2 . Where n indicates the output layer number. We present more in detail the architecture of a Feedforward Neural Network (FNN) hereafter.

heading: 1.1 NOTATIONS AND FEEDFORWARD NEURAL NETWORKS (FNN)
text: In this section, we recall the main details of FNNs, which will be used as a standard model for Deep Learning1. The notation used throughout this paper is similar to notations presented in LeCun et al. (1988). The main equation governing the FNN, a.k.a. the forward model is the following:
x (p) k = F ( Wk x (p) k−1 + βk ) , k ∈ [1..n] , p ∈ D. (2)
where D is the database, p designs one single sample from the database (e.g. one single image or audio recording), k is the layer number and n is the total number of layers in the network. The initial input for sample p is x(p)0 (e.g. the vectorized input image data). The state variable x
(p) is transformed at each layer k through a multiplication by a weight matrix Wk and an addition of a bias vector βk. The activation function F , which is typically a sigmoid or a Rectified Linear Unit (ReLU), is applied element wise on the resulting activation vector: a(p)k = Wk x (p) k−1 + βk, and serves to introduce non-linearity in the neural network.

heading: 1.2 THE LAGRANGIAN AND BACKPROPAGATION
text: The origin of backpropagation could be traced back to the early 1970s and could be derived by casting Deep Learning as a constrained optimization problem. This section is similar to LeCun et al. (1988). The main novelty in this section is the addition of a regularization term dependent on the state variable R(x) to the cost function `: `R = ` + R (x). The regularization has to verify the admissibility criteria which we define hereafter:
Definition 1: A regularization term is admissible if the second derivative of the augmented cost function w.r.t. to the state variable is separable and non-null, i.e.:
∀ (p, q) ∈ D, ∀ (k,m) ∈ [1..n]
( ∂`R
∂x (p) k x (q) m ) p,q,k,m ∝ δp=q. (3)
An example of such regularization is the following function:
R(x) = 1
2 ∑ p∈D ∑ k=1..n−1 〈 x (p) k ,Λ (p) k x (p) k 〉 (4)
1Convolutional Neural Networks (CNN) are a type of FNN. The concepts of this paper apply also to CNN.
( Λ
(p) k ) p,k are symmetric positive matrices. Such a choice might seem atypical since the regularization often concerns the network parameters, mainly the weights or the biases. This ”admissible” regularization was introduced for the sole purpose of guaranteeing a non-null partial derivative of the cost ` w.r.t the state variable x(p)k . This last property is essential to solve Equ. 16. We will show later that such a regularization minimizes also the norm of the gradient (Barrett & Dherin, 2020; Smith et al., 2021). In order to derive the backpropagation algorithm we introduce the following Lagrangian as per the notation of LeCun et al. (1988):
L(x,W ,β,b) = `R + ∑ p∈D ∑ k=1..n 〈 x (p) k − F ( Wk x (p) k−1 + βk ) ,b (p) k 〉 . (5)
The Lagrangian contains the original cost function, the admissible regularization term, and the product of the forward equation with the adjoint state vectors ( b
(p) k ) p,k . These adjoint vectors are de-
fined for each layer of the network and each sample of the database; however, their values are not fixed yet and will be chosen in a way that simplifies the gradient computation. If the state variable x verifies the forward equations, then the Lagrangian simplifies to the cost function:
L(x (W ,β) ,W ,β,b) = `R (x (W ,β)) . (6) The total derivative of the previous Lagrangian w.r.t. the weights or biases is the gradient and could be expressed in terms of partial derivatives as follows2: dL(x(W ,β),W ,β,b) dWk = ∂L∂Wk + ∑ p∈D ∑ m=1..n ( dx(p)m dWk )( ∂L ∂x (p) m ) = d`RdWk ,
dL(x(W ,β),W ,β,b) dβk = ∂L∂βk + ∑ p∈D ∑ m=1..n ( dx(p)m dβk )( ∂L
∂x (p) m
) = d`Rdβk .
(7)
The partial derivatives of the Lagrangian w.r.t to {Wk,βk}k=1..n are straightforward to compute. The Fréchet derivatives ( dx(p)m dWk , dx(p)m dβk ) p,k,m are non-trivial to evaluate and the core idea of back-
propagation is the selection of the adjoint state variables ( b
(p) k ) p,k which cancels the superfluous
terms. Such a simplification is achievable if:( ∂L ∂x
(p) k ) p∈D, k∈[1..n] = 0. (8)
Computing the previous partial derivative yields the following backpropagation equation: ∂L ∂x
(p) k
= b (p) k − 1k=1..n−1W T k+1∇F
( a
(p) k+1
) b
(p) k+1 +
∂`R
∂x (p) k
= 0. (9)
1E is the indicator function, and it is equal to one if the underlying condition E is true and null otherwise. The nabla operator ∇F is a diagonal square matrix with an element-wise derivative of its argument (along the diagonal). The resolution of the backpropagation system could be split into a boundary condition and a backward propagation system. Further details could be found in LeCun et al. (1988). The particular choice of the adjoint state vectors ( b
(p) k ) p∈D, k=1..n yields a simple
formula for the gradient of the cost function:Gk = ∂L ∂Wk = − ∑ p∇F ( a (p) k ) b (p) k x (p)T k−1 ,
gk = ∂L ∂βk
= − ∑ p∇F ( a (p) k ) b (p) k .
(10)
The backpropagation algorithm is fundamental for Deep Learning. It is simple and provides a straightforward solution to an otherwise tedious problem to solve. Unfortunately, developing an efficient second-order backpropagation method remains elusive. The current state-of-the-art is based on the use of the R−operator (Pearlmutter, 1994) to compute the product of the Hessian with a given vector without explicitly calculating or storing the Hessian (Martens, 2010; Dauphin et al., 2014; Agarwal et al., 2017). In the following section, we provide a different framework for backpropagation which allows the characterization of the second-order Newton direction.
2A denominator layout notation for derivation is used throughout this paper.

heading: 2 THE SIFRIAN AND THE EXACT NEWTON DIRECTION
text: It is legitimate to wonder if the Lagrangian multipliers method could be replicated for higher-order optimization. We hereafter provide a solution for the second-order case and build an equivalent Lagrangian which we call the Sifrian.

heading: 2.1 THE SIFRIAN
text: The Lagrangian is useful to derive the first-order gradient, and it combines the cost function and the forward model multiplied by adjoint variables. Our approach to extend the Lagrangian consists in creating a function with the forward model, the backward model, and the definition of the gradient. We choose to omit the loss function, as its derivation will generate the gradient again. Starting from this paradigm, we introduce the Sifrian3.
Definition 2: The Sifrian S of the feed-forward network is defined as follows: S(x,b,W ,β,G,g,γ, ζ,N ,η) = ∑
k,p
〈 x
(p) k − F
( a
(p) k
) ,γ
(p) k 〉 + 〈 b
(p) k − 1k=1..n−1W Tk+1∇F
( a
(p) k+1
) b
(p) k+1 + ∂`R ∂x
(p) k
, ζ (p) k 〉 + ∑
k
〈 Gk + ∑ p∇F ( a (p) k ) b (p) k x (p)T k−1 ,Nk 〉 + 〈 gk + ∑ p∇F ( a (p) k ) b (p) k ,ηk 〉 .
(11)
The Sifrian includes all the equations generated by the first order backpropagation, multiplied by four new adjoint parameters { γ
(p) k , ζ (p) k ,Nk,ηk
} . While the Lagrangian describes a saddle
point (Nocedal & Wright, 2006), the Sifrian, describes an equilibrium and it has a unique null value when all the forward, backward and gradients equations are verified i.e.:
S(x (W ,β) ,b (W ,β,x) ,W ,β,G (W ,β,x,b) ,g (W ,β,x,b) ,γ, ζ,N ,η) = 0. (12)
We call the previous quantity the equilibrated Sifrian, and it is critical to notice its invariance to total derivation w.r.t. either the weights or biases (i.e., remains null). Expressing the total derivative of the equilibrated Sifrian using partial derivatives yields the following two equations: dS dWk = ∂S∂Wk + ∑ m,p dx(p)m dWk ∂S ∂x (p) m + db(p)m dWk ∂S ∂b (p) m + ∑ m dGm dWk ∂S ∂Gm + dgmdWk ∂S ∂gm = 0,
dS dβk = ∂S∂βk + ∑ m,p dx(p)m dβk ∂S ∂x (p) m + db(p)m dβk ∂S ∂b (p) m + ∑ m dGm dβk ∂S ∂Gm + dgmdβk ∂S ∂gm = 0.
(13)
Similar to first order backpropagation, we select the new adjoints variables{ γ
(p) k , ζ (p) k ,Nk,ηk } p∈D, k=1..n to circumvent the calculation of any superfluous terms. We
make the following choice:
∀m,∀ q : ∂S ∂Wm = −Gm, ∂S ∂βm = −gm, ∂S ∂x (q) m = 0, ∂S ∂b (q) m = 0. (14)
In this case, after vectorizing the weight gradient Gk (and also Nk), and concatenating all their values across all layers, the Equ. 13 could be cast as the following block matrix system:[
dG dW dg dW dG dβ dg dβ ] [ N η ] = [ G g ] →H [ N η ] = [ G g ] . (15)
3The choice of the Sifrian appellation stems from the Arabic word ”Sifr”, which means zero. ”Rien” means also ”nothing” in French.
The matrix that appears in the previous system (H) is the Hessian. Hence the computation of the Newton direction could be achieved by solving the associated system of four equations in Equ. 14. The derivation of the four types of Sifrian equations4 is provided with further details in the Appendix (A.1). We report hereafter the equations for piece-wise affine activation functions (e.g., ReLU or Leaky ReLU). The general case of C2 activation function is included in the Appendix (A.1).
(1) ∑
p∇F (p) k
( γ
(p) k x (p)T k−1 + b (p) k x (p)T k−1 + 1k=2..nb (p) k ζ (p)T k−1
) = 0,
(2) ∑
p∇F (p) k
( γ
(p) k + b (p) k
) = 0,
(3)γ (p) k − 1k<nW Tk+1∇F (p) k+1γ (p) k+1 + ∑ m ∂2`R ∂x
(p) k ∂x (p) m
ζ (p) m + 1k<nN T k+1∇F (p) k+1b (p) k+1 = 0,
(4) ζ (p) k − 1k=2..n∇F (p) k Wkζ (p) k−1 +∇F (p) k ( Nkx (p) k−1 + ηk ) = 0.
(16)
Comparison with the R-operator: Characterizing the Hessian effect without explicitly computing its values has been proposed by Pearlmutter (1994), using an operator R which requires the preselection of a specific direction to study its transformation through the Hessian. The Sifrian formulation can achieve the same goal by relaxing the two gradient constraints from the Sifrian equations and preselecting the values of {N ,η}. The product requires a forward pass to compute ζ, a backward pass to compute γ. Finally, the sought-after product is provided by the gradient of the Sifrian w.r.t. to the weights and biases. The Sifrian formulation is derived differently from the R-operator. Nevertheless, it can still fulfill the same role, and it is more suitable for inversion, which we will conduct hereafter.

heading: 2.2 THE EXACT STOCHASTIC NEWTON (ESN)
text: The resolution of the previous system is intricate in the general case. However, if we consider a single sample, i.e., the stochastic case and a strictly monotonous activation function, then a closedform solution could be derived.
ESN Theorem: If the cost function is admissibly regularized, the activation function F is strictlymonotonous piece-wise affine, the first order adjoints ( b
(p) k ) p,k are non-null and the curvature of
the output layer {
∂2`R ∂x (p) n ∂x (p) n
} is invertible, then the minimal-rank stochastic Newton second-order
update has a closed-form solution: Nk = − ∇F ( a (p) k ) b (p) k∥∥∥∇F(a(p)k )b(p)k ∥∥∥2 2 ( ∂`R ∂x (p) k−1 )T ,
ηk = 1k=n
[ ∇F ( a (p) n )]−1 { ∂2`R
∂x (p) n ∂x (p) n }−1 ∂`R ∂x (p) n −Nkx(p)k−1.
(17)
Proof: Under the proposed assumptions, the Sifrian second equation simplifies to γ(p)k = −b (p) k , which starts a cascade of simplifications. The first equations becomes 1k=2..nζ (p) k−1 = 0. The last layer of the third equation yields: γ(p)n + ∂ 2`R
∂x (p) n ∂x (p) n ζ (p) n = 0. The remaining two equations become: 1k=1..n−1 ∂`R ∂x (p) k + 1k=1..n−1N T k+1∇F ( a (p) k+1 ) b (p) k+1 = 0, − {
∂2`R ∂x (p) n ∂x (p) n }−1 ∂`R ∂x (p) n +∇F ( a (p) k )( Nkx (p) k−1 + ηk ) = 0.
(18)
4The Sifrian yields a total of 2 |D|+ 2 equations with an equal number of unknowns.
The Newton weight is a matrix, and it is determined through its product with the adjoint vector. Inverting the Newton direction is a degenerate problem, but the minimal rank solution is unique, and it is given in Equ. 17. A major issue still remains: the results of the Newton weight do not extend to the first layer (Equ. 18). To circumvent this shortcoming, we introduce an additional layer before the input, which corresponds to the identity multiplication and a null bias. This layer remains unchanged throughout the optimization process, and we call it the ”white layer”. This pre-input extension allows the definition of an update rule for the first layer.
White layer concept The third equation of the Sifrian allows the characterization of the Newton weight Nk only between the 2nd and the output (nth) layers. Such a problem is fundamental and it is legitimate to wonder how is it possible for Hessian-Free methods to work for inversion. The answer lays within the regularization/damping (Marquardt, 1970) which covers this structural gap. The important role of the first layer for the Newton method has been highlighted recently by Wadia et al. (2021). In our case, and in order to extend the inversion result to the first layer, we just introduce an additional layer before the input, which has the same size as the input, and does not bring any modification, hence the appellation ”white layer”. The white layer corresponds to a matrix multiplication with the identity, a null bias, an identity activation function and is unchanged during the training.
Admissible regularization: The ESN solution relies primarily on the existence of the term( ∂`R ∂x
(p) k
) which would be null if not for the admissible regularization, whose effect revolves ul-
timately around minimizing the norm of the state variable ( x
(p) k ) p,k . Given that the cost function `
steers the weights/biases to minimize the mismatches, and consequently the adjoints ( b
(p) k ) p,k then
the admissible regularization is contributing to the gradient’s norm minimization (Equ. 10), which is a desirable features for generalization purposes (Barrett & Dherin, 2020; Smith et al., 2021).
The ”granular” solution: The admissible regularisation could be removed and this would lead to a ”granular” solution which updates only the bias of the output layer ηn. The granular solution is not practical and it is of no use for the neural network training.
Hypotheses relaxation: Some of the hypotheses of the ESN theorem could be relaxed: notably the invertibility of the output curvature: { ∂2`R
∂x (p) n ∂x (p) n
} . Moore-Penrose inverse could be used instead.
The output curvature is often positive definite and its inversion is simple for e.g. it is the identity for an l2 cost function.

heading: 3 ADJUSTMENT OF THE EXACT STOCHASTIC NEWTON (ESN)
text: The ESN solution is exact and simple to compute; however, it diverges if implemented as it is. In this section, we first address the non-convexity issue and adjust the solution to become a descent direction. We also explore the stochasticity effect stemming from the samples randomness during the training process and provide an enhancement to the ESN. For the remainder of this paper we assume also that regularization is a of the same form as Equ. 4. We start this section by correcting the concave directions of the ESN.

heading: 3.1 ESCAPING THE SADDLES: A DESCENT DIRECTION
text: The minimal rank solution (Equ. 17) has some similarities with the gradient expression (Equ. 10), yet it is not a descent direction. To further illustrate this, we consider the inner product of the Newton and the gradient expressions layer by layer. For the weights, the product is always positive:
〈Nk,Gk〉 = 〈 x (p) k−1,Λ (p) k−1x (p) k−1 〉 = ∥∥∥x(p)k−1∥∥∥2
Λ (p) k−1
≥ 0. (19)
However, most of the Newton biases corresponds to ascent directions :
〈ηk,gk〉 = − ∥∥∥x(p)k−1∥∥∥2
Λ (p) k−1
+ 1k=n ∥∥∥∥ ∂`R ∂x (p) n ∥∥∥∥2( ∂2`R
∂x (p) n ∂x (p) n )−1 . (20) The last layer bias nature depends on the values of the regularization and the error. One of the main advantages of the Sifrian formulation is that we can treat the weights and biases separately and we can exactly pinpoint to the origin of non-convexity: the bias update. We choose to build our saddle-free variation by just changing the signs of the ascent direction. This manipulation is straightforward for all layers except the output bias update. The problematic term corresponds to the ”granular” solution (i.e., no regularization), which overfits every sample. We chose to keep the ”granular” solution as is and switch the sign of the remaining term. This guarantees that the update is a descent direction5. The saddle-free ESN is the following: Nk = − ∇F ( a (p) k ) b (p) k∥∥∥∇F(a(p)k )b(p)k ∥∥∥2 2 ( ∂`R ∂x (p) k−1 )T ,
ηk = 1k=n
[ ∇F ( a (p) n )]−1 { ∂2`R
∂x (p) n ∂x (p) n }−1 ∂`R ∂x (p) n +Nkx (p) k−1.
(21)
The above convexity correction is significantly simpler than the currently existing methods such as the Gauss Newton truncation (e.g. Schraudolph, 2002; Martens, 2010) or the absolute Hessian (Dauphin et al., 2014). The convexity correction is the first step in improving the ESN. Another intrinsic issue in the above formulation is that we try to optimize a given sample without taking into consideration the yet ”unknown” coming samples. The essence of stochasticity in still missing. In the next section, we include such an aspect into the ESN method.

heading: 3.2 RANDOMIZATION
text: The ESN solution we have proposed so far is exact but remains deterministic since it focuses on training a single pre-selected sample p. An intrinsic source of stochasticity in online learning stems from the randomness of the streamed samples. As such, it is preferable to use random variables to describe the dynamics of ESN. The two fundamental equations of ESN after discarding the ”granular” solution, adding the white layer are: ( ∂`R ∂x (p) k−1 )T + y (p)T k Nk = 0, with y (p) k = ∇F ( a (p) k ) b (p) k ,
Nkx (p) k−1 + ηk = 0.
(22)
5This point could be further justified with a first order Taylor development of the loss function.
Let’s consider Ω = D as the set of possible outcomes of the training samples6 and randomize the sample index p → p. Aside from {N ,η} all the vectors in the above system become random variables. We solve the randomized Equ. 22 using the following least-squares criteria:
argmin{N ,η} ∫ Ω=D dΩ (p) ∥∥∥∥∥∥ ( ∂`R
∂x (p) k−1
)T + y
(p)T k Nk ∥∥∥∥∥∥ 2
2
+ ∥∥∥Nkx(p)k−1 + ηk∥∥∥2
2 . (23)
We changed the ESN into a least-squares problem, which is desirable, given the sizeable wellestablished literature about optimal unbiased estimators for this class of problems (Aitken, 1936; Kariya & Kurata, 2004). The minimisation could be done for each layer separately, which reduces the inverse problem dimensionality. The search for critical points is detailed in the Appendix (A.4), and yields the following randomized ESN system which involves a Sylvester equation (forNk) : Ep [ y (p) k y (p) k T ] Nk +NkVar ( x (p) k−1 ) = −Ep [ y (p) k ( ∂`R ∂x (p) k−1 )T] , NkEp [ x (p) k−1 ] + ηk = 0. (24)
We notice that the above system could considerably simplify if the state variables x(p)k−1 is normalized i.e. null mean and identity variance. This aspect has been extensively explored in Deep Learning through batch, weight, layer and group normalization (Ioffe & Szegedy, 2015; Salimans & Kingma, 2016; Ba et al., 2016b; Wu & He, 2018). Normalizing the state variable corresponds to the following transformation:
s (p) k ←
(√ Var ( x
(p) k
))−1 ( x
(p) k − Ep
[ x
(p) k
]) . (25)
The Sylvester equation for theNk could be entirely circumvented, by using the proposed VarianceInvariance Theorem described hereafter.
Variance-Invariance Theorem: In addition to the hypotheses of the ESN Thereom (Section 2.2). If ∀k ∈ [1..n] : Var ( x
(p) k−1
) 0, then the randomized ESN is invariant to the variances
Var ( x
(p) k−1 ) which can be substituted with the identity. The randomized ESN system becomes: Nk = −E [ I + y (p) k y (p)T k ]−1 E [ y (p) k ( ∂`R ∂x (p) k−1 )T] , ηk = NkE [ x (p) k−1 ] . (26)
Proof: The demonstration builds on the invariance of the Newton method to affine transformations and is detailed step by step in the Appendix (A.5)

heading: 3.3 DISCUSSION
text: Expectations and estimators: We used the expectation E extensively in this paper to derive the theoretical ESN solution. However, estimators would replace expectation during implementation. The performance and the computational cost would be directly related to the estimators choice: rank one, moving average, or even batch estimators are all possible with different compromises.
Notes on convergence: Several convergence results for the Newton method exist (e.g Bottou & Le Cun, 2005; Byrd et al., 2011; Meng et al., 2020). The ESN convergence fits well within the general framework developed by Sunehag et al. (2009), under some mild and non-restrictive conditions.
To batch or not to batch: After randomization, ESN could be used in batch mode. However, there are no guarantees that such a solution solves the ”batch” Sifrian equations. Therefore, we focus only on using spot or historical data.
6A σ−field is more appropriate, but we use only a random space Ω for the sake of simplicity.

heading: 4 RESULTS
text: We validate the ESN method for FNNs, using a classification task on MNIST dataset (Lecun et al., 1998), and also an autoencoder training with the CURVES dataset (Hinton & Salakhutdinov, 2006). Leaky Relu and He initializer (He et al., 2016; Glorot & Bengio, 2010) are used for all the training. Further implementation details are available in the Appendix (A.6). We report the results hereafter:

heading: 5 CONCLUSION
text: This paper showed that it is possible to derive an exact expression of the stochastic second-order update for Deep Learning via a second-order Lagrangian which we called Sifrian, under some regularization constraints. We revisited backpropagation using Lagrange multipliers method and developed a new framework to better understand the mechanics of the Newton method at the stochastic level. This framework yielded unexpected insights on how to correct the non-convexity, the importance of the first layer, and how to get a closed-form solution. We further proceeded to randomize the solution, to better align the training with its stochastic nature and avoid overfitting. We tested the solution using simple training models; and despite the known limitations of the Newton method, the ESN (Exact Stochastic Newton) solution was able to perform as well as adaptive methods in a non-convex setting. ESN holds great potential but still requires further refinement for batch training, parameters selection and estimators choice. Generalization of this work to other type of networks is also a research venue which needs further investigation.
REPRODUCIBILITY STATEMENT
The results of this paper are reproducible. Commented Matlab code is included in the supplementary material. All experiments are run on a workstation with eight cores 2.6 GHz CPU.

heading: A APPENDIX
text: A.1 C2 SIFRIAN EQUATIONS
The inner product used in the Sifrian is either the inner product of vectors or the trace product for matrices, which are equivalent after vectorization:
〈A,B〉 = trace ( BTA ) = vec (B) T vec (A) ,
〈x,y〉 = xTy = trace ( xyT ) .
The derivation is carried through perturbation factorization for e.g.:
S (Wk + δWk)− S (Wk) ∼ 〈
∂S ∂Wk , δWk
〉 .
We carry the derivation and we report hereafter the full Sifrian equations for a C2 activation functions: (1) − ∑ p∇F (p) k ( γ (p) k x (p)T k−1 + 1k=2..nb (p) k ζ (p)T k−1 ) + ∑ p∇2F (p) k b (p) k [ Nkx (p) k−1 + ηk − 1k=2..nWkζ (p) k−1 ] x (p)T k−1 = ∑ p∇F (p) k b (p) k x (p)T k−1 , (2) ∑ p∇2F (p) k b (p) k [ Nkx (p) k−1 + ηk − 1k=2..nWkζ (p) k−1 ] − ∑ p∇F (p) k ( γ (p) k + b (p) k ) = 0, (3) γ (p) k − 1k=1..n−1W Tk+1∇F (p) k+1γ (p) k+1 + ∑ m ∂2`R ∂x (p) k ∂x (p) m ζ (p) m ) +1k=1..n−1W T k+1∇2F (p) k+1b (p) k+1 [ Nk+1x (p) k + ηk+1 − 1k=1..n−1Wk+1ζ (p) k ] +1k=1..n−1N T k+1∇F (p) k+1b (p) k+1 = 0,
(4) ζ (p) k − 1k=2..n∇F (p) k Wkζ (p) k−1 +∇F (p) k [ Nkx (p) k−1 + ηk ] = 0.
A.2 ESN RESOLUTION FOR C2 ACTIVATION FUNCTIONS
The stochastic case for C2 functions has the same type of simplifications. Comparing the first and second equations:
(1) − ∑
p∇F (p) k
( γ
(p) k x (p)T k−1 + 1k=2..nb (p) k ζ (p)T k−1 ) + ∑
p∇2F (p) k b (p) k [ Nkx (p) k−1 + ηk − 1k=2..nWkζ (p) k−1 ] x (p)T k−1
= ∑
p∇F (p) k b (p) k x (p)T k−1 ,
(2) ∑
p∇2F (p) k b (p) k [ Nkx (p) k−1 + ηk − 1k=2..nWkζ (p) k−1 ] − ∑ p∇F (p) k ( γ (p) k + b (p) k ) = 0,
This yields that 1k=1..n−1ζp(k) = 0. The fourth equation state that [ Nkx (p) k−1 + ηk ] = 0 is always null except for the last layer. Let’s focus on the last layer. Three equations allows to characterize it: (2) b (p) n = ( ∇F (p)n )−1 ∇2F (p)n b(p)n [ Nnx (p) n−1 + ηn ] − γ(p)n , (3) γ (p) n + ∂2`R ∂x (p) n ∂x (p) n ζ (p) n = 0,
(4) ζ (p) n +∇F (p)n [ Nnx (p) n−1 + ηn ] = 0.
If we assume that the C2 has a non null second derivative everywhere then, the previous three equations yield a specific value of [ Nnx (p) n−1 + ηn ] = 0 and also ζ(p)n :
Except the last layer we have: γ(p)k = −b (p) k , and a similar system to the piecewise affine case could be derived.
A.3 EXAMPLE: A SIMPLE REGULARIZATION
In this section we give a concrete simple example to further explain the ESN method. In the case of an l2 cost function and regularization defined by:
R(x) = λ
2 ∑ p∈D ∑ k=1..n−1 〈 x (p) k ,x (p) k 〉 .
we have the following solution: Nk = − ∇F ( a (p) k ) b (p) k x (p)T k−1∥∥∥∇F(a(p)k )b(p)k ∥∥∥2 2 ,
ηk = 1k=n
[ ∇F ( a (p) n )]−1 ( d(p) − x(p)n ) −Nkx(p)k−1.
• The previous solution is straightforward to compute; the first order backpropagation provides all the elements.
• It is co-linear with the gradient except for the bias of the last layer.
• The adjoints b(p)k ) are impacted by choice of the regularization parameter λ. • The minimal rank Newton solution is not a descent direction. Most of the biases correspond
to an ascent direction. • The squared norm division is problematic and will lead to divergent solution if the norm
of the backpropagated errors is very small i.e. ∥∥∥∇F (a(p)k )b(p)k ∥∥∥2
2 → 0. A sufficient
condition to avoid divergence consists in defining layer by layer regularization such as: λk−1 ∝ ∥∥∥∇F (a(p)k )b(p)k ∥∥∥. • It is possible to remove the regularization i.e. λ → 0, this yields a ”granular” solution
which states that the second order optimal solution consists in just adjusting the last layer bias with a slope corrected error: [ ∇F ( a (p) n )]−1 ( d(p) − x(p)n ) . The granular solution is
not suitable for neural network training, as it overfits every sample.
A.4 PROOF: PROPOSITION RANDOMIZED ESN
We start from the randomized least-squares criteria.
argmin{N ,η} ∫ Ω=D dΩ (p) ∥∥∥∥∥∥ ( ∂`R
∂x (p) k−1
)T + y
(p)T k Nk ∥∥∥∥∥∥ 2
2
+ ∥∥∥Nkx(p)k−1 + ηk∥∥∥2
2 .
The optimum is a critical point i.e. : ∀k ( ∂
∂Nk , ∂ ∂ηk
) = (0, 0) .
Each layer could be optimized separately which is desirable feature that reduce the dimensionality of the problem. We get the following derivatives: ∂ ∂Nk → ∫ dΩ ( y (p) k ( ∂`R ∂x (p) k−1 )T + y (p) k y (p)T k Nk +Nkx (p) k−1x (p)T k−1 + ηkx (p)T k−1 ) = 0,
∂ ∂ηk → ∫ dΩ ( Nkx (p) k−1 + ηk ) = 0
Integration over Ω yields expectations and two sets of equations: E [ y (p) k ( ∂`R ∂x (p) k−1 )T] + E [ y (p) k y (p)T k ] Nk +NkE [ x (p) k−1x (p)T k−1 ] + ηkE [ x (p)T k−1 ] = 0, NkE [ x (p) k−1 ] + ηk = 0. We simplify the system by right multiplying the second equation with E [ x
(p) k−1
] and substracting,we
get:  E [ y (p) k y (p)T k ] Nk +NkVar ( x (p) k−1 ) = −E [ y (p) k ( ∂`R ∂x (p) k−1 )T] , NkE [ x (p) k−1 ] + ηk = 0.
The above expression corresponds to a Sylvester equation. The last two equations are the main result of the proposition.
A.5 INVARIANCE OF ESN TO NORMALISATION:
In this section we demonstrate that the ESN is invariant to normalisation which would allow for a simpler resolution of the Sylvester Equation. Such a result is expected for a Newton method, we do the demonstration step by step. First let’s start by reformulating the forward problem with the normalized state variable s:
s (p) k =
(√ Var ( x
(p) k
))−1 ( x
(p) k − Ep
[ x
(p) k ]) = Σ−1k ( x (p) k −mk ) .
The original state variable is an affine transformation of the normalized one:
x (p) k = Σks (p) k + mk.
The matrix Σk is symmetric and positive. We assume further that it is definite i.e. invertible in this context. The forward model becomes:
Σks (p) k + mk − F ( Wk ( Σk−1s (p) k−1 + mk−1 ) + βk ) .
The square-root of the variance and the expectation could be absorbed into the weights/biases:
 W ′k = WkΣk−1,
β′k = βk +Wkmk−1.
With the above notations, the forward model becomes:
Σks (p) k + mk = F ( W ′ks (p) k−1 + β ′ k ) .
We use a new adjoint u. The Lagrangian becomes: L(x,W , β,b) = `R + ∑ p∈D ∑ k=1..n 〈 Σks (p) k + mk − F ( W ′ks (p) k−1 + β ′ k ) ,u (p) k 〉 .
The derivative w.r.t. s(p)k should be null, which yields the following backpropagation equation.
∂L ∂s
(p) k
= Σku (p) k − 1k=1..n−1W ′ k+1
T∇F ( a
(p) k+1
) u
(p) k+1 + Σk
∂`R
∂x (p) k
= 0.
Given the definition of W ′k+1. The adjoint equation is the same as the original, hence u = b, and the first order adjoint is unmodified. The gradient is slightly modified: G′k = ∂L ∂W ′k = − ∑ p∇F ( a (p) k ) u (p) k s (p)T k−1 ,
g′k = ∂L ∂β′k
= − ∑ p∇F ( a (p) k ) u (p) k .
Let’s move to the Sifrian: S(s,u,W ′,β′,G′,g′,γ′, ζ′,N ′,η′) = ∑
k,p
〈 Σks (p) k + mk − F ( W ′ks (p) k−1 + β ′ k ) ,γ′k (p) 〉
+ 〈 Σku (p) k − 1k=1..n−1W ′k+1T∇F ( a (p) k+1 ) u (p) k+1 + Σk ∂`R ∂x
(p) k
, ζ′k (p) 〉 + ∑
k 〈 G′k + ∑ p∇F ( a (p) k ) u (p) k s (p)T k−1 ,N ′ k 〉 + 〈 g′k + ∑ p∇F ( a (p) k ) u (p) k ,η ′ k 〉 .
The Sifrian equations are:
(1) ∑
p∇F (p) k ( γ′k (p)s (p)T k−1 + u (p) k s (p)T k−1 + 1k=2..nu (p) k ζ ′ k−1 (p)T ) = 0,
(2) ∑
p∇F (p) k ( γ′k (p) + u (p) k ) = 0,
(3) Σkγ ′ k (p) − 1k<nW ′k+1T∇F (p) k+1γ ′ k+1
(p) + ∑
m Σk ∂2`R
∂x (p) k ∂x (p) m
Σmζ ′ m (p) + 1k<nN ′ k+1 T∇F (p)k+1u (p) k+1 = 0,
(4) Σkζ ′ k (p) − 1k=2..n∇F (p)k W ′kζ (p) k−1 +∇F (p) k ( N ′ks (p) k−1 + η ′ k ) = 0.
We follow the same steps of the ESN Theorem demonstration, add the white layer and we remove the granular solution. We finally get the following stochastic system: ( ∂`R ∂x (p) k )T Σk + ( ∇F ( a (p) k+1 ) u (p) k+1 )T N ′k+1 = 0,
N ′ks (p) k−1 + η ′ k = 0.
The randomization yields, E [ y (p) k y (p)T k ] N ′k +N ′ kVar ( s (p) k−1 ) = −E [ y (p) k ( ∂`R ∂x (p) k−1 )T] Σk, N ′kE [ s (p) k−1 ] + η′k = 0
The s is normalized hence:  Var ( s (p) k−1 ) = I, E [ s
(p) k−1
] = 0.
The randomized, ESN with the normalized state variable has the following solution: N ′k = −E [ I + y (p) k y (p)T k ]−1 E [ y (p) k ( ∂`R ∂x (p) k−1 )T] Σk,
η′k = 0
We revert to the original variables from Equ. A.5, add convexity correction and we get the following solution:  Nk = −E [ I + y (p) k y (p)T k ]−1 E [ y (p) k ( ∂`R ∂x (p) k−1 )T] , ηk = NkE [ x (p) k−1
] A.6 IMPLEMENTATION DETAILS:
In this section we specify technical details of our implementation. First, we remove the ”granular” terms from our solution, since they do overfit each sample by construction.: Nk = −E [ I + y (p) k y (p)T k ]−1 E [ y (p) k ( ∂`R ∂x (p) k−1 )T] , ηk = NkE [ x (p) k−1
] ESN is solved for each layer independently, which means that the inverted matrices are a scale smaller than the empirical Fisher matrix for e.g.. To further illustrate this point, let’s consider, for e.g., the feedforward network with layers sizes [784, 256, 32, 10] (MNIST classification): the empirical Fisher matrix would have a size of 210298x210298, while the biggest matrix inversion for ESN would be only 256x256. In practice, the computation requires initially two ex-
pectations ( E [ y
(p) k
( ∂`R
∂x (p) k−1
)T] ,E [ x
(p) k−1
]) which could be estimated using exponentially de-
caying moving averages (Similar to ADAM). The product with the non-centered precision matrix E [ I + y
(p) k y (p)T k
]−1 , could be approximated in several ways. We list hereafter at least three meth-
ods:
• Direct method: the non-centered covariance matrix could be estimated using an exponentially decaying moving average and inverted directly given the small sizes.
• Recursive Least-Squares (RLS): A sliding window could be used to update recursively the estimation ofN(k) using Sherman-Morrison formula with an initialization of the precision matrix at I . Such a method is a classic of online learning and is similar to Kalman filtering.
• Rank one method: instead of involving the historical data, the simplest/crudest estimation of the covariance matrix is based on the spot/instantaneous information. The precision matrix could be derived using Sherman-Morrison formula (applied once):
( I + y
(p) k y (p)T k
)−1 = I − y(p)k y(p)Tk 1 + ∥∥∥y(p)k ∥∥∥2 2  .
Probability distribution/weights: It is common to use weights in the context of least-squares regression. For Deep Learning, we could increase the importance of samples that are not well assimilated into the network, i.e., high errors(or vice-versa). This is possible for e.g. by choosing
probabilities/weights proportional to errors’ norms: w(p) ∼ ∥∥∥x(p)n − d(p)∥∥∥2
2 . The ESN regression
equations become: Nk = −E [ w(p)I + w(p)y (p) k y (p)T k ]−1 E [ w(p)y (p) k ( ∂`R ∂x (p) k−1 )T] , ηk = Nk E [ w(p)x (p) k−1 ] E[w(p)]
Effects of the admissible regularization: We stated in this paper that the admissible regularisation R(x) reduces the norm of the state variable x. To further illustrate this statistical effect, we will use MNIST classification as an example, with regularized and un-regularized (normal) stochastic gradient descent method. For the sake of simplicity we consider the following form of regularisation:
R (x) = 1
2 ∑ p∈D ∑ k=1..n−1 λ 〈 x (p) k ,x (p) k 〉 .
The architecture of the FNN is the same from the experimental section i.e. {784, 256, 32, 10}. We select a learning rate of µ = 0.001 and various values of λ = { 0, 10−2, 10−3, 10−4, 10−5 } .

Review summaries:
This paper claims that it is possible to compute the Newton method's update exactly for deep neural networks (multi-layer perceptrons). The motivation is that Newton's method, as a second-order optimizer that includes loss curvature information, should improve upon first-order optimizers, such as gradient descent, that use only first derivatives (gradient) of the loss. Second order methods tipically converge in a smaller number of iterations, but each update is expensive to compute and that's why they are not widely used in practice. It must be noted that Newton's method is almost never mentioned in the context of Neural Networks (NN), because it is only guaranteed to converge for convex loss functions (which is clearly not the case for NN). However, the authors claim they have some tricks to fix that. 

This paper proposes a stochastic second order method to train neural network under some specific regularisation criterion. The method is based on the Sifrian, an extension of the Lagrangian that splits the definition of the gradient of each layer's parameter as different constraints with their own multiplier. Solving the best direction from the Sifrian is complicated in the general case but can be done when considering only one sample. This allows for a stochastic algorithm to train the neural network. Finally, some limited  numerical experiments are conducted.