\section{Introduction}
\label{sec:intro}

Machine learning generally amounts to solving an optimization problem where a loss function has to be minimized. As the problems tackled are getting more and more complex (nonlinear, nonconvex, etc.), fewer efficient algorithms exist, and the best recourse seems to rely on iterative schemes that exploit first-order or second-order derivatives of the loss function to get successive improvements and converge towards a local minimum. This explains why variants of gradient descent are becoming increasingly ubiquitous in machine learning and have been made widely available in the main deep learning libraries, being the tool of choice to optimize deep neural networks. Other types of local algorithms exist when no derivatives are known \citep{sigaud2018policy}, but in this paper we assume that some derivatives are available and only consider first order gradient-based or second order Hessian-based methods.

Among these methods, vanilla gradient descent strongly benefits from its computational efficiency as it simply computes partial derivatives at each step of an iterative process. Though it is widely used, it is limited for two main reasons: it depends on arbitrary parameterizations and may diverge or converge very slowly if the step size is not properly tuned.
To address these issues, several lines of improvement exist. Here, we focus on two of them. On the one hand, 
first-order methods such as the natural gradient introduce particular metrics to restrict gradient steps and make them independent from parametrization choices \citep{amari98natural}. On the other hand, second-order methods use the Hessian matrix of the loss or its approximations to take into account its local curvature.

Both types of approaches enhance the vanilla gradient descent update, multiplying it by the inverse of a large matrix (of size $d^2$, where $d$ is the dimensionality of the parameter space). 
We propose a simple framework that unifies these first-order or second-order improvements of the gradient descent, 
and use it to study precisely the similarities and differences between the 6 aforementioned methods.
This general framework uses a first-order approximation of the loss and constrains the step with a quadratic norm. Therefore, 
each modification $\mymath{\delta \vth}$ of the vector of parameters $\mymath{\vec \theta}$ is computed via an optimization problem of the following form:
\begin{equation}
  \label{eq:framework}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T M(\mymath{\vec \theta}) \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
where $\nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})$ is the gradient of the loss $L(\mymath{\vec \theta})$, and $M(\mymath{\vec \theta})$ a symmetric positive-definite matrix.
The 6 methods differ by the matrix $M(\mymath{\vec \theta})$, which has an effect not only on the size of the steps, but also on the direction of the steps,
as illustrated in Figure~\ref{fig:metric_effect}.
\begin{figure}[!ht]
	\centering
	\includegraphics[width=0.67\textwidth]{metric_effect.pdf}
    \caption{Different metrics affect both the gradient step size and direction.
    Here, a different $\mymath{\delta \vth}$ is obtained with $M=I$ or $M$ an arbitrary symmetric positive-definite matrix.}
	\label{fig:metric_effect}
\end{figure}

The solution of the minimization problem~\eqref{eq:framework} has the following form (see Appendix~\ref{ap:solution}): 
$$\mymath{\delta \vth} = -\alpha M(\mymath{\vec \theta})^{-1} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta}).$$

In Section~\ref{sec:first}, we show how the vanilla gradient descent, the classical Gauss-Newton method and the 
natural gradient descent method fit into the proposed framework. It can be noted that these 3 approaches constrain the steps in a way that 
is independent from the loss function.
In Section~\ref{sec:second}, we consider approaches that depend on the loss, namely the gradient covariance matrix approach, Newton's method 
and the generalized Gauss-Newton method, and show that they also fit into the framework.
Table~\ref{summary} summarizes the different values of $M(\mymath{\vec \theta})$ for all 6 approaches. 

\begin{table}
\begin{center}
\begin{tabular}{c|l}
 $M(\mymath{\vec \theta})$ & Corresponding approach\\
  \hline
  & \\
 $I$ & vanilla gradient descent \\[2pt]
 ${{\rm I\!E}}{}_{\mymath{\vec s}}{\left[J_{\mymath{\vec x}}(\mymath{\vec \theta})^TJ_{\mymath{\vec x}}(\mymath{\vec \theta})\right]} + \lambda I$ & classical Gauss-Newton \\[4pt]
 ${{\rm I\!E}}{}_\mymath{\vec s}\left[ \mymath{\nabla_{\vth}} \log( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) ) \mymath{\nabla_{\vth}}\log( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) )^T \right] + \lambda I$ & natural gradient (with empirical Fisher matrix) \\[4pt]
 ${{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right] + \lambda I$ & gradient covariance matrix \\[4pt]
 $H(\mymath{\vec \theta}) + \lambda I$ & Newton's method\\[2pt]
 ${{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_{\mymath{\vec y}}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right] + \lambda I$ & generalized Gauss-Newton\\[2pt]
 \hline
\end{tabular}
\end{center}
 \caption{\label{summary} The matrices $M(\mymath{\vec \theta})$ associated to 6 popular variants of the gradient descent, when interpreted
 as instances of the optimization problem~\eqref{eq:framework}.
 See Section~\ref{sec:not} for the definitions of the notations.}
\end{table}

Providing a unifying view of several first-order and second-order variants of the gradient descent, 
the framework presented in this paper
makes the connections between different approaches more obvious,
and can hopefully give new insights on these connections and help clarifying
some of the literature. Finally, we believe that it can facilitate the selection between these 
methods given a specific problem.

\section{Problem statement and notations}
\label{sec:not}

\begin{tabular}{c|l}
 Notation & Description\\
 \hline
  & \\
 $\mymath{\vec s} = (\mymath{\vec x}, \mymath{\vec y})$ & a sample, with an input variable $\mymath{\vec x}$, and an output variable $\mymath{\vec y}$\\[2pt]

 $L(\mymath{\vec \theta})$ & the scalar loss to minimize, $\mymath{\vec \theta}$ being the vector of parameters\\[2pt]
 $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ & p.d.f. estimating the conditional distribution of the output variable\\[2pt]
 $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ & a finite-dimensional representation of the distribution $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$\\[2pt]
  $J_{\mymath{\vec x}}(\mymath{\vec \theta})$ & Jacobian of the function $\mymath{\vec \theta} \mapsto \mymath{\vec h_{\vth}}(\mymath{\vec x})$\\[2pt]
 ${{\rm I\!E}}{}_{\mymath{\vec s}}[\cdot]$ & expectation over the samples\\[2pt]
 $l_\mymath{\vec \theta}(\mymath{\vec s}) = l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))$ & the atomic loss of which $L$ is the average over the samples: $L(\mymath{\vec \theta}) = {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ l_\mymath{\vec \theta}(\mymath{\vec s}) \right]}$\\[2pt]
 $\mymath{\delta \vth}$ & small update of $\mymath{\vec \theta}$ computed at every iteration \\[2pt]
 $(\cdot)^T$ & transpose operator \\[2pt]
 $\| \vec v \|$ & Euclidean norm of the vector $\vec v$: $\| \vec v \| = \sqrt{\vec v^T \vec v}$ \\[2pt]
 ${{\rm I\!E}}{}_{a \sim p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})}[f(a)]$ & expected value of $f(a)$, when $a$ follows the distribution  $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$\\[2pt]
 $\nabla_\mymath{\vec \theta} f(\mymath{\vec \theta})$ & gradient w.r.t. $\mymath{\vec \theta}$: $\nabla_\mymath{\vec \theta} f(\mymath{\vec \theta}) = \frac{\partial f(\mymath{\vec \theta})}{\partial \mymath{\vec \theta}}$ \\[2pt]
  $M_{CGN}(\mymath{\vec \theta})$ & Classical Gauss-Newton matrix: $
M_{CGN}(\mymath{\vec \theta}) = {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[J_{\mymath{\vec x}}(\mymath{\vec \theta})^TJ_{\mymath{\vec x}}(\mymath{\vec \theta})\right]}
$\\[2pt]
 $\mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta})$ & Fisher information matrix of $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$\\[2pt]
 $F(\mymath{\vec \theta})$ & empirical Fisher matrix: $ F(\mymath{\vec \theta}) = {{\rm I\!E}}{}_\mymath{\vec s}\left[ \mymath{\nabla_{\vth}} \log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)} \mymath{\nabla_{\vth}}\log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)}^T \right]$\\[2pt] 
 $KL\big( p_1, p_2 \big)$ & Kullback-Leibler divergence between the probability distributions $p_1$ and $p_2$\\[2pt]
 $H(\mymath{\vec \theta})$ & Hessian of the loss $L$, defined by $\left[H(\mymath{\vec \theta})\right]_{i,j} = \frac{\partial^2 L}{\partial \mymath{\vec \theta}_i \partial \mymath{\vec \theta}_j}(\mymath{\vec \theta})$\\[2pt]
 $\mathcal{H}_{\mymath{\vec y}}(\mymath{\vec h_{\vth}}(\mymath{\vec x}))$ & Hessian of the function $\vec h \mapsto l(\mymath{\vec y}, \vec h)$ at $\vec h = \mymath{\vec h_{\vth}}(\mymath{\vec x})$\\[2pt]
 \hline
\end{tabular}

We consider a context of regression analysis in which, based on samples $\vec{s} = (\mymath{\vec x}, \mymath{\vec y})$, the objective is to estimate the conditional distribution of $\mymath{\vec y}$ given $\mymath{\vec x}$.
More formally, this conditional distribution is estimated by a parametrized 
probability density function (p.d.f.) $p_\mymath{\vec \theta}(\mymath{\vec y} | \mymath{\vec x})$, and the goal of the learning is to progressively optimize the vector $\mymath{\vec \theta}$ to improve the accuracy of this probability estimation. We furthermore assume that the 
p.d.f. $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$ can be represented by a finite-dimensional vector $\mymath{\vec h_{\vth}}(\mymath{\vec x})$. For instance, in many applications, $p_\mymath{\vec \theta}(\cdot | \mymath{\vec x})$ is
a multivariate Gaussian distribution, and in this case the vector $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ would typically contain the mean and covariance matrix components.

The accuracy of $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$ is measured via a loss function 
$L$ estimated over a dataset of samples $\mymath{\vec s}$. $L$ depends only on $\theta$ and we assume that it is 
expressed as the expected value (over the sample distribution) of an atomic loss, a loss \emph{per sample} $l_\mymath{\vec \theta}(\mymath{\vec s})$:
$$
L(\mymath{\vec \theta}) = {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ l_\mymath{\vec \theta}(s) \right]}.
$$
In practice, the expected value over the samples is estimated with an empirical mean over a batch
$\mathcal{B} = (\mymath{\vec s}_1, \mymath{\vec s}_2, \dots, \mymath{\vec s}_N)$, so the loss actually used can be written 
$\hat{L}_\mathcal{B}(\mymath{\vec \theta}) = \frac{1}{N}\sum_i{ l_\mymath{\vec \theta}(\mymath{\vec s}_i) }$, the gradient of 
which being directly expressible from the gradient of the atomic loss (w.r.t. $\mymath{\vec \theta}$). In the remainder of the paper, 
we keep expressions based on the expected value $ {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ \cdot \right]}$, knowing that
at every iteration it is replaced by an empirical mean over a (new) batch of samples.

The dependency of the atomic loss to $\mymath{\vec \theta}$ is via $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$, so we can also express it as a function of the finite-dimensional representation $\mymath{\vec h_{\vth}}(\mymath{\vec x})$:
$$
l_\mymath{\vec \theta}(\mymath{\vec s}) = l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x})).
$$

Remark: choosing a concrete context of regression analysis helps 
us to simplify the notations, and give examples, but the results obtained are not specific to this setting, as the 6 gradient descent variants considered can also be useful for other types of learning tasks, leading to similar expressions that can be brought back to our general framework.  

\section{Vanilla, classical Gauss-Newton and natural gradient descent}
\label{sec:first}

All variants of the gradient descent are iterative numerical optimization methods: they start 
with a random initial \mymath{\vec \theta} and attempt to decrease the value of $L(\mymath{\vec \theta})$ over iterations by adding a 
small increment vector \mymath{\delta \vth} to \mymath{\vec \theta} at each step. The core of all these algorithms 
is to determine the direction and magnitude of \mymath{\delta \vth}.

\subsection{Vanilla gradient descent}
\label{sec:vanilla}
The so-called ``vanilla'' gradient descent is a first-order method that relies on a first-order Taylor approximation of the loss function $L$:
\begin{equation}
\label{eq:loss_taylor_first_order}
L(\mymath{\vec \theta} + \mymath{\delta \vth}) \simeq L(\mymath{\vec \theta}) + \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}.
\end{equation}

At each iteration, the objective is the minimization of $\nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}$ with the variable $\mymath{\delta \vth}$.
If the gradient is non-zero, the value of this term is unbounded below: it suffices for instance
to set $\mymath{\delta \vth} = -\alpha \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})$ with $\alpha$ arbitrarily large.
As a result, constraints are needed to 
avoid making excessively large steps. 
In vanilla approaches, the Euclidean metric ($\|\mymath{\delta \vth}\| = \sqrt{\mymath{\delta \vth}^T \mymath{\delta \vth}}$) is used to bound the increments $\mymath{\delta \vth}$. 
The optimization problem solved at every step of the scheme is: 
\begin{equation}
  \label{eq:vanilla}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth} \\ 
    \mymath{\delta \vth}^T \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
\noindent where $\epsilon$ is a user-defined upper bound. It is 
the most trivial instance of the general framework~(\ref{eq:framework}). As 
shown in Appendix~\ref{ap:solution}, the solution
of this problem is $\mymath{\delta \vth} = -\alpha \mymath{\nabla_{\vth}} L(\mymath{\vec \theta})$, with $\alpha = \frac{\epsilon}{\|\nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})\|}$.

To set the size of the step, instead of tuning $\epsilon$, the most common approach is 
to use the expression $-\alpha \mymath{\nabla_{\vth}} L(\mymath{\vec \theta})$ and directly tune $\alpha$, which is called the \emph{learning rate}.
An interesting property with this approach is that, as $\mymath{\vec \theta}$ gets closer to an optimum, the norm of the gradient $\|\nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})\|$ decreases, so 
the $\epsilon$ corresponding to the fixed $\alpha$ decreases as well.
This means that the steps tend to become smaller and smaller, which is a 
necessary property to make asymptotic convergence possible.

\subsection{Classical Gauss-Newton}
\label{sec:cgn}

As mentioned in Section~\ref{sec:not}, the atomic loss function $l_\mymath{\vec \theta}(\mymath{\vec s}) = l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))$ depends indirectly on the parameters $\mymath{\vec \theta}$ via the vector $\mymath{\vec h_{\vth}}(\mymath{\vec x})$, which is a finite-dimensional representation of the p.d.f. $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$. So it may be more meaningful to measure the modifications arising from an update $\mymath{\delta \vth}$ by looking at the change on $\mymath{\vec h_{\vth}}(\mymath{\vec x})$, not simply on $\mymath{\delta \vth}$ as with the vanilla gradient descent approach.
The constraint $\mymath{\delta \vth}^T \mymath{\delta \vth} \leq \epsilon^2$ acts as if all components of $\mymath{\delta \vth}$ 
had the same importance, which is not necessarily the case. Some components of $\mymath{\vec \theta}$ might have much smaller effects on $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ 
than others, and this will not be taken into account with the vanilla gradient descent method, which typically performs badly with unbalanced parametrizations.
Measuring and bounding the change on the vector $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ makes the updates independent from the way $\mymath{\vec h_{\vth}}$ is parametrized. To do this, a natural choice is to bound the expected squared Euclidean distance between $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ and $\mymath{\vec h_{\vth + \delth}}(\mymath{\vec x})$:
$$
{{\rm I\!E}}{}_{\mymath{\vec s}}{\left[\|\mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}) - \mymath{\vec h_{\vth}}(\mymath{\vec x})\|^2\right]}\leq \epsilon^2.
$$
Using again a first-order approximation, we have 
$\mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}) - \mymath{\vec h_{\vth}}(\mymath{\vec x}) \simeq J_{\mymath{\vec x}}(\mymath{\vec \theta})\mymath{\delta \vth}$, 
where $J_{\mymath{\vec x}}(\mymath{\vec \theta})$ is the Jacobian of the function $\mymath{\vec \theta} \mapsto \mymath{\vec h_{\vth}}(\mymath{\vec x})$.
The constraint can be rewritten:
$$
{{\rm I\!E}}{}_{\mymath{\vec s}}{\left[\|J_{\mymath{\vec x}}(\mymath{\vec \theta})\mymath{\delta \vth}\|^2\right]} = \mymath{\delta \vth}^T{{\rm I\!E}}{}_{\mymath{\vec s}}{\left[J_{\mymath{\vec x}}(\mymath{\vec \theta})^TJ_{\mymath{\vec x}}(\mymath{\vec \theta})\right]}\mymath{\delta \vth} \leq \epsilon^2,
$$
resulting in the optimization problem:
\begin{equation}
  \label{eq:jacobian}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[J_{\mymath{\vec x}}(\mymath{\vec \theta})^TJ_{\mymath{\vec x}}(\mymath{\vec \theta})\right]} \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
which fits into the general framework~\eqref{eq:framework} if the matrix $M_{CGN}(\mymath{\vec \theta}) = {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[J_{\mymath{\vec x}}(\mymath{\vec \theta})^TJ_{\mymath{\vec x}}(\mymath{\vec \theta})\right]}$ is symmetric positive-definite. 

\paragraph{Damping.}
The structure of the matrix $M_{CGN}(\mymath{\vec \theta})$ makes it symmetric and 
positive semi-definite, but not necessarily definite-positive. To ensure the definite-positiveness, 
a regularization or damping term $\lambda I$ can be added, resulting
in the constraint $\mymath{\delta \vth}^T \big(M_{CGN}(\mymath{\vec \theta}) + \lambda I \big) \mymath{\delta \vth} \leq \epsilon^2$, which can be
rewritten: 
$$
\mymath{\delta \vth}^T M_{CGN}(\mymath{\vec \theta}) \mymath{\delta \vth} + \lambda \mymath{\delta \vth}^T \mymath{\delta \vth} \leq \epsilon^2.
$$
We see that this kind of damping, often called Tikhonov damping 
\citep{martens2012training}, regularizes the constraint with
a term proportional to the squared Euclidean norm of $\mymath{\delta \vth}$.
It must be noted that with a regularization term, the constraint 
is not independent to the parametrization in $\mymath{\vec \theta}$ anymore. And if a large value of $\lambda$ is chosen (which also usually requires
increasing $\epsilon$), the method becomes similar 
to the vanilla gradient descent.

\paragraph{A more common definition of the classical Gauss-Newton method concerns atomic losses expressed as squared errors.} 
We assume that the atomic loss $l_\mymath{\vec \theta}(\mymath{\vec s})$ is defined as follows: 
$$l_\mymath{\vec \theta}(\mymath{\vec s}) = \frac{1}{2}\|\mymath{\Delta_{\vth}}(\mymath{\vec s})\|^2,$$
where $\Delta_\mymath{\vec \theta}(\mymath{\vec s})$ is a vector-valued function.
Functions of the form $\Delta_\mymath{\vec \theta}(\mymath{\vec s}) = \mymath{\vec y} - f_\mymath{\vec \theta}(\mymath{\vec x})$ are typical examples in the context of regression.

Denoting by $\mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})$ the Jacobian of $\mymath{\vec \theta} \mapsto \mymath{\Delta_{\vth}}(\mymath{\vec s})$: 
$\mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta}) = \left(\frac{\partial \Delta_{\mymath{\vec \theta}_i}(\mymath{\vec s})}{\partial \mymath{\vec \theta}_j}\right)_{i,j}$, we have:
\begin{equation*}
\begin{split}
l_{\mymath{\vec \theta} + \delta\mymath{\vec \theta}}(\mymath{\vec s}) & = \frac{1}{2}
\left( \Delta_\mymath{\vec \theta}(\mymath{\vec s}) + \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta}) \mymath{\delta \vth} + O(\mymath{\delta \vth}^2) \right)^T  \left( \Delta_\mymath{\vec \theta}(\mymath{\vec s}) + \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})\mymath{\delta \vth} + O(\mymath{\delta \vth}^2) \right)\\
& = l_\mymath{\vec \theta}(\mymath{\vec s}) + \Delta_\mymath{\vec \theta}(\mymath{\vec s})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})\mymath{\delta \vth} + \frac{1}{2}\mymath{\delta \vth}^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})\mymath{\delta \vth} + O(\mymath{\delta \vth}^3).
\end{split}
\end{equation*}
$\mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T\Delta_\mymath{\vec \theta}(\mymath{\vec s})$ is the gradient of the loss $l$ in $\mymath{\vec \theta}$, so the equality can be rewritten:
$$
l_{\mymath{\vec \theta} + \delta\mymath{\vec \theta}}(\mymath{\vec s}) = l_\mymath{\vec \theta}(\mymath{\vec s}) + \mymath{\nabla_{\vth}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \mymath{\delta \vth} + \frac{1}{2}\mymath{\delta \vth}^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})\mymath{\delta \vth} + O(\mymath{\delta \vth}^3).
$$
By averaging over the samples, we get:
$$
L(\mymath{\vec \theta} + \mymath{\delta \vth}) = L(\mymath{\vec \theta}) +  \mymath{\nabla_{\vth}} L(\mymath{\vec \theta})^T\mymath{\delta \vth} + \frac{1}{2}\mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta}) \right]} \mymath{\delta \vth} + O(\mymath{\delta \vth}^3).
$$
The classical Gauss-Newton method is often motivated by the minimization of this second-order Taylor approximation (see \cite{bottou2018optimization}). 
Assuming that ${{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta}) \right]}$ is positive-definite, 
as shown in Appendix~\ref{ap:quad} the minimum is reached with
$$
\mymath{\delta \vth} = -{{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta})^T \mathcal{J}^\Delta_{\mymath{\vec s}}(\mymath{\vec \theta}) \right]}^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}).
$$
As in the update obtained with the optimization problem~(\ref{eq:jacobian}), the matrix with a structure of type Jacobian transpose-times-Jacobian is characteristic of the classical Gauss-Newton approach. To ensure positive-definiteness, damping can be added in the exact same way. 
The derivation that lead to~(\ref{eq:jacobian}) shows that this kind of update does not only make sense with a squared error-type of loss, so in some sense it is 
a generalization of the context in which a classical Gauss-Newton approach may be useful. 
If the dependency of the loss to $\mymath{\vec \theta}$ is naturally expressed via a finite-dimensional vector $\vec v(\mymath{\vec \theta})$ (e.g. $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ or $\mymath{\Delta_{\vth}}(\mymath{\vec s})$ in the above cases), then 
measuring the quantity $\|\vec v(\mymath{\vec \theta} + \mymath{\delta \vth}) - \vec v(\mymath{\vec \theta})\|$ to evaluate the magnitude of the modifications induced by $\mymath{\delta \vth}$ is likely to be more meaningful than using the vanilla approach (i.e. simply measuring $\|(\mymath{\vec \theta} + \mymath{\delta \vth}) - \mymath{\vec \theta}\| = \|\mymath{\delta \vth} \| $).

\paragraph{Learning rate.}
As shown in Appendix~\ref{ap:solution}, the general framework~\eqref{eq:framework}
has a unique solution $\mymath{\delta \vth}=-\alpha M(\mymath{\vec \theta})^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta})$, 
with $\alpha = \frac{\epsilon}{\sqrt{\mymath{\nabla_{\vth}} L(\mymath{\vec \theta})^T M(\mymath{\vec \theta})^{-1}\mymath{\nabla_{\vth}} L(\mymath{\vec \theta})}}$. The classical Gauss-Newton approach corresponds to $M(\mymath{\vec \theta}) = M_{CGN}(\mymath{\vec \theta}) + \lambda I$, or $M(\mymath{\vec \theta}) = M_{CGN}(\mymath{\vec \theta})$ if we ignore the damping.
With the approach based on the second-order approximation of the loss expressed as a squared error, the resulting update has a form $\mymath{\delta \vth}=- M(\mymath{\vec \theta})^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta})$,
which is similar to the above expression except that $\alpha = 1$. However, this theoretical difference in $\alpha$ (referred to as the learning rate in Section~\ref{sec:vanilla}) is not very significant in practice since its value is usually redefined separately. It is very common to use heuristics to set $\alpha$ to smaller values so as to increase the stability of the iterative method. Indeed,
in the proposed framework, if $\epsilon$ is constant, the updates are not getting smaller and smaller, which means that no convergence is possible.
Another example of motivation for the redefinition of $\alpha$ is that when $M(\mymath{\vec \theta})$ is a very large matrix, $M(\mymath{\vec \theta})^{-1}$ is often estimated via drastic approximations. In that case, it can be preferable to only keep the update direction of the solution of~\eqref{eq:framework}, and then perform a line search to find a value of $\alpha$ for which it is reverified that the corresponding step size is reasonable.
This line search is an important component of the popular reinforcement learning
algorithm TRPO \citep{schulman2015trust}.

\subsection{Natural gradient}
\label{sec:nat}

To go one step further in terms of independence to the parametrization, it is possible to measure directly the
change from $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$ to $p_{\mymath{\vec \theta} + \mymath{\delta \vth}}( \cdot | \mymath{\vec x})$ with a metric on probability density functions.
This way, the updates do not even depend on the choice of
finite-dimensional representation via $\mymath{\vec h_{\vth}}$.
Amari \citep{amari1997neural, amari98natural} proposed and 
 popularized the notion of natural gradient, which is based on a matrix called the 
Fisher information matrix, defined for the p.d.f $p_\mymath{\vec \theta}( \cdot | \mymath{\vec x})$ by: 
$$
\mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) = {{\rm I\!E}}{}_{a\sim p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})}{\left[ \mymath{\nabla_{\vth}} \log{\left( p_\mymath{\vec \theta}(a|\mymath{\vec x}) \right)} \mymath{\nabla_{\vth}}\log{\left( p_\mymath{\vec \theta}(a|\mymath{\vec x}) \right)}^T \right]}.
$$
It can be used to measure a ``distance'' $d\ell$ between two infinitesimally close 
probability distributions $ p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ and $p_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\cdot|\mymath{\vec x})$ as follows:
$$
d\ell^2(p_\mymath{\vec \theta}(\cdot|\mymath{\vec x}), p_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\cdot|\mymath{\vec x})) = \mymath{\delta \vth}^T \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \mymath{\delta \vth}.
$$
Averaging over the samples, we extrapolate a measure of distance between $\mymath{\vec \theta}$ and $\mymath{\vec \theta} + \mymath{\delta \vth}$:
$$
DL^2(\mymath{\vec \theta}, \mymath{\vec \theta} + \mymath{\delta \vth}) = \mymath{\delta \vth}^T {{\rm I\!E}}{}_\mymath{\vec s}\left[ \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \right] \mymath{\delta \vth},
$$
where ${{\rm I\!E}}{}_\mymath{\vec s}\left[ \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \right] = {{\rm I\!E}}{}_\mymath{\vec s}\left[ {{\rm I\!E}}{}_{a\sim p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})}{\left[ \mymath{\nabla_{\vth}} \log{\left( p_\mymath{\vec \theta}(a|\mymath{\vec x}) \right)} \mymath{\nabla_{\vth}}\log{\left( p_\mymath{\vec \theta}(a|\mymath{\vec x}) \right)}^T \right]} \right]$ is 
the averaged Fisher information matrix.
It is common to approximate ${{\rm I\!E}}{}_\mymath{\vec s}\left[ {{\rm I\!E}}{}_{a \sim p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})}[\cdot] \right]$ with the 
empirical mean over the samples,
which reduces the above expression to
$$
DL^2(\mymath{\vec \theta}, \mymath{\vec \theta} + \mymath{\delta \vth}) \approx \mymath{\delta \vth}^T {{\rm I\!E}}{}_\mymath{\vec s}\left[ \mymath{\nabla_{\vth}} \log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)} \mymath{\nabla_{\vth}}\log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)}^T \right] \mymath{\delta \vth}.
$$
$ {{\rm I\!E}}{}_\mymath{\vec s}\left[ \mymath{\nabla_{\vth}} \log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)} \mymath{\nabla_{\vth}}\log{\left( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) \right)}^T \right]$ is called the \emph{empirical}
Fisher matrix \citep{martens2014new}. We denote it by $F(\mymath{\vec \theta})$.
Putting an upper bound on $\mymath{\delta \vth}^T F(\mymath{\vec \theta}) \mymath{\delta \vth}$ results in the following optimization problem:
\begin{equation}
  \label{eq:natural}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T F(\mymath{\vec \theta}) \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
which yields \emph{natural} gradient steps of the form
$$
\mymath{\delta \vth} = -\alpha F(\mymath{\vec \theta})^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}),
$$
provided that $F(\mymath{\vec \theta})$ is invertible. $F(\mymath{\vec \theta})$ is always positive semi-definite. 
Therefore, as in Section~\ref{sec:cgn} with the classical Gauss-Newton approach, a damping term
can be added to ensure invertibility (but again, by doing so the independence to the parametrization is lost).
The Fisher information matrix is in some sense uniquely defined by the property of
invariance to reparametrization of the metric it induces \citep{cencov1982}, and it can be obtained 
from many different derivations. But a particularly interesting fact is that 
$d\ell^2(p_\mymath{\vec \theta}(\cdot|\mymath{\vec x}), p_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\cdot|\mymath{\vec x}))$ corresponds 
to the second-order approximation of the Kullback-Leibler divergence $KL(p_\mymath{\vec \theta}(\cdot|\mymath{\vec x}),p_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\cdot|\mymath{\vec x}))$
\citep{kullback97, akimoto2013objective}.
Hence, the terms $\mymath{\delta \vth}^T \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \mymath{\delta \vth}$ and $\mymath{\delta \vth}^T F(\mymath{\vec \theta}) \mymath{\delta \vth}$ share some of the properties of the Kullback-Leibler divergence. For instance, when the variance of the probability distribution $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ decreases, 
the same parameter modification $\mymath{\delta \vth}$ tends to result in increasingly large 
measures $\mymath{\delta \vth}^T \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \mymath{\delta \vth}$ (see Figure~\ref{fig:kl}). 

\begin{figure}[H]
\includegraphics[width=\textwidth]{kl_new.pdf}
    \caption{The same parameter change (here, a constant shift of the mean to the right) yields a larger Kullback-Leibler divergence when the variance is small.}
	\label{fig:kl}
	\end{figure}
    
Consequently, if the bound $\epsilon^2$ of Equation~(\ref{eq:natural}) is kept constant, the possible modifications of \mymath{\vec \theta} become somehow smaller 
when the variance of the parametrized distribution decreases. Thus the natural gradient iterations slow down when the variance becomes small, which is a desirable 
property when keeping some amount of variability is important.
Typically, in the context of reinforcement learning, this variability can be related to exploration, and it should not vanish early.
This is one of the reasons why several reinforcement learning algorithms benefit from the use of natural gradient 
steps \citep{PetersNAC,schulman2015trust,wu2017scalable}.

\paragraph{Relation between natural gradient and classical Gauss-Newton approaches.}
Let us consider a very simple case where 
$p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ is a multivariate normal distribution with fixed covariance matrix $\Sigma = \beta^2 I$.
The only variable parameter on which the distribution $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ depends is its mean $\vec \mu_\mymath{\vec \theta}$, so we can use it as
representation of the distribution itself and write
$$
\mymath{\vec h_{\vth}}(\mymath{\vec x}) = \vec \mu_\mymath{\vec \theta}.
$$
It can be shown that the Kullback-Leibler divergence between two normal distributions
of equal variance and different means is proportional 
to the squared Euclidean distance between the means.
More precisely, the Kullback-Leibler divergence between 
$p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$ and $p_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\cdot|\mymath{\vec x})$ is equal to 
$\frac{1}{2\beta^2}\| \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}) - \mymath{\vec h_{\vth}}(\mymath{\vec x})\|^2$.
For small values of $\mymath{\delta \vth}$, this expression is
approximately equal to the measure 
obtained with the true Fisher information
matrix: 
$$
\frac{1}{2\beta^2}\| \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}) - \mymath{\vec h_{\vth}}(\mymath{\vec x})\|^2 \approx \mymath{\delta \vth}^T \mathcal{I}_\mymath{\vec x}(\mymath{\vec \theta}) \mymath{\delta \vth}.
$$

Bounding the average over the samples of the right term is the motivation of the natural gradient descent method.
Besides, we have seen in Section~\ref{sec:cgn} that the
classical Gauss-Newton method can be considered as a way to bound 
${{\rm I\!E}}{}_{\mymath{\vec s}}[\| \mymath{\vec h_{\vth + \delth}}(\mymath{\vec s}) - \mymath{\vec h_{\vth}}(\mymath{\vec x})\|^2]$, which is equal to the average of the left term over the samples,
up to a multiplicative constant.
Hence, even though both methods introduce slightly different approximations,
we can conclude that in this context, 
the classical Gauss-Newton and natural gradient descent methods are very similar.
This property is used in \cite{pascanu2013revisiting} to perform a
natural gradient descent on deterministic neural networks, by interpreting
their outputs as the mean of a conditional Gaussian distribution with fixed variance.

\section{Gradient covariance matrix, Newton's method and generalized Gauss-Newton}
\label{sec:second}

The approaches seen in Section~\ref{sec:first} all fit the general framework~\eqref{eq:framework}:
\begin{equation*}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T M(\mymath{\vec \theta}) \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation*}
with matrices $M(\mymath{\vec \theta})$ that do not depend on the loss function. But since
the loss is typically based on quantities that are relevant for the task to achieve,
it can be a good idea to exploit it to constrain the steps.
We present 3 approaches that fit into the same framework but with matrices $M(\mymath{\vec \theta})$ 
that do depend on the loss, namely the gradient covariance matrix method, 
Newton's method, and the generalized Gauss-Newton method.

\subsection{Gradient covariance matrix}

The simplest way to use the loss to measure the
magnitude of a change due to parameter modifications
is to consider the expected squared difference between 
$l_\mymath{\vec \theta}(\mymath{\vec s})$ and $l_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\mymath{\vec s})$:
$$
{{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \big(l_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\mymath{\vec s}) - l_\mymath{\vec \theta}(\mymath{\vec s})\big)^2 \right].
$$
For a single sample \mymath{\vec s}, changing slightly the object $\mymath{\vec \theta}$ 
does not necessarily modify the atomic loss $l_\mymath{\vec \theta}(\mymath{\vec s})$, 
but in many cases it can be assumed that this loss becomes 
different for at least some of the samples, yielding a positive 
value for  
$
{{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \big(l_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\mymath{\vec s}) - l_\mymath{\vec \theta}(\mymath{\vec s})\big)^2 \right]
$ which quantifies in some sense the amount of change introduced by $\mymath{\delta \vth}$ with 
respect to the objective. It is often a meaningful measure as it usually depends
on the most relevant features for the task to achieve.
Let us replace $l_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\mymath{\vec s})$ by a first-order approximation: 
$$ l_{\mymath{\vec \theta}+\mymath{\delta \vth}}(\mymath{\vec s}) \simeq l_\mymath{\vec \theta}(\mymath{\vec s}) + \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \mymath{\delta \vth}.$$
The above expectation simplifies to
$$
{{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \big(\nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \mymath{\delta \vth}\big)^2 \right]
= \mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right] \mymath{\delta \vth}.
$$
${{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right]$ is called the 
gradient covariance matrix \citep{bottou2008tradeoffs}.
It can also be called the outer product metric~\citep{ollivier:hal-00857982}.
Putting a bound on $\mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right] \mymath{\delta \vth}$, 
the iterated optimization becomes: 
\begin{equation}
  \label{eq:covgrad}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth} \\ 
    \mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right] \mymath{\delta \vth}  \leq \epsilon^2.
  \end{array}
  \right.
\end{equation}

It results in updates of the form: 
$$
\mymath{\delta \vth} = -\alpha {{\rm I\!E}}{}_{\mymath{\vec s}}\left[\nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right]^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}).
$$
Again, a regularization term may be added to ensure the invertibility of the matrix.
\paragraph{Link with the natural gradient.}
Let us assume that the atomic loss on a sample $\mymath{\vec s} = (\mymath{\vec x}, \mymath{\vec y})$ is 
the negative log-likelihood (a very common choice):
$$
l_\mymath{\vec \theta}(\mymath{\vec s}) = -\log( p_\mymath{\vec \theta}(\mymath{\vec y}|\mymath{\vec x}) ).
$$
It follows that the empirical Fisher matrix, as defined in Section~\ref{sec:nat}, is equal 
to ${{\rm I\!E}}{}_{\mymath{\vec s}}\left[ \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s}) \nabla_{\mymath{\vec \theta}} l_\mymath{\vec \theta}(\mymath{\vec s})^T \right]$,
which is exactly the definition of the gradient covariance matrix. Therefore, in this 
case, the two approaches are identical. Several algorithms use this identity for the natural gradient
computation, e.g. \cite{george2018fast}.

\subsection{Newton's method}
\label{sec:Newton}

Let us consider now a second-order approximation of the loss: 
$$
L(\mymath{\vec \theta} + \mymath{\delta \vth}) \approx L(\mymath{\vec \theta}) + \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})^T \mymath{\delta \vth} + \frac{1}{2}\mymath{\delta \vth}^T H(\mymath{\vec \theta}) \mymath{\delta \vth} ,
$$
where $H(\mymath{\vec \theta})$ is the Hessian matrix: $\left[H(\mymath{\vec \theta})\right]_{i,j} = \frac{\partial^2 L}{\partial \mymath{\vec \theta}_i \partial \mymath{\vec \theta}_j}(\mymath{\vec \theta})$.
Although there are obvious counterexamples, one can argue that the first-order approximation, i.e. $L(\mymath{\vec \theta} + \mymath{\delta \vth}) \approx L(\mymath{\vec \theta}) + \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})^T \mymath{\delta \vth}$ 
(which is used as minimization objective for gradient descents), is most likely good
as long as the second-order term $\frac{1}{2}\mymath{\delta \vth}^T H(\mymath{\vec \theta}) \mymath{\delta \vth}$ remains small. Therefore, 
it makes sense to directly put an upper bound on this quantity to restrict $\mymath{\delta \vth}$, as follows:
$$
\mymath{\delta \vth}^T H(\mymath{\vec \theta}) \mymath{\delta \vth} \leq \epsilon^2.
$$
This constraint defines a trust region, i.e. a neighborhood of $\mymath{\vec \theta}$
in which the first-order approximation of $L(\mymath{\vec \theta} + \mymath{\delta \vth})$ is supposed to be reasonably accurate. The trust region is bounded and well defined if the matrix $H(\mymath{\vec \theta})$ is symmetric positive-definite. 
However, $H(\mymath{\vec \theta})$ is symmetric but not even necessarily positive semi-definite, unlike the matrices 
obtained with the previous 
approaches. Therefore the required damping to make it definite-positive may be larger than with other methods.
It leads to the following optimization problem solved at every iteration:
\begin{equation}
  \label{eq:hess}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth} \\ 
    \mymath{\delta \vth}^T (H(\mymath{\vec \theta}) + \lambda I) \mymath{\delta \vth}  \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
and to updates of the form:
$$
\mymath{\delta \vth} = -\alpha (H(\mymath{\vec \theta}) + \lambda I)^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}).
$$
\paragraph{The more usual derivation of Newton's method.} The same update 
direction is obtained by directly minimizing the damped second-order approximation:
$$
L(\mymath{\vec \theta}) + \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})^T \mymath{\delta \vth} + \frac{1}{2}\mymath{\delta \vth}^T (H(\mymath{\vec \theta}) + \lambda I) \mymath{\delta \vth}.
$$
When $(H(\mymath{\vec \theta}) + \lambda I)$ is symmetric positive-definite, as shown in Appendix~\ref{ap:quad} the minimum of this 
expression is obtained for: 
$$
\mymath{\delta \vth} = - (H(\mymath{\vec \theta}) + \lambda I)^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}).
$$

\subsection{Generalized Gauss-Newton}

$L(\mymath{\vec \theta})$ is equal to ${{\rm I\!E}}{}_{\mymath{\vec s}}{\left[ l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x})) \right]}$: it
does not depend directly on $\mymath{\vec \theta}$ but on the outputs of $\mymath{\vec h_{\vth}}$, which are vectors of finite dimension.
Posing $\delta \vec h = \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}) - \mymath{\vec h_{\vth}}(\mymath{\vec x})$, a second-order Taylor expansion 
of $l(\mymath{\vec y}, \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x}))$ can be written:
$$
l(\mymath{\vec y}, \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x})) = l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x})) + \frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T \delta \vec h
+ \frac{1}{2}\delta \vec h^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) \delta \vec h + O(\delta \vec h^3),
$$
where $\mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x}))$ is the Hessian matrix of the atomic loss $l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))$
with respect to variations of $\mymath{\vec h_{\vth}}(\mymath{\vec x})$, and $\frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}$
is the gradient of $l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))$ w.r.t. variations of $\mymath{\vec h_{\vth}}(\mymath{\vec x})$.
Using the equality $\delta \vec h = J_{\mymath{\vec x}}(\mymath{\vec \theta}) \mymath{\delta \vth} + O(\mymath{\delta \vth}^2)$ (with $J_{\mymath{\vec x}}(\mymath{\vec \theta})$ the Jacobian of the function $\mymath{\vec \theta} \mapsto \mymath{\vec h_{\vth}}(\mymath{\vec x})$):
\begin{equation*}
\begin{split}
l(\mymath{\vec y}, \mymath{\vec h_{\vth + \delth}}(\mymath{\vec x})) = & \ l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x})) + \frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T J_{\mymath{\vec x}}(\mymath{\vec \theta}) \mymath{\delta \vth}
+ \frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T O(\mymath{\delta \vth}^2) \\
& + \frac{1}{2}\mymath{\delta \vth}^T J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}( \mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \mymath{\delta \vth} + O(\mymath{\delta \vth}^3).
\end{split}
\end{equation*}
The generalized Gauss-Newton approach is an approximation that consists in dropping
the term $\frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T O(\mymath{\delta \vth}^2)$. 
Averaging over the samples yields:
$$
L(\mymath{\vec \theta} + \mymath{\delta \vth}) \approx L(\mymath{\vec \theta}) + {{\rm I\!E}}{}_{\mymath{\vec s}}\left[\frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T J_{\mymath{\vec x}}(\mymath{\vec \theta})\right]^T \mymath{\delta \vth} + 
\frac{1}{2}\mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right] \mymath{\delta \vth}.
$$
Noticing that
${{\rm I\!E}}{}_{\mymath{\vec s}}\left[\frac{\partial l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))}{\partial \vec h}^T J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right] = \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})$,
it results in the following approximation:
$$
L(\mymath{\vec \theta} + \mymath{\delta \vth}) \approx L(\mymath{\vec \theta}) + \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})^T \mymath{\delta \vth} + 
\frac{1}{2}\mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right] \mymath{\delta \vth}.
$$
As for Newton's method, the usual way to derive the generalized Gauss-Newton method 
is to directly minimize this expression (see \cite{martens2014new}), but we can
also put a bound on the quantity $\mymath{\delta \vth}^T {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right] \mymath{\delta \vth}$ to define a trust region for the validity of the first-order approximation (as in Section~\ref{sec:Newton}),
provided that ${{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right]$ is symmetric positive-definite.
If the loss $l(\mymath{\vec y}, \mymath{\vec h_{\vth}}(\mymath{\vec x}))$ is convex in $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ (which is often true), 
the matrix is at least positive semi-definite, so a small damping term suffices to make it positive-definite.
If a non-negligible portion of the matrices $J_{\mymath{\vec x}}(\mymath{\vec \theta})$ are full rank,
the damping term may be added to $\mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x}))$ rather than to the full matrix. 
See \cite{martens2012training} for an extensive discussion on different options for damping 
and their benefits and drawbacks. With the damping on the full matrix,
the optimization problem to solve at every iteration becomes:
\begin{equation}
  \label{eq:GGN}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_\mymath{\vec \theta} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T \left( {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right]  + \lambda I\right) \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation}
resulting in updates of the form: 
$$
\mymath{\delta \vth} = - \alpha \left( {{\rm I\!E}}{}_{\mymath{\vec s}}\left[ J_{\mymath{\vec x}}(\mymath{\vec \theta})^T \mathcal{H}_\mymath{\vec y}(\mymath{\vec h_{\vth}}(\mymath{\vec x})) J_{\mymath{\vec x}}(\mymath{\vec \theta}) \right]  + \lambda I\right)^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}).
$$

\section{Summary and conclusion}
\label{sec:summ}

In sections~\ref{sec:first} and \ref{sec:second} we motivated and derived 6 different ways to compute parameter updates, that can all be interpreted as 
solving an optimization problem of this type:
\begin{equation*}
  \left\{
  \begin{array}{l}
    \min_{\mymath{\delta \vth}} \nabla_{\mymath{\vec \theta}} L(\mymath{\vec \theta})^T \mymath{\delta \vth}\\ 
    \mymath{\delta \vth}^T M(\mymath{\vec \theta}) \mymath{\delta \vth} \leq \epsilon^2,
  \end{array}
  \right.
\end{equation*}
resulting in updates of the form:
$$
\mymath{\delta \vth} = - \alpha M(\mymath{\vec \theta})^{-1} \mymath{\nabla_{\vth}} L(\mymath{\vec \theta}),
$$

$M(\mymath{\vec \theta})$ being symmetric positive-definite. The quadratic term
of the inequality corresponds to a specific metric defined by $M(\mymath{\vec \theta})$
used to measure the magnitude of the modification induced by $\mymath{\delta \vth}$. 
To evaluate this magnitude, the focus can simply be on the norm of $\mymath{\delta \vth}$, 
or on the effect of $\mymath{\delta \vth}$ on the loss, or on the effect
of $\mymath{\delta \vth}$ on $\mymath{\vec h_{\vth}}(\mymath{\vec x})$ or on $p_\mymath{\vec \theta}(\cdot|\mymath{\vec x})$, resulting in various approaches,
with various definitions of $M(\mymath{\vec \theta})$. In a context of probabilistic regression, we gave 6 examples that 
correspond to popular variants of the gradient descent, summarized in Table~\ref{summary}. All methods except the natural gradient can be declined to deterministic cases.
Unifying several first-order or second-order variants of the gradient descent method
enabled us to reveal links between these different approaches,
and contexts in which some of them are equivalent.
The proposed framework gives a compact overview
of common variants of the gradient descent, and
hopefully can help
choosing adequately between them depending on the problem to solve.
Perhaps, it can also help designing new variants or combining
existing ones to obtain new desired features.

\section*{Acknowledgements}

This research was partially supported by the French National Research Agency (ANR), Project \mbox{ANR-18-CE33-0005} HUSKI.


