# Adapting Newton's Method to Neural Networks through a Summary of Higher-Order Derivatives

Pierre Wolinski

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Wolinski. Adapting Newton's Method to Neural Networks through a Summary of Higher-Order Derivatives. 2024. hal-04600082

HAL Id: hal-04600082

https://hal.science/hal-04600082

Preprint submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Adapting Newton's Method to Neural Networks through a Summary of Higher-Order Derivatives

Pierre Wolinski¹

¹Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, France
pierre.wolinski@universite-paris-saclay.fr

February 6, 2024

## Abstract

We consider a gradient-based optimization method applied to a function $\mathcal{L}$ of a vector of variables $\theta$, in the case where $\theta$ is represented as a tuple of tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$. This framework encompasses many common use-cases, such as training neural networks by gradient descent.

First, we propose a computationally inexpensive technique providing higher-order information on $\mathcal{L}$, especially about the interactions between the tensors $\mathbf{T}_s$, based on automatic differentiation and computational tricks. Second, we use this technique at order 2 to build a second-order optimization method which is suitable, among other things, for training deep neural networks of various architectures.

This second-order method leverages the partition structure of $\theta$ into tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$, in such a way that it requires neither the computation of the Hessian of $\mathcal{L}$ according to $\theta$, nor any approximation of it. The key part consists in computing a smaller matrix interpretable as a “Hessian according to the partition”, which can be computed exactly and efficiently. In contrast to many existing practical second-order methods used in neural networks, which perform a diagonal or block-diagonal approximation of the Hessian or its inverse, the method we propose does not neglect interactions between layers.

Our code is available on GitHub: https://github.com/p-wol/GroupedNewton.

## 1 Introduction

The appealing theoretical properties of Newton's method have led to numerous attempts to adapt it to neural network optimization. Therefore, the study of the Hessian of a loss according to many parameters has become an area of research in itself, leading to a large number of methods to approximate it accurately with a small computational cost.

**Newton's method applied to neural network optimization.** When it comes to neural networks, Newton's method suffers from several problems. Some are *technical*, such as building an accurate and computationally efficient method to estimate the Hessian. But some of them are *essential*, in the sense that they cannot be solved solely by a perfect knowledge of the full Hessian. For instance, several works (Sagun et al., 2018) have shown that many eigenvalues of the Hessian are close to zero when training neural networks, making it impossible to use Newton's method in practice, even if the Hessian is perfectly known.

Therefore, in this work, we do not aim to build a technique to estimate the Hessian accurately and efficiently. We propose to take a step aside and focus on two related goals. First, we aim to access efficiently to higher-order information in order to use it for optimization. To do so, we leverage the available computational tools, that is, automatic differ-
# Adapting Newton's Method to Neural Networks through a Summary of Higher-Order Derivatives

Pierre Wolinski¹

¹Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, France
pierre.wolinski@universite-paris-saclay.fr

February 6, 2024

## Abstract

We consider a gradient-based optimization method applied to a function $\mathcal{L}$ of a vector of variables $\theta$, in the case where $\theta$ is represented as a tuple of tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$. This framework encompasses many common use-cases, such as training neural networks by gradient descent.

First, we propose a computationally inexpensive technique providing higher-order information on $\mathcal{L}$, especially about the interactions between the tensors $\mathbf{T}_s$, based on automatic differentiation and computational tricks. Second, we use this technique at order 2 to build a second-order optimization method which is suitable, among other things, for training deep neural networks of various architectures.

This second-order method leverages the partition structure of $\theta$ into tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$, in such a way that it requires neither the computation of the Hessian of $\mathcal{L}$ according to $\theta$, nor any approximation of it. The key part consists in computing a smaller matrix interpretable as a “Hessian according to the partition”, which can be computed exactly and efficiently. In contrast to many existing practical second-order methods used in neural networks, which perform a diagonal or block-diagonal approximation of the Hessian or its inverse, the method we propose does not neglect interactions between layers.

Our code is available on GitHub: https://github.com/p-wol/GroupedNewton.

## 1 Introduction

The appealing theoretical properties of Newton's method have led to numerous attempts to adapt it to neural network optimization. Therefore, the study of the Hessian of a loss according to many parameters has become an area of research in itself, leading to a large number of methods to approximate it accurately with a small computational cost.

**Newton's method applied to neural network optimization.** When it comes to neural networks, Newton's method suffers from several problems. Some are *technical*, such as building an accurate and computationally efficient method to estimate the Hessian. But some of them are *essential*, in the sense that they cannot be solved solely by a perfect knowledge of the full Hessian. For instance, several works (Sagun et al., 2018) have shown that many eigenvalues of the Hessian are close to zero when training neural networks, making it impossible to use Newton's method in practice, even if the Hessian is perfectly known.

Therefore, in this work, we do not aim to build a technique to estimate the Hessian accurately and efficiently. We propose to take a step aside and focus on two related goals. First, we aim to access efficiently to higher-order information in order to use it for optimization. To do so, we leverage the available computational tools, that is, automatic differ-
entiation and partition of the set of parameters into
tensors stored on a GPU. Second, we aim to build
a second-order optimization method with properties
similar to the properties of Newton's method (but
inevitably weaker). Notably, we want a method in-
variant by layer-wise affine reparameterizations of the
model, providing a Hessian-inspired matrix showing
the interactions between each pair of tensors, along
with layer-wise step sizes. Naturally, the computa-
tional cost should remain reasonable, and the prob-
lem of close-to-zero eigenvalues of Hessian has to be
dealt with.

**First contribution: extracting higher-order information.** Formally, we study a loss $\mathcal{L}$ to minimize according to a vector of parameters $\theta \in \mathbb{R}^p$, which can be represented as a tuple of tensors $(\mathbf{T}_1, \cdots, \mathbf{T}_S)$. In a multilayer perceptron with $L$ layers, the $S = 2L$ tensors $(\mathbf{T}_s)_{1\le s\le 2L}$ are the tensors of weights and the vectors of biases of each layer. In that case, $S \ll P$. Within this framework, we propose a technique *summarizing* the order-$d$ derivative of the loss, which is a tensor belonging to $\mathbb{R}^{P^d}$, into a tensor belonging to $\mathbb{R}^{S^d}$, which is significantly smaller and easier to compute.

**Second contribution: a scalable second-order optimization method.** Then, we make use of the preceding technique at orders 2 and 3 to build a second-order optimization method. Formally, the method presented here and Newton's method look alike: in both cases, a linear system $\mathbf{H}_0\mathbf{x} = \mathbf{g}_0$ has to be solved (according to $\mathbf{x}$), where $\mathbf{g}_0$ and $\mathbf{H}_0$ contain respectively first-order and second-order information about $\mathcal{L}$. Despite this formal resemblance, the difference is enormous: with Newton's method, $\mathbf{H}_0$ is equal to the Hessian $\mathbf{H}$ of $\mathcal{L}$ of size $P \times P$, while with ours, $\mathbf{H}_0$ is equal to a matrix $\bar{\mathbf{H}}$ of size $S \times S$. Thus, $\bar{\mathbf{H}}$ is undoubtedly smaller and easier to compute than $\mathbf{H}$ when $S \ll P$. Nevertheless, since $\bar{\mathbf{H}}$ is a dense matrix, it still contains information about the interactions between the tensors $\mathbf{T}_s$ when they are used in $\mathcal{L}$. This point is crucial, since most second-order optimization methods applied in neural networks use a simplified version of the Hessian (or

its inverse), usually a diagonal or block-diagonal ap-
proximation, which ignores the interactions between
layers. Finally, we propose an anisotropic version of
Nesterov's cubic regularization (Nesterov and Polyak,
2006), which uses order-3 information to regularize $\bar{\mathbf{H}}$
and avoid instabilities when computing $\bar{\mathbf{H}}^{-1}\bar{\mathbf{g}}$.

As a preview of the experimental results, we re-
port in Figure 1 the *summary of the Hessian* $\bar{\mathbf{H}}$ and
its inverse $\bar{\mathbf{H}}^{-1}$ at two stages of training of LeNet-5:
interactions between different layers exist and are not
negligible. And, to achieve the proof-of-concept, we
also provide in Section 5 a short series of experiments
showing that this method is able to achieve training
losses comparable to those of Adam and K-FAC.

<figure><img src="image_1.png" /><figcaption>Figure 1: LeNet-5 trained by SGD on CIFAR-10. Submatrices of $\bar{\mathbf{H}}$ (first row) and $\bar{\mathbf{H}}^{-1}$ (second row), where focus is on interactions: weight-weight, weight-bias, bias-bias of the different layers. We show it at initialization and before the 73rd epoch (best valid. loss). Legend: white: close-to-zero value; dark red: large positive value; dark blue: negative value, large in norm.</figcaption></figure>

**Structure of the paper.** First, we show the
context and motivation of our work in Section 2.
Then, we provide in Sections 3 and 4 standalone
presentations of the contributions, respectively the
higher-order information extraction technique and
the second-order optimization method. In Section
5, we present experimental results showing that the
entiation and partition of the set of parameters into
tensors stored on a GPU. Second, we aim to build
a second-order optimization method with properties
similar to the properties of Newton's method (but
inevitably weaker). Notably, we want a method in-
variant by layer-wise affine reparameterizations of the
model, providing a Hessian-inspired matrix showing
the interactions between each pair of tensors, along
with layer-wise step sizes. Naturally, the computa-
tional cost should remain reasonable, and the prob-
lem of close-to-zero eigenvalues of Hessian has to be
dealt with.

**First contribution: extracting higher-order information.** Formally, we study a loss $\mathcal{L}$ to minimize according to a vector of parameters $\theta \in \mathbb{R}^p$, which can be represented as a tuple of tensors $(\mathbf{T}_1, \cdots, \mathbf{T}_S)$. In a multilayer perceptron with $L$ layers, the $S = 2L$ tensors $(\mathbf{T}_s)_{1\le s\le 2L}$ are the tensors of weights and the vectors of biases of each layer. In that case, $S \ll P$. Within this framework, we propose a technique *summarizing* the order-$d$ derivative of the loss, which is a tensor belonging to $\mathbb{R}^{P^d}$, into a tensor belonging to $\mathbb{R}^{S^d}$, which is significantly smaller and easier to compute.

**Second contribution: a scalable second-order optimization method.** Then, we make use of the preceding technique at orders 2 and 3 to build a second-order optimization method. Formally, the method presented here and Newton's method look alike: in both cases, a linear system $\mathbf{H}_0\mathbf{x} = \mathbf{g}_0$ has to be solved (according to $\mathbf{x}$), where $\mathbf{g}_0$ and $\mathbf{H}_0$ contain respectively first-order and second-order information about $\mathcal{L}$. Despite this formal resemblance, the difference is enormous: with Newton's method, $\mathbf{H}_0$ is equal to the Hessian $\mathbf{H}$ of $\mathcal{L}$ of size $P \times P$, while with ours, $\mathbf{H}_0$ is equal to a matrix $\bar{\mathbf{H}}$ of size $S \times S$. Thus, $\bar{\mathbf{H}}$ is undoubtedly smaller and easier to compute than $\mathbf{H}$ when $S \ll P$. Nevertheless, since $\bar{\mathbf{H}}$ is a dense matrix, it still contains information about the interactions between the tensors $\mathbf{T}_s$ when they are used in $\mathcal{L}$. This point is crucial, since most second-order optimization methods applied in neural networks use a simplified version of the Hessian (or

its inverse), usually a diagonal or block-diagonal ap-
proximation, which ignores the interactions between
layers. Finally, we propose an anisotropic version of
Nesterov's cubic regularization (Nesterov and Polyak,
2006), which uses order-3 information to regularize $\bar{\mathbf{H}}$
and avoid instabilities when computing $\bar{\mathbf{H}}^{-1}\bar{\mathbf{g}}$.

As a preview of the experimental results, we re-
port in Figure 1 the *summary of the Hessian* $\bar{\mathbf{H}}$ and
its inverse $\bar{\mathbf{H}}^{-1}$ at two stages of training of LeNet-5:
interactions between different layers exist and are not
negligible. And, to achieve the proof-of-concept, we
also provide in Section 5 a short series of experiments
showing that this method is able to achieve training
losses comparable to those of Adam and K-FAC.

<figure><img src="image_3.png" /><figcaption>Figure 1: LeNet-5 trained by SGD on CIFAR-10.
Submatrices of $\bar{\mathbf{H}}$ (first row) and $\bar{\mathbf{H}}^{-1}$ (second row), where focus is on interactions: weight-weight, weight-bias, bias-bias of the different layers. We show it at initialization and before the 73rd epoch (best valid. loss). Legend: white: close-to-zero value; dark red: large positive value; dark blue: negative value, large in norm.</figcaption></figure>

**Structure of the paper.** First, we show the
context and motivation of our work in Section 2.
Then, we provide in Sections 3 and 4 standalone
presentations of the contributions, respectively the
higher-order information extraction technique and
the second-order optimization method. In Section
5, we present experimental results showing that the
developed methods are usable in practice. Finally, we discuss the results in Section 6.

# 2 Context and motivation

## 2.1 Higher-order information

Extracting the maximum higher-order information about a loss at a minimal computational cost and using it to improve optimization is not a novel idea. This is typically what is done by Dangel (2023), despite it does not go beyond the second-order derivative. In that research direction, the Hessian-vector product technique (Pearlmutter, 1994) is decisive, since it allows to compute at a very small cost the projection of higher-order derivatives in given directions (see Appendix A). Regarding the derivatives of order 3 and beyond, Nesterov's cubic regularization of Newton's method (Nesterov and Polyak, 2006) uses some information of order 3 to avoid too large training steps. Incidentally, we develop an anisotropic variation of it in Section 4.2.

Overall, the current need for information contained in higher-order derivatives is still limited. Anyway, we hope that making it available in large models would lead to the emergence of new uses.

## 2.2 Second-order methods

The Hessian $\mathbf{H}$ of the loss $\mathcal{L}$ according to the vector of parameters $\theta$ is known to contain useful information about $\mathcal{L}$. Above all, the Hessian is used to develop second-order optimization algorithms. Let us denote by $\theta_t$ the value of $\theta$ at time step $t$, $\mathbf{g}_t \in \mathbb{R}^P$ the gradient of $\mathcal{L}$ at step $t$ and $\mathbf{H}_t$ its Hessian at step $t$. One of the most widely known second-order optimization method is Newton's method, whose step is (Nocedal and Wright, 1999, Chap. 3.3):

$$ \theta_{t+1} := \theta_t - \mathbf{H}_t^{-1} \mathbf{g}_t. $$

Under some conditions, including strong convexity of $\mathcal{L}$, Newton's method is known to have a quadratic convergence rate (Nocedal and Wright, 1999, Th. 3.7), which makes this method very appealing. Besides, there exist other methods making use of

second-order information, without requiring the full computation of the Hessian. For instance, Cauchy's steepest descent (Cauchy, 1847) is a variation of the usual gradient descent, where the step size is tuned by extracting very little information from the Hessian:

$$ \theta_{t+1} := \theta_t - \eta_t^* \mathbf{g}_t, \quad \text{where} \quad \eta_t^* := \frac{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t}{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t}, $$

where the value of $\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t$ can be obtained with little computational cost (see Appendix A). However, when optimizing a quadratic function $f$ with Cauchy's steepest descent, $f(\theta_t)$ is known to decrease at a rate $(\frac{\lambda_{max}-\lambda_{min}}{\lambda_{max}+\lambda_{min}})^2$, where $\lambda_{max}$ and $\lambda_{min}$ are respectively the largest and the smallest eigenvalues of the Hessian of $f$ (Luenberger and Ye, 2008, Chap. 8.2, Th. 2). If the Hessian of $f$ is strongly anisotropic, then this rate is close to one and optimization is slow. For a comparison of both methods, one may refer to (Gill et al., 1981; Luenberger and Ye, 2008; Nocedal and Wright, 1999).

Finally, one can guess that there is some space between Newton's method, which requires to know the full Hessian, and Cauchy's steepest descent, which requires only minimal and computationally cheap information about the Hessian. The optimization method presented in Section 4 explores precisely this in-between space.

**Quasi-Newton methods.** In many high-dimensional situations, the computation of the Hessian $\mathbf{H}_t$, as well as the inversion of the linear system $\mathbf{g}_t = \mathbf{H}_t \mathbf{x}_t$, is computationally intensive. Quasi-Newton methods are designed to avoid any direct computation of the Hessian, and make an extensive use of the gradients and of finite difference methods to approximate the direction of $\mathbf{H}_t^{-1} \mathbf{g}_t$. A list of common quasi-Newton methods may be found in (Nocedal and Wright, 1999, Chap. 8).

As explained in Nocedal and Wright (1999), since it is easy to compute the Hessian by using Automatic Differentiation (AutoDiff), quasi-Newton methods tend to lose their interest. Nevertheless, they should remain useful in situations where such computation is too difficult. Besides, the availability of a tool as
developed methods are usable in practice. Finally, we discuss the results in Section 6.

# 2 Context and motivation

## 2.1 Higher-order information

Extracting the maximum higher-order information about a loss at a minimal computational cost and using it to improve optimization is not a novel idea. This is typically what is done by Dangel (2023), despite it does not go beyond the second-order derivative. In that research direction, the Hessian-vector product technique (Pearlmutter, 1994) is decisive, since it allows to compute at a very small cost the projection of higher-order derivatives in given directions (see Appendix A). Regarding the derivatives of order 3 and beyond, Nesterov's cubic regularization of Newton's method (Nesterov and Polyak, 2006) uses some information of order 3 to avoid too large training steps. Incidentally, we develop an anisotropic variation of it in Section 4.2.

Overall, the current need for information contained in higher-order derivatives is still limited. Anyway, we hope that making it available in large models would lead to the emergence of new uses.

## 2.2 Second-order methods

The Hessian $\mathbf{H}$ of the loss $\mathcal{L}$ according to the vector of parameters $\theta$ is known to contain useful information about $\mathcal{L}$. Above all, the Hessian is used to develop second-order optimization algorithms. Let us denote by $\theta_t$ the value of $\theta$ at time step $t$, $\mathbf{g}_t \in \mathbb{R}^P$ the gradient of $\mathcal{L}$ at step $t$ and $\mathbf{H}_t$ its Hessian at step $t$. One of the most widely known second-order optimization method is Newton's method, whose step is (Nocedal and Wright, 1999, Chap. 3.3):

$$ \theta_{t+1} := \theta_t - \mathbf{H}_t^{-1} \mathbf{g}_t. $$

Under some conditions, including strong convexity of $\mathcal{L}$, Newton's method is known to have a quadratic convergence rate (Nocedal and Wright, 1999, Th. 3.7), which makes this method very appealing. Besides, there exist other methods making use of

second-order information, without requiring the full computation of the Hessian. For instance, Cauchy's steepest descent (Cauchy, 1847) is a variation of the usual gradient descent, where the step size is tuned by extracting very little information from the Hessian:

$$ \theta_{t+1} := \theta_t - \eta_t^* \mathbf{g}_t, \quad \text{where} \quad \eta_t^* := \frac{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t}{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t}, $$

where the value of $\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t$ can be obtained with little computational cost (see Appendix A). However, when optimizing a quadratic function $f$ with Cauchy's steepest descent, $f(\theta_t)$ is known to decrease at a rate $(\frac{\lambda_{\max}-\lambda_{\min}}{\lambda_{\max}+\lambda_{\min}})^2$, where $\lambda_{\max}$ and $\lambda_{\min}$ are respectively the largest and the smallest eigenvalues of the Hessian of $f$ (Luenberger and Ye, 2008, Chap. 8.2, Th. 2). If the Hessian of $f$ is strongly anisotropic, then this rate is close to one and optimization is slow. For a comparison of both methods, one may refer to (Gill et al., 1981; Luenberger and Ye, 2008; Nocedal and Wright, 1999).

Finally, one can guess that there is some space between Newton's method, which requires to know the full Hessian, and Cauchy's steepest descent, which requires only minimal and computationally cheap information about the Hessian. The optimization method presented in Section 4 explores precisely this in-between space.

**Quasi-Newton methods.** In many high-dimensional situations, the computation of the Hessian $\mathbf{H}_t$, as well as the inversion of the linear system $\mathbf{g}_t = \mathbf{H}_t \mathbf{x}$, is computationally intensive. Quasi-Newton methods are designed to avoid any direct computation of the Hessian, and make an extensive use of the gradients and of finite difference methods to approximate the direction of $\mathbf{H}_t^{-1} \mathbf{g}_t$. A list of common quasi-Newton methods may be found in (Nocedal and Wright, 1999, Chap. 8).

As explained in Nocedal and Wright (1999), since it is easy to compute the Hessian by using Automatic Differentiation (AutoDiff), quasi-Newton methods tend to lose their interest. Nevertheless, they should remain useful in situations where such computation is too difficult. Besides, the availability of a tool as
powerful as AutoDiff should push us to invent methods avoiding the computation of the full Hessian.

**Block-diagonal approximations of the Hessian.**
In the development of neural network training algorithms, several attempts have been made to use the Hessian to build second-order optimization methods with good properties. A first rough attempt has been made by Wang and Lin (1998). In this work, the Hessian matrix is divided into blocks, following the division of the network into layers, and its off-diagonal blocks are removed. From another perspective, Olivier (2015) has kept this block-diagonal structure, but performed an additional approximation on the remaining blocks. In both cases, the interactions between the different layers are neglected.

**Kronecker-Factored Approximate Curvature (K-FAC).**
This technique of approximation of the Hessian has been proposed in (Martens and Grosse, 2015) in the context of neural network training. K-FAC leverages the specific architecture of neural networks to output an approximation of the true Hessian, which is usually too costly to compute. However, the resulting Hessian suffers from several issues. First, the main approximation is quite rough, since "[it assumes] statistical independence between products [...] of unit activities and products [...] of unit input derivatives" (Martens and Grosse, 2015, Sec. 3.1). Second, even with an approximation of the full Hessian, one still has to invert it, which is computationally intensive even for small networks. To overcome this difficulty, Martens and Grosse (2015) make a block-diagonal or block-tridiagonal approximation of the inverse of the Hessian, which kills many of the interactions between the parameters or the layers. However, this approximation is less rough than a direct block-diagonal approximation of the Hessian.

**Summarizing the Hessian.**
One should also mention Lu et al. (2018), which proposes to approximate the Hessian with a matrix composed of blocks in which all the coefficient are identical. Thus, the Hessian can be compressed into a smaller matrix, which looks like the summary of Hessian matrix $\mathbf{H}$ used

in Section 4. In a completely different setup, Yuan et al. (2022) proposes a “Sketched Newton-Raphson”, which is driven by the same spirit as the method presented in Section 4: instead of dealing with a complicated large matrix, one should “project” it on spaces of lower dimension.

**Optimization methods invariant by affine reparameterization.**
Several optimization methods based on the Hessian, such as Newton's method, have an optimization step invariant by affine reparameterization of $\theta$ (Amari, 1998) (Nesterov, 2003, Chap. 4.1.2). Specifically, when using Newton's method, it is equivalent to optimize $\mathcal{L}$ according to $\theta$ and according to $\tilde{\theta} = A\theta + B$ ($A \in \mathbb{R}^{P \times P}$ invertible, $B \in \mathbb{R}^P$). This affine-invariance property holds also when the function $\mathcal{L}$ to minimize is a negative log-likelihood, and one chooses to minimize $\theta$ with the natural gradient method (Amari, 1998). This method also requires, at some point, the computation of the Hessian of $\mathcal{L}$.

## 2.3 Motivation

**What are we really looking for?**
The methods presented in Section 2.2, which all attempt to construct the Hessian matrix $\mathbf{H}$ or its inverse $\mathbf{H}^{-1}$, do the same implicit assumption: we want to use a Newton-like optimization method, so we must have access to $\mathbf{H}$ or $\mathbf{H}^{-1}$. This assumption is certainly correct when the loss to optimize is strongly convex. But, when the loss is non convex and very complicated, for instance when training a neural network, this assumption lacks justification. Worse, it has been shown empirically that, at the end of the training of a neural network, the eigenvalues of the Hessian are concentrated around zero (Sagun et al., 2018), with only a few large positive eigenvalues.

Therefore, Newton's method itself does not seem to be recommended in neural network training, so we may not need to compute the full Hessian at all, which would relieve us from a tedious, if not impossible, task.<fnref n="1" />

<footnote n="1">¹"No one is bound to do the impossible."</footnote>
powerful as AutoDiff should push us to invent methods avoiding the computation of the full Hessian.

**Block-diagonal approximations of the Hessian.**
In the development of neural network training algorithms, several attempts have been made to use the Hessian to build second-order optimization methods with good properties. A first rough attempt has been made by Wang and Lin (1998). In this work, the Hessian matrix is divided into blocks, following the division of the network into layers, and its off-diagonal blocks are removed. From another perspective, Olivier (2015) has kept this block-diagonal structure, but performed an additional approximation on the remaining blocks. In both cases, the interactions between the different layers are neglected.

**Kronecker-Factored Approximate Curvature (K-FAC).**
This technique of approximation of the Hessian has been proposed in (Martens and Grosse, 2015) in the context of neural network training. K-FAC leverages the specific architecture of neural networks to output an approximation of the true Hessian, which is usually too costly to compute. However, the resulting Hessian suffers from several issues. First, the main approximation is quite rough, since "[it assumes] statistical independence between products [...] of unit activities and products [...] of unit input derivatives" (Martens and Grosse, 2015, Sec. 3.1). Second, even with an approximation of the full Hessian, one still has to invert it, which is computationally intensive even for small networks. To overcome this difficulty, Martens and Grosse (2015) make a block-diagonal or block-tridiagonal approximation of the inverse of the Hessian, which kills many of the interactions between the parameters or the layers. However, this approximation is less rough than a direct block-diagonal approximation of the Hessian.

**Summarizing the Hessian.**
One should also mention Lu et al. (2018), which proposes to approximate the Hessian with a matrix composed of blocks in which all the coefficient are identical. Thus, the Hessian can be compressed into a smaller matrix, which looks like the summary of Hessian matrix **H** used

in Section 4. In a completely different setup, Yuan et al. (2022) proposes a “Sketched Newton-Raphson”, which is driven by the same spirit as the method presented in Section 4: instead of dealing with a complicated large matrix, one should “project” it on spaces of lower dimension.

**Optimization methods invariant by affine reparameterization.**
Several optimization methods based on the Hessian, such as Newton's method, have an optimization step invariant by affine reparameterization of $\theta$ (Amari, 1998) (Nesterov, 2003, Chap. 4.1.2). Specifically, when using Newton's method, it is equivalent to optimize $\mathcal{L}$ according to $\theta$ and according to $\tilde{\theta} = A\theta + B$ ($A \in \mathbb{R}^{P \times P}$ invertible, $B \in \mathbb{R}^P$). This affine-invariance property holds also when the function $\mathcal{L}$ to minimize is a negative log-likelihood, and one chooses to minimize $\theta$ with the natural gradient method (Amari, 1998). This method also requires, at some point, the computation of the Hessian of $\mathcal{L}$.

## 2.3 Motivation

**What are we really looking for?**
The methods presented in Section 2.2, which all attempt to construct the Hessian matrix **H** or its inverse **H**⁻<fnref n="1" />, do the same implicit assumption: we want to use a Newton-like optimization method, so we must have access to **H** or **H**⁻<fnref n="1" />. This assumption is certainly correct when the loss to optimize is strongly convex. But, when the loss is non convex and very complicated, for instance when training a neural network, this assumption lacks justification. Worse, it has been shown empirically that, at the end of the training of a neural network, the eigenvalues of the Hessian are concentrated around zero (Sagun et al., 2018), with only a few large positive eigenvalues.

Therefore, Newton's method itself does not seem to be recommended in neural network training, so we may not need to compute the full Hessian at all, which would relieve us from a tedious, if not impossible, task.<fnref n="1" />

<footnote n="1">¹"No one is bound to do the impossible."</footnote>
**Importance of the interactions between layers.**

Also, we would like to show the importance of keeping the off-diagonal blocks in the Hessian or its inverse. Several empirical works have shown that the role and the behavior of each layer must be considered along its interactions between the other layers. We give two examples. First, Zhang et al. (2022) has shown that, at the end of their training, many networks exhibit a strange feature: some (but not all) layers can be reinitialized to their initial value with little damage to the performance. Second, Kornblith et al. (2019) has compared the similarity between the representations of the data after each layer. It has been remarked that the similarity matrix of the layers may change qualitatively when changing the number of layers (Kornblith et al., 2019, Fig. 3). Among all, these empirical results motivate our search for mathematical objects showing how layers interact.

**Per-layer scaling of the learning rates.**

A whole line of research is concerned with building a well-grounded method to find a good scaling of the distribution of initialization of the parameters, along with a good scaling for the learning rates, which may be chosen layer-wise. One may refer to the paper introducing the Neural Tangent Kernels (Jacot et al., 2018), in which a layer-wise scaling for the weights has been proposed and theoretically grounded. Also, in the “feature learning” line of works, (Yang and Hu, 2021) proposes a general relation between several scalings related to weight initialization and training. Therefore, there is an interest in finding a scalable and theoretically grounded method to build per-layer learning rates.

**Unleashing the power of AutoDiff.** Nowadays, several libraries provide easy-to-use automatic differentiation packages, allowing the user to compute numerically the gradient of a function, and even higher-order derivatives.<fnref n="2" /> Disregarding the computational cost, the full Hessian could theoretically be computed numerically without approximation. In order to make this computation feasible in practice, one should aim

for an easier goal: instead of computing the full Hessian, one may consider computing a smaller matrix, consisting of projections of the Hessian.

Moreover, we may hope that such projections would “squeeze” the close-to-zero eigenvalues of the Hessian, in such a way that the eigenvalues of the projected matrix would be mostly not too close to zero.

# 3 Summarizing higher-order information

Let us consider the minimization of a loss function $\mathcal{L}: \mathbb{R}^P \to \mathbb{R}$ according to a variable $\theta \in \mathbb{R}^P$.

**Full computation of the derivatives.** The object containing all the information of order $d$ about $\mathcal{L}$ is a $d$-linear form, which can be represented as a tensor of order $d$ as follows:

$$ \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta) : \mathbb{R}^P \times \dots \times \mathbb{R}^P \to \mathbb{R} \\ (\mathbf{u}^1, \dots, \mathbf{u}^d) \mapsto \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}^1, \dots, \mathbf{u}^d]. $$

This tensor belongs to $\mathbb{R}^{P^d}$ and contains $P^d$ scalars. Even when considering the symmetries, it is computationally too demanding to compute it exactly for $d \ge 2$ in most cases. For instance, it is not even possible to compute numerically the full Hessian of $\mathcal{L}$ according to the parameters of a small neural network, i.e., with $P = 10^5$ and $d = 2$, the Hessian contains $P^d = 10^{10}$ scalars.

**Terms of the Taylor expansion.** At the opposite, one can obtain cheap higher-order information about $\mathcal{L}$ at $\theta$ by considering a specific direction $\mathbf{u} \in \mathbb{R}^P$. The Taylor expansion of $\mathcal{L}(\theta + \mathbf{u})$ gives:

$$ \mathcal{L}(\theta + \mathbf{u}) = \mathcal{L}(\theta) + \sum_{d=1}^{D} \frac{1}{d!} \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}, \dots, \mathbf{u}] + o(\|\mathbf{u}\|^D). $$

<footnote n="2">²With PyTorch: `torch.autograd.grad`.</footnote>
**Importance of the interactions between layers.**
Also, we would like to show the importance of keeping the off-diagonal blocks in the Hessian or its inverse. Several empirical works have shown that the role and the behavior of each layer must be considered along its interactions between the other layers. We give two examples. First, Zhang et al. (2022) has shown that, at the end of their training, many networks exhibit a strange feature: some (but not all) layers can be reinitialized to their initial value with little damage to the performance. Second, Kornblith et al. (2019) has compared the similarity between the representations of the data after each layer. It has been remarked that the similarity matrix of the layers may change qualitatively when changing the number of layers (Kornblith et al., 2019, Fig. 3). Among all, these empirical results motivate our search for mathematical objects showing how layers interact.

**Per-layer scaling of the learning rates.** A whole line of research is concerned with building a well-grounded method to find a good scaling of the distribution of initialization of the parameters, along with a good scaling for the learning rates, which may be chosen layer-wise. One may refer to the paper introducing the Neural Tangent Kernels (Jacot et al., 2018), in which a layer-wise scaling for the weights has been proposed and theoretically grounded. Also, in the “feature learning” line of works, (Yang and Hu, 2021) proposes a general relation between several scalings related to weight initialization and training. Therefore, there is an interest in finding a scalable and theoretically grounded method to build per-layer learning rates.

**Unleashing the power of AutoDiff.** Nowadays, several libraries provide easy-to-use automatic differentiation packages, allowing the user to compute numerically the gradient of a function, and even higher-order derivatives.<fnref n="2" /> Disregarding the computational cost, the full Hessian could theoretically be computed numerically without approximation. In order to make this computation feasible in practice, one should aim

for an easier goal: instead of computing the full Hessian, one may consider computing a smaller matrix, consisting of projections of the Hessian.

Moreover, we may hope that such projections would “squeeze” the close-to-zero eigenvalues of the Hessian, in such a way that the eigenvalues of the projected matrix would be mostly not too close to zero.

# 3 Summarizing higher-order information

Let us consider the minimization of a loss function $\mathcal{L}: \mathbb{R}^P \to \mathbb{R}$ according to a variable $\theta \in \mathbb{R}^P$.

**Full computation of the derivatives.** The object containing all the information of order $d$ about $\mathcal{L}$ is a $d$-linear form, which can be represented as a tensor of order $d$ as follows:

$$ \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta) : \mathbb{R}^P \times \dots \times \mathbb{R}^P \to \mathbb{R} \\ (\mathbf{u}^1, \dots, \mathbf{u}^d) \mapsto \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}^1, \dots, \mathbf{u}^d]. $$

This tensor belongs to $\mathbb{R}^{P^d}$ and contains $P^d$ scalars. Even when considering the symmetries, it is computationally too demanding to compute it exactly for $d \ge 2$ in most cases. For instance, it is not even possible to compute numerically the full Hessian of $\mathcal{L}$ according to the parameters of a small neural network, i.e., with $P = 10^5$ and $d = 2$, the Hessian contains $P^d = 10^{10}$ scalars.

**Terms of the Taylor expansion.** At the opposite, one can obtain cheap higher-order information about $\mathcal{L}$ at $\theta$ by considering a specific direction $\mathbf{u} \in \mathbb{R}^P$. The Taylor expansion of $\mathcal{L}(\theta + \mathbf{u})$ gives:

$$ \mathcal{L}(\theta + \mathbf{u}) = \mathcal{L}(\theta) + \sum_{d=1}^{D} \frac{1}{d!} \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}, \dots, \mathbf{u}] + o(\|\mathbf{u}\|^D). $$

<footnote n="2">²With PyTorch: `torch.autograd.grad`.</footnote>
The terms of the Taylor expansion contain higher-order information about $\mathcal{L}$ in the direction $\mathbf{u}$. Notably, they can be used to predict how $\mathcal{L}(\theta)$ would change when translating $\theta$ in the direction of $\mathbf{u}$. Additionally, computing the first $D$ terms has a complexity of order $D \times P$, which is manageable even for large models. The trick allowing for such a low complexity, the *Hessian-vector product*, has been proposed by *Pearlmutter* (1994) and is recalled in Appendix A.

**An intermediate solution.** Now, let us assume that, in the practical implementation of a gradient-based method of optimization of $\mathcal{L}(\theta)$, $\theta$ is represented by a tuple of tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$. So, each Taylor term can be expressed as:

$$ 
\begin{aligned}
& \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}, \dots, \mathbf{u}] \\
&= \sum_{s_1=1}^{S} \dots \sum_{s_d=1}^{S} \frac{\partial^d \mathcal{L}}{\partial \mathbf{T}_{s_1} \dots \partial \mathbf{T}_{s_d}}(\theta)[\mathbf{U}^{s_1}, \dots, \mathbf{U}^{s_d}] \\
&= \mathbf{D}_{\theta}^{(d)}(\mathbf{u})[\mathbf{1}_S, \dots, \mathbf{1}_S],
\end{aligned}
$$

where $\mathbf{1}_S \in \mathbb{R}^S$ is a vector full of ones, the tuple of tensors $(\mathbf{U}^1, \dots, \mathbf{U}^S)$ represents $\mathbf{u}$,<fnref n="3" /> and $\mathbf{D}_{\theta}^{(d)}(\mathbf{u}) \in \mathbb{R}^{S^d}$ is a tensor of order $d$ with size $S$ in every dimension with values:

$$ (\mathbf{D}_{\theta}^{(d)}(\mathbf{u}))_{s_1, \dots, s_d} = \frac{\partial^d \mathcal{L}}{\partial \mathbf{T}_{s_1} \dots \partial \mathbf{T}_{s_d}}(\theta)[\mathbf{U}^{s_1}, \dots, \mathbf{U}^{s_d}]. $$

Thus, $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ is equivalently a $d$-linear form on $\mathbb{R}^S$ and a tensor of order $d$ and size $S$ in every dimension. Moreover, the trick of *Pearlmutter* (1994) applies also to the computation of $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$, which is then much less costly to compute than the Hessian.

**Properties of $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$.** We show a comparison between the three techniques in Table 1. If $S$ is small enough, computing $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ becomes feasible for $d \ge 2$. For usual multilayer perceptrons with $L$ layers, there is one tensor of weights and one vector of biases per layer, so $S = 2L$. This allows to compute $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ in practice for $d = 2$ even when $L \approx 20$.

<table><caption>Table 1: Comparison between three techniques extracting higher-order information about $\mathcal{L}$: size of the result and complexity of the computation.</caption><thead><tr><th></th><th>size</th><th>complexity</th></tr></thead><tbody><tr><td>Full derivative <sup>∂<sup>d</sup>ζ</sup>/<sub>∂θ<sup>d</sup></sub>(θ)</td><td>P<sup>d</sup></td><td>P<sup>d</sup></td></tr><tr><td>Taylor term <b>D</b><sub>θ</sub><sup>(d)</sup>(<b>u</b>)[1<sub>S</sub>, ..., ]</td><td>1</td><td>d × P</td></tr><tr><td>Tensor <b>D</b><sub>θ</sub><sup>(d)</sup>(<b>u</b>)</td><td>S<sup>d</sup></td><td>S<sup>d-1</sup> × P</td></tr></tbody></table>

Therefore, the tensors $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ extract more information than the naive Taylor terms, while keeping a reasonable computational cost. Moreover, through their off-diagonal elements, they give access to information about the interactions between the tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$ when there are processed in the function $\mathcal{L}$.

# 4 A scalable second-order optimization method

## 4.1 Presentation of the method

The method presented here consists in partitioning the set of indices of parameters $\{1, \dots, P\}$ into $S$ subsets $(\mathcal{I}_s)_{1 \le s \le S}$, attribute for all $1 \le s \le S$ the same learning rate $\eta_s$ to the parameters $(\theta_p)_{p \in \mathcal{I}_s}$, and find the vector of learning rates $\boldsymbol{\eta} = (\eta_1, \dots, \eta_S)$ optimizing the decrease of the loss $\mathcal{L}$ for the current training step $t$, by using its order-2 Taylor approximation.<fnref n="4" /> Formally, given a direction $\mathbf{u}_t \in \mathbb{R}^P$ in the parameter space (typically, $\mathbf{u}_t = \mathbf{g}_t$, the gradient) and $\mathbf{U}_t := \text{Diag}(\mathbf{u}_t) \in \mathbb{R}^{P \times P}$, we consider the training step:

$$ \boldsymbol{\theta}_{t+1} := \boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}_t, $$

which is a training step in a direction based on $\mathbf{u}_t$, distorted by a subset-wise step size $\boldsymbol{\eta}_t$. Then, we minimize the order-2 Taylor approximation $\boldsymbol{\Delta}_2$ of

<footnote n="3">³($\mathbf{U}^1, \dots, \mathbf{U}^S$) is to $\mathbf{u}$ as $(\mathbf{T}_1, \dots, \mathbf{T}_S)$ is to $\theta$.</footnote>

<footnote n="4">⁴With the notation of Section 3, $\mathcal{I}_s$ is the set of indices $p$ of the parameters $\theta_p$ belonging to the tensor $\mathbf{T}_s$, so the scalars $(\theta_p)_{p \in \mathcal{I}_s}$ correspond exactly to the scalars belonging to $\mathbf{T}_s$. Consequently, everything is as if a specific learning rate $\eta_s$ is assigned to each tensor $\mathbf{T}_s$.</footnote>
The terms of the Taylor expansion contain higher-order information about $\mathcal{L}$ in the direction $\mathbf{u}$. Notably, they can be used to predict how $\mathcal{L}(\theta)$ would change when translating $\theta$ in the direction of $\mathbf{u}$. Additionally, computing the first $D$ terms has a complexity of order $D \times P$, which is manageable even for large models. The trick allowing for such a low complexity, the *Hessian-vector product*, has been proposed by *Pearlmutter* (1994) and is recalled in Appendix A.

**An intermediate solution.** Now, let us assume that, in the practical implementation of a gradient-based method of optimization of $\mathcal{L}(\theta)$, $\theta$ is represented by a tuple of tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$. So, each Taylor term can be expressed as:

$$ 
\begin{aligned}
& \frac{\partial^d \mathcal{L}}{\partial \theta^d}(\theta)[\mathbf{u}, \dots, \mathbf{u}] \\
&= \sum_{s_1=1}^{S} \dots \sum_{s_d=1}^{S} \frac{\partial^d \mathcal{L}}{\partial \mathbf{T}_{s_1} \dots \partial \mathbf{T}_{s_d}}(\theta)[\mathbf{U}^{s_1}, \dots, \mathbf{U}^{s_d}] \\
&= \mathbf{D}_{\theta}^{(d)}(\mathbf{u})[\mathbf{1}_S, \dots, \mathbf{1}_S],
\end{aligned}
$$

where $\mathbf{1}_S \in \mathbb{R}^S$ is a vector full of ones, the tuple of tensors $(\mathbf{U}^1, \dots, \mathbf{U}^S)$ represents $\mathbf{u}$,<fnref n="3" /> and $\mathbf{D}_{\theta}^{(d)}(\mathbf{u}) \in \mathbb{R}^{S^d}$ is a tensor of order $d$ with size $S$ in every dimension with values:

$$ (\mathbf{D}_{\theta}^{(d)}(\mathbf{u}))_{s_1, \dots, s_d} = \frac{\partial^d \mathcal{L}}{\partial \mathbf{T}_{s_1} \dots \partial \mathbf{T}_{s_d}}(\theta)[\mathbf{U}^{s_1}, \dots, \mathbf{U}^{s_d}]. $$

Thus, $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ is equivalently a $d$-linear form on $\mathbb{R}^S$ and a tensor of order $d$ and size $S$ in every dimension. Moreover, the trick of *Pearlmutter* (1994) applies also to the computation of $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$, which is then much less costly to compute than the Hessian.

**Properties of $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$.** We show a comparison between the three techniques in Table 1. If $S$ is small enough, computing $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ becomes feasible for $d \ge 2$. For usual multilayer perceptrons with $L$ layers, there is one tensor of weights and one vector of biases per layer, so $S = 2L$. This allows to compute $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ in practice for $d = 2$ even when $L \approx 20$.

<table><caption>Table 1: Comparison between three techniques extracting higher-order information about $\mathcal{L}$: size of the result and complexity of the computation.</caption><thead><tr><th></th><th>size</th><th>complexity</th></tr></thead><tbody><tr><td>Full derivative <sup>∂<sup>d</sup>ζ</sup>/<sub>∂θ<sup>d</sup></sub>(θ)</td><td>P<sup>d</sup></td><td>P<sup>d</sup></td></tr><tr><td>Taylor term D<sub>θ</sub><sup>(d)</sup>(u)[1<sub>S</sub>, ..., ]</td><td>1</td><td>d × P</td></tr><tr><td>Tensor D<sub>θ</sub><sup>(d)</sup>(u)</td><td>S<sup>d</sup></td><td>S<sup>d-1</sup> × P</td></tr></tbody></table>

Therefore, the tensors $\mathbf{D}_{\theta}^{(d)}(\mathbf{u})$ extract more information than the naive Taylor terms, while keeping a reasonable computational cost. Moreover, through their off-diagonal elements, they give access to information about the interactions between the tensors $(\mathbf{T}_1, \dots, \mathbf{T}_S)$ when there are processed in the function $\mathcal{L}$.

# 4 A scalable second-order optimization method

## 4.1 Presentation of the method

The method presented here consists in partitioning the set of indices of parameters $\{1, \dots, P\}$ into $S$ subsets $(\mathcal{I}_s)_{1 \le s \le S}$, attribute for all $1 \le s \le S$ the same learning rate $\eta_s$ to the parameters $(\theta_p)_{p \in \mathcal{I}_s}$, and find the vector of learning rates $\boldsymbol{\eta} = (\eta_1, \dots, \eta_S)$ optimizing the decrease of the loss $\mathcal{L}$ for the current training step $t$, by using its order-2 Taylor approximation.<fnref n="4" /> Formally, given a direction $\mathbf{u}_t \in \mathbb{R}^P$ in the parameter space (typically, $\mathbf{u}_t = \mathbf{g}_t$, the gradient) and $\mathbf{U}_t := \text{Diag}(\mathbf{u}_t) \in \mathbb{R}^{P \times P}$, we consider the training step:

$$ \boldsymbol{\theta}_{t+1} := \boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}_t, $$

which is a training step in a direction based on $\mathbf{u}_t$, distorted by a subset-wise step size $\boldsymbol{\eta}_t$. Then, we minimize the order-2 Taylor approximation $\boldsymbol{\Delta}_2$ of

<footnote n="3">³($\mathbf{U}^1, \dots, \mathbf{U}^S$) is to $\mathbf{u}$ as $(\mathbf{T}_1, \dots, \mathbf{T}_S)$ is to $\theta$.</footnote>

<footnote n="4">⁴With the notation of Section 3, $\mathcal{I}_s$ is the set of indices $p$ of the parameters $\theta_p$ belonging to the tensor $\mathbf{T}_s$, so the scalars $(\theta_p)_{p \in \mathcal{I}_s}$ correspond exactly to the scalars belonging to $\mathbf{T}_s$. Consequently, everything is as if a specific learning rate $\eta_s$ is assigned to each tensor $\mathbf{T}_s$.</footnote>
$\mathcal{L}(\theta_{t+1}) - \mathcal{L}(\theta_t):$

$$ \Delta_2(\eta_t) := -\mathbf{g}_t^T \mathbf{U}_t \mathbf{I}_{P:S} \eta_t + \frac{1}{2} \eta_t^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \eta_t, $$

which gives:

$$ \theta_{t+1} = \theta_t - \mathbf{U}_t \mathbf{I}_{P:S} \eta_t^*, \quad (1) $$

$$ \eta_t^* := (\mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S})^{-1} \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t, $$

where $\mathbf{I}_{S:P} \in \mathbb{R}^{S \times P}$ is the partition matrix, verifying $(\mathbf{I}_{S:P})_{sp} = 1$ if $p \in \mathcal{I}_s$ and 0 otherwise, and $\mathbf{I}_{P:S} := \mathbf{I}_{S:P}^T$. Alternatively, $\eta_t^*$ can be written:

$$ \eta_t^* = \bar{\mathbf{H}}_t^{-1} \bar{\mathbf{g}}_t, \quad \bar{\mathbf{H}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \quad \mathbf{g}_t \in \mathbb{R}^{S \times S}, \\ \bar{\mathbf{g}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t \quad \mathbf{g}_t \in \mathbb{R}^S. $$

Details are provided in Appendix B.

With the notation of Section 3, $\bar{\mathbf{H}}_t = \mathbf{D}_{\theta_t}^{(2)}(\mathbf{u}_t)$ and $\bar{\mathbf{g}}_t = \mathbf{D}_{\theta_t}^{(1)}(\mathbf{u}_t)$. Incidentally, computing $\bar{\mathbf{H}}$ is of complexity $SP$, and solving the system $\bar{\mathbf{H}}\mathbf{x} = \bar{\mathbf{g}}$ is of complexity $S^2$.

## 4.2 Regularizing $\bar{\mathbf{H}}$ by using order-3 information

The method proposed in Section 4.1 requires to compute $\eta^* = \bar{\mathbf{H}}^{-1}\bar{\mathbf{g}}$. Usually, inverting such a linear system at every step is considered as hazardous and unstable. That is why, when using Newton's method, instead of computing a direction of descent $\mathbf{u} := \mathbf{H}^{-1}\mathbf{g}$, it is very common to add a regularization term: $\mathbf{u}_{\lambda} := (\mathbf{H} + \lambda\mathbf{I})^{-1}\mathbf{g}$ (Nocedal and Wright, 1999, Chap. 6.3).

However, this regularization technique suffers from several drawbacks. The most obvious is computational: finding the best hyperparameter $\lambda$ requires either many runs of the full optimization process, either one run with a technique adjusting $\lambda$ on the fly. Anyway, there is an inevitable extra computational cost.

Besides of that, the theoretical ground of such a regularization technique is not fully satisfactory. Basically, the main problem is not having a matrix $\bar{\mathbf{H}}$ with close-to-zero eigenvalues: after all, if the loss landscape is very flat in a specific direction, it is better to make a large training step. The problem lies

in the order-2 approximation of the loss made in the training step (1), as well as in Newton's method: instead of optimizing the true decrease of the loss, we optimize the decrease of its order-2 approximation. Thus, the practical question is: does this approximation model faithfully the loss at the current point $\theta_t$, in a region encompassing also the next point $\theta_{t+1}$?

To answer that question, one has to take into account order-3 information, and regularize $\bar{\mathbf{H}}$ in such a way that the resulting update remains in an area around $\theta_t$ where the cubic term of the Taylor approximation is negligible. In practice, we propose an anisotropic version of Nesterov's cubic regularization Nesterov and Polyak (2006).

**Anisotropic Nesterov cubic regularization.** By using the technique presented in Section 3, the diagonal coefficients $(D_1, \dots, D_S)$ of $\mathbf{D}_{\theta}^{(3)}(\mathbf{u}) \in \mathbb{R}^{S \times S \times S}$ are available with little computational cost. Let:

$$ \mathbf{D} := \mathrm{Diag}(|D_1|^{1/3}, \dots, |D_S|^{1/3}) \in \mathbb{R}^S. $$

We modify the method of Nesterov and Polyak (2006) by integrating an anisotropic factor $\mathbf{D}$ into the cubic term. Thus, our goal is to minimize according to $\eta$ the function $T$:

$$ T(\eta) := -\eta^T \bar{\mathbf{g}} + \frac{1}{2} \eta \bar{\mathbf{H}} \eta + \frac{\lambda_{\mathrm{int}}}{6} \|\mathbf{D}\eta\|^3, $$

where $\lambda_{\mathrm{int}}$ is the internal damping coefficient, which can be used to tune the strength of the cubic regularization. Under conditions detailed in Appendix D, this minimization problem is equivalent to finding a solution $\eta_*$ such that:

$$ \eta_* = \left( \bar{\mathbf{H}} + \frac{\lambda_{\mathrm{int}}}{2} \|\mathbf{D}\eta_*\| \mathbf{D}^2 \right)^{-1} \bar{\mathbf{g}}, \quad (2) $$

which is no more than a regularized version of (1). Finally, this multi-dimensional minimization problem boils down to a scalar root finding problem (see Appendix D).

## 4.3 Properties

The final method is a combination of the training step (1) with regularization (2):
$\mathcal{L}(\theta_{t+1}) - \mathcal{L}(\theta_t):$

$$ \Delta_2(\eta_t) := -\mathbf{g}_t^T \mathbf{U}_t \mathbf{I}_{P:S} \eta_t + \frac{1}{2} \eta_t^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \eta_t, $$

which gives:

$$ \theta_{t+1} = \theta_t - \mathbf{U}_t \mathbf{I}_{P:S} \eta_t^*, \quad (1) $$

$$ \eta_t^* := (\mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S})^{-1} \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t, $$

where $\mathbf{I}_{S:P} \in \mathbb{R}^{S \times P}$ is the partition matrix, verifying $(\mathbf{I}_{S:P})_{sp} = 1$ if $p \in \mathcal{I}_s$ and 0 otherwise, and $\mathbf{I}_{P:S} := \mathbf{I}_{S:P}^T$. Alternatively, $\eta_t^*$ can be written:

$$ \begin{aligned} \eta_t^* &= \bar{\mathbf{H}}_t^{-1} \bar{\mathbf{g}}_t, & \bar{\mathbf{H}}_t &:= \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} & \in \mathbb{R}^{S \times S}, \\ \bar{\mathbf{g}}_t &:= \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t & & \in \mathbb{R}^S. \end{aligned} $$

Details are provided in Appendix B.

With the notation of Section 3, $\bar{\mathbf{H}}_t = \mathbf{D}_{\theta_t}^{(2)}(\mathbf{u}_t)$ and $\bar{\mathbf{g}}_t = \mathbf{D}_{\theta_t}^{(1)}(\mathbf{u}_t)$. Incidentally, computing $\bar{\mathbf{H}}$ is of complexity $SP$, and solving the system $\bar{\mathbf{H}}\mathbf{x} = \bar{\mathbf{g}}$ is of complexity $S^2$.

## 4.2 Regularizing $\bar{\mathbf{H}}$ by using order-3 information

The method proposed in Section 4.1 requires to compute $\eta^* = \bar{\mathbf{H}}^{-1}\bar{\mathbf{g}}$. Usually, inverting such a linear system at every step is considered as hazardous and unstable. That is why, when using Newton's method, instead of computing a direction of descent $\mathbf{u} := \mathbf{H}^{-1}\mathbf{g}$, it is very common to add a regularization term: $\mathbf{u}_{\lambda} := (\mathbf{H} + \lambda\mathbf{I})^{-1}\mathbf{g}$ (Nocedal and Wright, 1999, Chap. 6.3).

However, this regularization technique suffers from several drawbacks. The most obvious is computational: finding the best hyperparameter $\lambda$ requires either many runs of the full optimization process, either one run with a technique adjusting $\lambda$ on the fly. Anyway, there is an inevitable extra computational cost.

Besides of that, the theoretical ground of such a regularization technique is not fully satisfactory. Basically, the main problem is not having a matrix $\bar{\mathbf{H}}$ with close-to-zero eigenvalues: after all, if the loss landscape is very flat in a specific direction, it is better to make a large training step. The problem lies

in the order-2 approximation of the loss made in the training step (1), as well as in Newton's method: instead of optimizing the true decrease of the loss, we optimize the decrease of its order-2 approximation. Thus, the practical question is: does this approximation model faithfully the loss at the current point $\theta_t$, in a region encompassing also the next point $\theta_{t+1}$?

To answer that question, one has to take into account order-3 information, and regularize $\bar{\mathbf{H}}$ in such a way that the resulting update remains in an area around $\theta_t$ where the cubic term of the Taylor approximation is negligible. In practice, we propose an anisotropic version of Nesterov's cubic regularization Nesterov and Polyak (2006).

**Anisotropic Nesterov cubic regularization.** By using the technique presented in Section 3, the diagonal coefficients $(D_1, \dots, D_S)$ of $\mathbf{D}_{\theta}^{(3)}(\mathbf{u}) \in \mathbb{R}^{S \times S \times S}$ are available with little computational cost. Let:

$$ \mathbf{D} := \mathrm{Diag}(|D_1|^{1/3}, \dots, |D_S|^{1/3}) \in \mathbb{R}^S. $$

We modify the method of Nesterov and Polyak (2006) by integrating an anisotropic factor $\mathbf{D}$ into the cubic term. Thus, our goal is to minimize according to $\eta$ the function $T$:

$$ T(\eta) := -\eta^T \bar{\mathbf{g}} + \frac{1}{2} \eta \bar{\mathbf{H}} \eta + \frac{\lambda_{\mathrm{int}}}{6} \|\mathbf{D}\eta\|^3, $$

where $\lambda_{\mathrm{int}}$ is the internal damping coefficient, which can be used to tune the strength of the cubic regularization. Under conditions detailed in Appendix D, this minimization problem is equivalent to finding a solution $\eta_*$ such that:

$$ \eta_* = \left( \bar{\mathbf{H}} + \frac{\lambda_{\mathrm{int}}}{2} \|\mathbf{D}\eta_*\| \mathbf{D}^2 \right)^{-1} \bar{\mathbf{g}}, \quad (2) $$

which is no more than a regularized version of (1). Finally, this multi-dimensional minimization problem boils down to a scalar root finding problem (see Appendix D).

## 4.3 Properties

The final method is a combination of the training step (1) with regularization (2):
**Method 4.1.** Training step $\theta_{t+1} = \theta_t - U_t I_{P:S} \eta_t^*$, where $\eta_t$ is the solution with the largest norm $\|D\eta\|$ of the equation:

$$ \eta = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|D\eta\| D^2 \right)^{-1} \bar{\mathbf{g}}. $$

Method 4.1 has several interesting properties.

**Encompassing Newton's method and Cauchy's steepest descent.** Without the cubic regularization ($\lambda_{\text{int}} = 0$), Newton's method is recovered when using the discrete partition, that is, $S = P$ with $\mathcal{I}_s = \{s\}$ for all $s$, and Cauchy's steepest descent is recovered when using the trivial partition, that is, $S = 1$ with $\mathcal{I}_1 = \{1, \dots, P\}$. See Appendix C for more details.

**No need to compute or approximate the full Hessian.** The full computation of the Hessian $\mathbf{H}_t \in \mathbb{R}^{P \times P}$ is not required. Instead, one only needs to compute the $S \times S$ matrix $\bar{\mathbf{H}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S}$, which can be done efficiently by computing $\mathbf{u}^T \mathbf{H}_t \mathbf{v}$ for a number $S \times S$ of pairs of well-chosen directions $(\mathbf{u}, \mathbf{v}) \in \mathbb{R}^P \times \mathbb{R}^P$. This property is useful especially when $S \ll P$. When optimizing a neural network with $L = 10$ layers and $P = 10^6$ parameters, one can naturally partition the set of parameters into $S = 2L$ subsets, each one containing either all the weights or all the biases of each of the $L$ layers. In this situation, one has to solve a linear system of size $2L = 20$ at each step, which is much more reasonable than solving a linear system with $P = 10^6$ equations. We call this natural partition of the parameters of a neural network the canonical partition.

**No need to solve a very large linear system.** Using Equations (1) or (2) only requires solving a linear system of $S$ equations, instead of $P$ in Newton's method. With the cubic regularization, only a constant term is added to the complexity, since it is a matter of scalar root finding.

**The interactions between different tensors are not neglected.** The matrix $\bar{\mathbf{H}}_t$, which simulates

the Hessian $\mathbf{H}_t$, is basically dense, which means that it does not exhibit a diagonal or block-diagonal structure. So, the interactions between subsets of parameters are taken into account when performing optimization steps. In the framework of neural networks with the canonical partition, it means that interactions between layers are taken into account during optimization, including even when the layers are far from each other. This is a great advantage compared to many existing approximations of the Hessian or its inverse, which are diagonal or block-diagonal.

**Invariance by subset-wise affine reparameterization.** As shown in Appendix E, under a condition on the directions $\mathbf{u}_t$, the trajectory of optimization of a model trained by Method 4.1 is invariant by affine reparameterization of the sub-vectors of parameters $\theta_{\mathcal{I}_s} := \text{vec}(\{\theta_p : p \in \mathcal{I}_s\})$.<fnref n="5" /> More precisely, let $(\alpha_s)_{1 \le s \le S}$ and $(\beta_s)_{1 \le s \le S}$ be respectively a sequence of nonzero scalings and a sequence of offsets, and $\tilde{\theta}$ such that, for all $1 \le s \le S$, $\tilde{\theta}_{\mathcal{I}_s} = \alpha_s \theta_{\mathcal{I}_s} + \beta_s$. Then, the training trajectory of the model is the same with both parameterizations $\theta$ and $\tilde{\theta}$. This property is interesting in the case of neural networks, where one can either use the usual parameterization or the NTK parameterization, which consists in a layer-wise scaling of the parameters.

Compared to the standard regularization $\bar{\mathbf{H}} + \lambda \mathbf{I}$ and Nesterov's cubic regularization, the anisotropic Nesterov regularization does not break the property of invariance by subset-wise scaling of the parameters of (1). Among all, this is due to our choice to keep only the diagonal coefficients of $\mathbf{D}_{\theta}^{(3)}(\mathbf{u})$ while discarding the others. Notably, the off-diagonal coefficients contain cross-derivatives, which would be difficult to include in an invariant training step.

# 5 Experiments

## 5.1 Empirical computation of $\bar{\mathbf{H}}$ and $\eta$

As recalled in Section 2, many works perform a diagonal, block-diagonal or block-tridiagonal (Martens

<footnote n="5">⁵This holds typically if $\mathbf{u}_t$ is the gradient or a moving average of the gradients (momentum).</footnote>
**Method 4.1.** Training step $\theta_{t+1} = \theta_t - U_t I_{P:S} \eta_t^*$, where $\eta_t$ is the solution with the largest norm $\|D\eta\|$ of the equation:

$$ \eta = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|D\eta\| D^2 \right)^{-1} \bar{\mathbf{g}}. $$

Method 4.1 has several interesting properties.

**Encompassing Newton's method and Cauchy's steepest descent.** Without the cubic regularization ($\lambda_{\text{int}} = 0$), Newton's method is recovered when using the discrete partition, that is, $S = P$ with $\mathcal{I}_s = \{s\}$ for all $s$, and Cauchy's steepest descent is recovered when using the trivial partition, that is, $S = 1$ with $\mathcal{I}_1 = \{1, \dots, P\}$. See Appendix C for more details.

**No need to compute or approximate the full Hessian.** The full computation of the Hessian $\mathbf{H}_t \in \mathbb{R}^{P \times P}$ is not required. Instead, one only needs to compute the $S \times S$ matrix $\bar{\mathbf{H}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S}$, which can be done efficiently by computing $\mathbf{u}^T \mathbf{H}_t \mathbf{v}$ for a number $S \times S$ of pairs of well-chosen directions $(\mathbf{u}, \mathbf{v}) \in \mathbb{R}^P \times \mathbb{R}^P$. This property is useful especially when $S \ll P$. When optimizing a neural network with $L = 10$ layers and $P = 10^6$ parameters, one can naturally partition the set of parameters into $S = 2L$ subsets, each one containing either all the weights or all the biases of each of the $L$ layers. In this situation, one has to solve a linear system of size $2L = 20$ at each step, which is much more reasonable than solving a linear system with $P = 10^6$ equations. We call this natural partition of the parameters of a neural network the canonical partition.

**No need to solve a very large linear system.** Using Equations (1) or (2) only requires solving a linear system of $S$ equations, instead of $P$ in Newton's method. With the cubic regularization, only a constant term is added to the complexity, since it is a matter of scalar root finding.

**The interactions between different tensors are not neglected.** The matrix $\bar{\mathbf{H}}_t$, which simulates

the Hessian $\mathbf{H}_t$, is basically dense, which means that it does not exhibit a diagonal or block-diagonal structure. So, the interactions between subsets of parameters are taken into account when performing optimization steps. In the framework of neural networks with the canonical partition, it means that interactions between layers are taken into account during optimization, including even when the layers are far from each other. This is a great advantage compared to many existing approximations of the Hessian or its inverse, which are diagonal or block-diagonal.

**Invariance by subset-wise affine reparameterization.** As shown in Appendix E, under a condition on the directions $\mathbf{u}_t$, the trajectory of optimization of a model trained by Method 4.1 is invariant by affine reparameterization of the sub-vectors of parameters $\theta_{\mathcal{I}_s} := \text{vec}(\{\theta_p : p \in \mathcal{I}_s\})$.<fnref n="5" /> More precisely, let $(\alpha_s)_{1 \le s \le S}$ and $(\beta_s)_{1 \le s \le S}$ be respectively a sequence of nonzero scalings and a sequence of offsets, and $\tilde{\theta}$ such that, for all $1 \le s \le S$, $\tilde{\theta}_{\mathcal{I}_s} = \alpha_s \theta_{\mathcal{I}_s} + \beta_s$. Then, the training trajectory of the model is the same with both parameterizations $\theta$ and $\tilde{\theta}$. This property is interesting in the case of neural networks, where one can either use the usual parameterization or the NTK parameterization, which consists in a layer-wise scaling of the parameters.

Compared to the standard regularization $\bar{\mathbf{H}} + \lambda \mathbf{I}$ and Nesterov's cubic regularization, the anisotropic Nesterov regularization does not break the property of invariance by subset-wise scaling of the parameters of (1). Among all, this is due to our choice to keep only the diagonal coefficients of $\mathbf{D}_{\theta}^{(3)}(\mathbf{u})$ while discarding the others. Notably, the off-diagonal coefficients contain cross-derivatives, which would be difficult to include in an invariant training step.

# 5 Experiments

## 5.1 Empirical computation of $\bar{\mathbf{H}}$ and $\eta$

As recalled in Section 2, many works perform a diagonal, block-diagonal or block-tridiagonal (Martens

<footnote n="5">⁵This holds typically if $\mathbf{u}_t$ is the gradient or a moving average of the gradients (momentum).</footnote>
and Grosse, 2015) approximation of the Hessian or its inverse.

Since a summary $\bar{\mathbf{H}}$ of the Hessian and its inverse $\bar{\mathbf{H}}^{-1}$ are available and all their off-diagonal coefficients have been computed and kept, one can to check if these coefficient are indeed negligible.

**Setup.** We have trained LeNet-5 and VGG-11$^{6}$ on CIFAR-10 by SGD with momentum. Before each epoch, we compute the full-batch gradient, denoted by $\mathbf{u}$, which we use as a direction to compute $\bar{\mathbf{H}}$, again in full-batch. We report submatrices of $\bar{\mathbf{H}}$ and $\bar{\mathbf{H}}^{-1}$ at initialization and at the epoch where the validation loss is the best in Figure 1 (LeNet) and Figure 2 (VGG-11').

For the sake of readability, $\bar{\mathbf{H}}$ has been divided into blocks: a weight-weight block $\bar{\mathbf{H}}_{WW}$, a bias-bias block $\bar{\mathbf{H}}_{BB}$, and a weight-bias block $\bar{\mathbf{H}}_{WB}$. They represent the interactions between the layers: for instance, $(\bar{\mathbf{H}}_{WB})_{l_1 l_2}$ represents the interaction between the tensor of weights of layer $l_1$ and the vector of biases of layer $l_2$. Naturally, $l_1$ and $l_2$ are not necessarily different.

**Results on $\bar{\mathbf{H}}$.** First, the block-diagonal approximation of the Hessian is indeed very rough, while the block-diagonal approximation of the inverse Hessian seems to be more reasonable (at least in these setups), which has already been shown by Martens and Grosse (2015). Second, long-range interactions between layers seem to exist, both at initialization and after several epochs. For instance, in LeNet, all the layers (except the first one) seem to interact together at initialization (Fig. 1). In the matrix $\bar{\mathbf{H}}^{-1}$ computed on VGG, the block of the last 3 layers interact strongly and the block of the last 6 layers also interact, but a little less.

According to these observations, a neural networks should also be regarded as a whole, in which layers can hardly be studied independently from each other.

As far as we know, this result is the first scalable illustration of interactions between layers far from each other, based on second-order information.

<figure><img src="image_2.png" /><figcaption>Figure 2: VGG-11' trained by SGD on CIFAR-10. Submatrices of $\bar{\mathbf{H}}$ (first row) and $\bar{\mathbf{H}}^{-1}$ (second row), where focus is on interactions: weight-weight, weight-bias, bias-bias of the different layers, at initialization and before the 4th epoch (best validation loss).</figcaption></figure>

**Results on $\eta_*$**. The evolution of the learning rates $\eta_*$ computed according to (2) in LeNet and VGG has been reported in Figure 3. First, the learning rates computed for the biases are larger than for the weights. Second, even when considering only the weights, the computed $\eta_*$ can differ by several order of magnitude. Finally, the first two layers of LeNet (which are convolutional) have smaller $\eta_*$ than the last three (which are fully-connected). Conversely, in VGG, the weights of the last (convolutional) layers have a smaller $\eta_*$ than in the first ones.

## 5.2 Training experiments

In this section, we show a proof-of-concept of the optimization method 4.1, on simple vision tasks and with medium-sized neural networks. All the implementation details are available in Appendix F. Notably, we have introduced a damping factor $\lambda_1$, which leads to the following modification of the training step (1):

$$ \theta_{t+1} = \theta_t - \lambda_1 \mathbf{U}_t \mathbf{I}_{P:S} \eta_t^* $$

**Setup.** We consider 4 image classification setups:

<footnote>$^{6}$VGG-11' is a variant of VGG-11 with only one fully-connected layer at the end, instead of 3.</footnote>
and Grosse, 2015) approximation of the Hessian or its inverse.

Since a summary $\bar{\mathbf{H}}$ of the Hessian and its inverse $\bar{\mathbf{H}}^{-1}$ are available and all their off-diagonal coefficients have been computed and kept, one can to check if these coefficient are indeed negligible.

**Setup.** We have trained LeNet-5 and VGG-11$^{6}$ on CIFAR-10 by SGD with momentum. Before each epoch, we compute the full-batch gradient, denoted by $\mathbf{u}$, which we use as a direction to compute $\bar{\mathbf{H}}$, again in full-batch. We report submatrices of $\bar{\mathbf{H}}$ and $\bar{\mathbf{H}}^{-1}$ at initialization and at the epoch where the validation loss is the best in Figure 1 (LeNet) and Figure 2 (VGG-11').

For the sake of readability, $\bar{\mathbf{H}}$ has been divided into blocks: a weight-weight block $\bar{\mathbf{H}}_{WW}$, a bias-bias block $\bar{\mathbf{H}}_{BB}$, and a weight-bias block $\bar{\mathbf{H}}_{WB}$. They represent the interactions between the layers: for instance, $(\bar{\mathbf{H}}_{WB})_{l_1 l_2}$ represents the interaction between the tensor of weights of layer $l_1$ and the vector of biases of layer $l_2$. Naturally, $l_1$ and $l_2$ are not necessarily different.

**Results on $\bar{\mathbf{H}}$.** First, the block-diagonal approximation of the Hessian is indeed very rough, while the block-diagonal approximation of the inverse Hessian seems to be more reasonable (at least in these setups), which has already been shown by Martens and Grosse (2015). Second, long-range interactions between layers seem to exist, both at initialization and after several epochs. For instance, in LeNet, all the layers (except the first one) seem to interact together at initialization (Fig. 1). In the matrix $\bar{\mathbf{H}}^{-1}$ computed on VGG, the block of the last 3 layers interact strongly and the block of the last 6 layers also interact, but a little less.

According to these observations, a neural networks should also be regarded as a whole, in which layers can hardly be studied independently from each other.

As far as we know, this result is the first scalable illustration of interactions between layers far from each other, based on second-order information.

<figure><img src="image_4.png" /><figcaption>Figure 2: VGG-11' trained by SGD on CIFAR-10. Submatrices of $\bar{\mathbf{H}}$ (first row) and $\bar{\mathbf{H}}^{-1}$ (second row), where focus is on interactions: weight-weight, weight-bias, bias-bias of the different layers, at initialization and before the 4th epoch (best validation loss).</figcaption></figure>

**Results on $\eta_*$**. The evolution of the learning rates $\eta_*$ computed according to (2) in LeNet and VGG has been reported in Figure 3. First, the learning rates computed for the biases are larger than for the weights. Second, even when considering only the weights, the computed $\eta_*$ can differ by several order of magnitude. Finally, the first two layers of LeNet (which are convolutional) have smaller $\eta_*$ than the last three (which are fully-connected). Conversely, in VGG, the weights of the last (convolutional) layers have a smaller $\eta_*$ than in the first ones.

## 5.2 Training experiments

In this section, we show a proof-of-concept of the optimization method 4.1, on simple vision tasks and with medium-sized neural networks. All the implementation details are available in Appendix F. Notably, we have introduced a damping factor $\lambda_1$, which leads to the following modification of the training step (1):

$$ \theta_{t+1} = \theta_t - \lambda_1 \mathbf{U}_t \mathbf{I}_{P:S} \eta_t^* $$

**Setup.** We consider 4 image classification setups:

<footnote>$^{6}$VGG-11' is a variant of VGG-11 with only one fully-connected layer at the end, instead of 3.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 3: Setup: LeNet, VGG-11' trained by SGD on CIFAR-10.</figcaption></figure>

Learning rates $\eta_*$ computed according to (2), specific to each tensor of weights and tensor of biases of each layer. For each epoch $k \in \{10, 30, 50, 70, 90\}$, the reported value has been averaged over the epochs $[k-10, k+9]$ to remove the noise.

* **MLP**: multilayer perceptron trained on MNIST with layers of sizes 1024, 200, 100, 10, with tanh activation function;

* **LeNet**: LeNet-5 (LeCun et al., 1998) model trained on CIFAR-10 with 2 convolutional layers of sizes 6, 16, and 3 fully connected layers of sizes 120, 84, 10;

* **VGG**: VGG-11' trained on CIFAR-10. VGG-11' is a variant of VGG-11 (Simonyan and Zisserman, 2014) with only one fully-connected layer at the end, instead of 3, with ELU activation function (Clevert et al., 2015), without batch-norm;

* **BigMLP**: multilayer perceptron trained on CIFAR-10, with 20 layers of size 1024 and one classification layer of size 10, with ELU activation function.

And we have tested 3 optimization methods:

* **Adam**: the best learning rate has been selected by grid-search;

* **K-FAC**: the best learning rate and damping have been selected by grid-search;

* **NewtonSummary** (ours): the best $\lambda_1$ and $\lambda_{int}$ have been selected by grid search.

**Results.** The evolution of the training loss has been plotted in Figure 4 for each of the 3 optimization methods, for 5 different seeds. In each series of experiments, training is successful, while being slow or unstable at some points. Anyway, the minimum training loss achieved by Method 4.1 (NewtonSummary) is comparable to the minimum training loss achieved by K-KAC or Adam in all the series except with BigMLP, whose training is extremely slow.

Some runs have been stopped early due to failures caused by a very large step size (large $\eta_*$). Actually, we did not use any safeguard, such as a regularization term $\lambda I$ added to $\mathbf{H}$, or clipping the learning rates exceeding a given value. This choice has been made in order to maintain the property of invariance by reparameterization, and to avoid increasing the number of hyperparameters to tune.

<figure><img src="image_4.png" /><figcaption>Figure 4: Training curves in different setups. The reported loss is the negative log-likelihood computed on the training set.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Setup: LeNet, VGG-11' trained by SGD on CIFAR-10.</figcaption></figure>

Learning rates $\eta_*$ computed according to (2), specific to each tensor of weights and tensor of biases of each layer. For each epoch $k \in \{10, 30, 50, 70, 90\}$, the reported value has been averaged over the epochs $[k-10, k+9]$ to remove the noise.

* **MLP**: multilayer perceptron trained on MNIST with layers of sizes 1024, 200, 100, 10, with tanh activation function;

* **LeNet**: LeNet-5 (LeCun et al., 1998) model trained on CIFAR-10 with 2 convolutional layers of sizes 6, 16, and 3 fully connected layers of sizes 120, 84, 10;

* **VGG**: VGG-11' trained on CIFAR-10. VGG-11' is a variant of VGG-11 (Simonyan and Zisserman, 2014) with only one fully-connected layer at the end, instead of 3, with ELU activation function (Clevert et al., 2015), without batch-norm;

* **BigMLP**: multilayer perceptron trained on CIFAR-10, with 20 layers of size 1024 and one classification layer of size 10, with ELU activation function.

And we have tested 3 optimization methods:

* **Adam**: the best learning rate has been selected by grid-search;

* **K-FAC**: the best learning rate and damping have been selected by grid-search;

* **NewtonSummary** (ours): the best $\lambda_1$ and $\lambda_{int}$ have been selected by grid search.

**Results.** The evolution of the training loss has been plotted in Figure 4 for each of the 3 optimization methods, for 5 different seeds. In each series of experiments, training is successful, while being slow or unstable at some points. Anyway, the minimum training loss achieved by Method 4.1 (NewtonSummary) is comparable to the minimum training loss achieved by K-KAC or Adam in all the series except with BigMLP, whose training is extremely slow.

Some runs have been stopped early due to failures caused by a very large step size (large $\eta_*$). Actually, we did not use any safeguard, such as a regularization term $\lambda I$ added to $\mathbf{H}$, or clipping the learning rates exceeding a given value. This choice has been made in order to maintain the property of invariance by reparameterization, and to avoid increasing the number of hyperparameters to tune.

<figure><img src="image_6.png" /><figcaption>Figure 4: Training curves in different setups. The reported loss is the negative log-likelihood computed on the training set.</figcaption></figure>
**Extension to very large models.** Since the matrix $\tilde{H}$ can be computed numerically as long as $S$ remains relatively small, then, for very large models, this method may become unpractical. However, Method 4.1 is flexible enough to be adapted to such models: one can regroup tensors “of the same kind” in order to build a coarser partition of the parameters, and consequently obtain a small $S$, which is exactly what is needed to compute $\tilde{H}$ and invert it. The main difficulty would then be to find a good partition of the parameters, by grouping all the tensors that “look alike”. We provide an example in Appendix G, with a very deep multilayer perceptron.

# 6 Discussion

**Convergence rate.** Method 4.1 does not come with any convergence rate. Given the convergence rates of Newton's method and Cauchy's steepest descent, we may expect to find some in-between convergence rates. Since Cauchy's steepest method is vulnerable to a highly anisotropic Hessian, it would be valuable to know how much this weakness is overcome with our method.

**Practicality.** Despite the interesting properties of Method 4.1 (scalability, invariance by reparameterization, evaluation of long-range interactions between layers), we have proposed nothing more than a proof-of-concept. As showed in the experimental section, this method is subject to instabilities during training, which is certainly not surprising for a second-order method, but it is not acceptable for the end user. Therefore, some additional tricks should be added to improve the stability of training, which is a common practice, but comes usually with additional hyper-parameters to tune. Besides, this method has to be improved to reduce the duration of each epoch, which is longer than with K-FAC.

# 7 Acknowledgments

This work was granted access to the HPC resources of IDRIS under the allocation 2023-AD011013762R1
**Extension to very large models.** Since the matrix $\tilde{H}$ can be computed numerically as long as $S$ remains relatively small, then, for very large models, this method may become unpractical. However, Method 4.1 is flexible enough to be adapted to such models: one can regroup tensors “of the same kind” in order to build a coarser partition of the parameters, and consequently obtain a small $S$, which is exactly what is needed to compute $\tilde{H}$ and invert it. The main difficulty would then be to find a good partition of the parameters, by grouping all the tensors that “look alike”. We provide an example in Appendix G, with a very deep multilayer perceptron.

made by GENCI.

# 6 Discussion

**Convergence rate.** Method 4.1 does not come with any convergence rate. Given the convergence rates of Newton's method and Cauchy's steepest descent, we may expect to find some in-between convergence rates. Since Cauchy's steepest method is vulnerable to a highly anisotropic Hessian, it would be valuable to know how much this weakness is overcome with our method.

**Practicality.** Despite the interesting properties of Method 4.1 (scalability, invariance by reparameterization, evaluation of long-range interactions between layers), we have proposed nothing more than a proof-of-concept. As showed in the experimental section, this method is subject to instabilities during training, which is certainly not surprising for a second-order method, but it is not acceptable for the end user. Therefore, some additional tricks should be added to improve the stability of training, which is a common practice, but comes usually with additional hyper-parameters to tune. Besides, this method has to be improved to reduce the duration of each epoch, which is longer than with K-FAC.

# 7 Acknowledgments

This work was granted access to the HPC resources of IDRIS under the allocation 2023-AD011013762R1
## References

Amari, S.-I. (1998). Natural gradient works efficiently in learning. *Neural computation*, 10(2):251–276.

Cauchy, A.-L. (1847). Méthode générale pour la résolution des systèmes d’équations simultanées. *Comptes rendus hebdomadaires des séances de l’Académie des sciences, Paris*, 25:536–538.

Clevert, D.-A., Unterthiner, T., and Hochreiter, S. (2015). Fast and accurate deep network learning by exponential linear units (elus). *arXiv preprint arXiv:1511.07289*.

Dangel, F. J. (2023). *Backpropagation beyond the gradient*. PhD thesis, Universität Tübingen.

Gill, P. E., Murray, W., and Wright, M. H. (1981). *Practical optimization*. Academic Press, San Diego.

Jacot, A., Gabriel, F., and Hongler, C. (2018). Neural tangent kernel: Convergence and generalization in neural networks. *Advances in Neural Information Processing Systems*, 31.

Kornblith, S., Norouzi, M., Lee, H., and Hinton, G. (2019). Similarity of neural network representations revisited. In *International conference on machine learning*, pages 3519–3529. PMLR.

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324.

Lu, Y., Harandi, M., Hartley, R., and Pascanu, R. (2018). Block mean approximation for efficient second order optimization. *arXiv preprint arXiv:1804.05484*.

Luenberger, D. G. and Ye, Y. (2008). *Linear and Nonlinear Programming*. Springer, fourth edition.

Martens, J. and Grosse, R. (2015). Optimizing neural networks with Kronecker-factored approximate curvature. In *International conference on machine learning*, pages 2408–2417. PMLR.

Nesterov, Y. (2003). *Introductory lectures on convex optimization: A basic course*, volume 87. Springer Science & Business Media.

Nesterov, Y. and Polyak, B. T. (2006). Cubic regularization of Newton method and its global performance. *Mathematical Programming*, 108(1):177–205.

Nocedal, J. and Wright, S. J. (1999). *Numerical optimization*. Springer.

Ollivier, Y. (2015). Riemannian metrics for neural networks I: feedforward networks. *arXiv preprint arXiv:1303.0818*.

Pearlmutter, B. A. (1994). Fast exact multiplication by the Hessian. *Neural computation*, 6(1):147–160.

Sagun, L., Evci, U., Guney, V. U., Dauphin, Y., and Bottou, L. (2018). Empirical analysis of the Hessian of over-parametrized neural networks. In *International Conference on Learning Representations*.

Simonyan, K. and Zisserman, A. (2014). Very deep convolutional networks for large-scale image recognition. *arXiv preprint arXiv:1409.1556*.

Wang, Y.-J. and Lin, C.-T. (1998). A second-order learning algorithm for multilayer networks based on block Hessian matrix. *Neural Networks*, 11(9):1607–1622.

Yang, G. and Hu, E. J. (2021). Tensor programs iv: Feature learning in infinite-width neural networks. In *International Conference on Machine Learning*, pages 11727–11737. PMLR.

Yuan, R., Lazaric, A., and Gower, R. M. (2022). Sketched Newton-Raphson. *SIAM Journal on Optimization*, 32(3):1555–1583.

Zhang, C., Bengio, S., and Singer, Y. (2022). Are all layers created equal? *The Journal of Machine Learning Research*, 23(1):2930–2957.
## References

Amari, S.-I. (1998). Natural gradient works efficiently in learning. *Neural computation*, 10(2):251–276.

Cauchy, A.-L. (1847). Méthode générale pour la résolution des systèmes d’équations simultanées. *Comptes rendus hebdomadaires des séances de l’Académie des sciences, Paris*, 25:536–538.

Clevert, D.-A., Unterthiner, T., and Hochreiter, S. (2015). Fast and accurate deep network learning by exponential linear units (elus). *arXiv preprint arXiv:1511.07289*.

Dangel, F. J. (2023). *Backpropagation beyond the gradient*. PhD thesis, Universität Tübingen.

Gill, P. E., Murray, W., and Wright, M. H. (1981). *Practical optimization*. Academic Press, San Diego.

Jacot, A., Gabriel, F., and Hongler, C. (2018). Neural tangent kernel: Convergence and generalization in neural networks. *Advances in Neural Information Processing Systems*, 31.

Kornblith, S., Norouzi, M., Lee, H., and Hinton, G. (2019). Similarity of neural network representations revisited. In *International conference on machine learning*, pages 3519–3529. PMLR.

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324.

Lu, Y., Harandi, M., Hartley, R., and Pascanu, R. (2018). Block mean approximation for efficient second order optimization. *arXiv preprint arXiv:1804.05484*.

Luenberger, D. G. and Ye, Y. (2008). *Linear and Nonlinear Programming*. Springer, fourth edition.

Martens, J. and Grosse, R. (2015). Optimizing neural networks with Kronecker-factored approximate curvature. In *International conference on machine learning*, pages 2408–2417. PMLR.

Nesterov, Y. (2003). *Introductory lectures on convex optimization: A basic course*, volume 87. Springer Science & Business Media.

Nesterov, Y. and Polyak, B. T. (2006). Cubic regularization of Newton method and its global performance. *Mathematical Programming*, 108(1):177–205.

Nocedal, J. and Wright, S. J. (1999). *Numerical optimization*. Springer.

Ollivier, Y. (2015). Riemannian metrics for neural networks I: feedforward networks. *arXiv preprint arXiv:1303.0818*.

Pearlmutter, B. A. (1994). Fast exact multiplication by the Hessian. *Neural computation*, 6(1):147–160.

Sagun, L., Evci, U., Guney, V. U., Dauphin, Y., and Bottou, L. (2018). Empirical analysis of the Hessian of over-parametrized neural networks. In *International Conference on Learning Representations*.

Simonyan, K. and Zisserman, A. (2014). Very deep convolutional networks for large-scale image recognition. *arXiv preprint arXiv:1409.1556*.

Wang, Y.-J. and Lin, C.-T. (1998). A second-order learning algorithm for multilayer networks based on block Hessian matrix. *Neural Networks*, 11(9):1607–1622.

Yang, G. and Hu, E. J. (2021). Tensor programs iv: Feature learning in infinite-width neural networks. In *International Conference on Machine Learning*, pages 11727–11737. PMLR.

Yuan, R., Lazaric, A., and Gower, R. M. (2022). Sketched Newton-Raphson. *SIAM Journal on Optimization*, 32(3):1555–1583.

Zhang, C., Bengio, S., and Singer, Y. (2022). Are all layers created equal? *The Journal of Machine Learning Research*, 23(1):2930–2957.
## A Fast computation of the terms of the Taylor expansion

In this appendix, we show how to use the trick of Pearlmutter (1994) to compute the terms of the Taylor expansion of $\mathcal{L}$. We recall that we want to compute:

$$ \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u}) := \frac{\partial^d \mathcal{L}}{\partial \boldsymbol{\theta}^d}(\boldsymbol{\theta})[\mathbf{u}, \dots, \mathbf{u}] \in \mathbb{R}. $$

We use the following recursion formula:

$$ \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d+1)}(\mathbf{u}) = \left( \frac{\partial \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})}{\partial \boldsymbol{\theta}} \right)^T \mathbf{u}. $$

Therefore, at each step $d$, we only have to compute the gradient of a scalar $\tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})$ according to $\boldsymbol{\theta} \in \mathbb{R}^P$, and compute a dot product in the space $\boldsymbol{\theta} \in \mathbb{R}^P$. So, computing $\tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})$ has a complexity proportional to $d \times P$, and does not require the computation of the full tensor $\frac{\partial^d \mathcal{L}}{\partial \boldsymbol{\theta}^d}(\boldsymbol{\theta}) \in \mathbb{R}^{P^d}$.

## B Derivation of the second-order method

We consider an update of $\boldsymbol{\theta}$ with one learning rate $\eta_s$ for each subset $\mathcal{I}_s$ of parameters. Let $\mathbf{I}_{S:P} \in \mathbb{R}^{S \times P}$ be the partition matrix, verifying $(\mathbf{I}_{S:P})_{sp} = 1$ if $p \in \mathcal{I}_s$ and 0 otherwise, and $\mathbf{I}_{P:S} := \mathbf{I}_{S:P}^T$. We consider an update based on a given direction $\mathbf{u}_t$ and we define $\mathbf{U}_t := \text{Diag}(\mathbf{u}_t)$:

$$ \boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}, $$

where $\boldsymbol{\eta} = (\eta_1, \dots, \eta_S) \in \mathbb{R}^S$.

The second-order approximation of $\mathcal{L}$ gives:

$$
\begin{align*}
\mathcal{L}(\boldsymbol{\theta}_{t+1}) &= \mathcal{L}(\boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \frac{\partial \mathcal{L}}{\partial \boldsymbol{\theta}}(\boldsymbol{\theta}_t) \\
&\quad + \frac{1}{2} \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{\theta}^2}(\boldsymbol{\theta}_t) \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t \\
&\quad + \frac{1}{2} \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \bar{\mathbf{g}}_t + \frac{1}{2} \boldsymbol{\eta}^T \bar{\mathbf{H}}_t \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2),
\end{align*}
$$

where:

$$ \bar{\mathbf{g}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t \quad \in \mathbb{R}^S, $$

$$ \bar{\mathbf{H}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \quad \in \mathbb{R}^{S \times S}. $$

Now, we omit the $o(\|\boldsymbol{\eta}\|^2)$ term and we want to minimize according to $\boldsymbol{\eta}$ the variation of the loss:

$$ \begin{align*}
\mathcal{L}(\boldsymbol{\theta}_{t+1}) - \mathcal{L}(\boldsymbol{\theta}_t) \\
\approx \boldsymbol{\Delta}_2(\boldsymbol{\eta}) := \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \bar{\mathbf{g}}_t + \frac{1}{2} \boldsymbol{\eta}^T \bar{\mathbf{H}}_t \boldsymbol{\eta}.
\end{align*} $$

We have:

$$ \frac{\partial \boldsymbol{\Delta}_2}{\partial \boldsymbol{\eta}} = -\bar{\mathbf{g}}_t + \bar{\mathbf{H}}_t \boldsymbol{\eta}, $$

which is zero if, and only if:

$$ \bar{\mathbf{g}}_t = \bar{\mathbf{H}}_t \boldsymbol{\eta}. $$

If this linear system can be inverted, one may choose:

$$ \boldsymbol{\eta} = \boldsymbol{\eta}_t^* := \bar{\mathbf{H}}_t^{-1} \bar{\mathbf{g}}_t. $$

## C Link with Cauchy's steepest descent and Newton's method

**Cauchy's steepest descent.** Let us consider the trivial partition: $S = 1$, $\mathcal{I}_1 = \{1, \dots, P\}$. So, $\mathbf{I}_{S:P} = (1, \dots, 1)^T$. Therefore, the training step is:

$$
\begin{align*}
\boldsymbol{\theta}_{t+1} &:= \boldsymbol{\theta}_t - \mathbf{G}_t \mathbf{1}_S (\mathbf{1}_S^T \mathbf{G}_t \mathbf{H}_t \mathbf{G}_t \mathbf{1}_S)^{-1} \mathbf{1}_S^T \mathbf{G}_t \mathbf{g}_t \\
&= \boldsymbol{\theta}_t - \mathbf{g}_t \frac{\mathbf{g}_t^T \mathbf{g}_t}{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t},
\end{align*}
$$
## A Fast computation of the terms of the Taylor expansion

In this appendix, we show how to use the trick of Pearlmutter (1994) to compute the terms of the Taylor expansion of $\mathcal{L}$. We recall that we want to compute:

$$ \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u}) := \frac{\partial^d \mathcal{L}}{\partial \boldsymbol{\theta}^d}(\boldsymbol{\theta})[\mathbf{u}, \dots, \mathbf{u}] \in \mathbb{R}. $$

We use the following recursion formula:

$$ \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d+1)}(\mathbf{u}) = \left( \frac{\partial \tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})}{\partial \boldsymbol{\theta}} \right)^T \mathbf{u}. $$

Therefore, at each step $d$, we only have to compute the gradient of a scalar $\tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})$ according to $\boldsymbol{\theta} \in \mathbb{R}^P$, and compute a dot product in the space $\boldsymbol{\theta} \in \mathbb{R}^P$. So, computing $\tilde{\mathbf{D}}_{\boldsymbol{\theta}}^{(d)}(\mathbf{u})$ has a complexity proportional to $d \times P$, and does not require the computation of the full tensor $\frac{\partial^d \mathcal{L}}{\partial \boldsymbol{\theta}^d}(\boldsymbol{\theta}) \in \mathbb{R}^{P^d}$.

## B Derivation of the second-order method

We consider an update of $\boldsymbol{\theta}$ with one learning rate $\eta_s$ for each subset $\mathcal{I}_s$ of parameters. Let $\mathbf{I}_{S:P} \in \mathbb{R}^{S \times P}$ be the partition matrix, verifying $(\mathbf{I}_{S:P})_{sp} = 1$ if $p \in \mathcal{I}_s$ and 0 otherwise, and $\mathbf{I}_{P:S} := \mathbf{I}_{S:P}^T$. We consider an update based on a given direction $\mathbf{u}_t$ and we define $\mathbf{U}_t := \text{Diag}(\mathbf{u}_t)$:

$$ \boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}, $$

where $\boldsymbol{\eta} = (\eta_1, \dots, \eta_S) \in \mathbb{R}^S$.

The second-order approximation of $\mathcal{L}$ gives:

$$
\begin{align*}
\mathcal{L}(\boldsymbol{\theta}_{t+1}) &= \mathcal{L}(\boldsymbol{\theta}_t - \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta}) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \frac{\partial \mathcal{L}}{\partial \boldsymbol{\theta}}(\boldsymbol{\theta}_t) \\
&\quad + \frac{1}{2} \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{\theta}^2}(\boldsymbol{\theta}_t) \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t \\
&\quad + \frac{1}{2} \boldsymbol{\eta}^T \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2) \\
&= \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \bar{\mathbf{g}}_t + \frac{1}{2} \boldsymbol{\eta}^T \bar{\mathbf{H}}_t \boldsymbol{\eta} + o(\|\boldsymbol{\eta}\|^2),
\end{align*}
$$

where:

$$ \bar{\mathbf{g}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{g}_t \quad \in \mathbb{R}^S, $$

$$ \bar{\mathbf{H}}_t := \mathbf{I}_{S:P} \mathbf{U}_t \mathbf{H}_t \mathbf{U}_t \mathbf{I}_{P:S} \quad \in \mathbb{R}^{S \times S}. $$

Now, we omit the $o(\|\boldsymbol{\eta}\|^2)$ term and we want to minimize according to $\boldsymbol{\eta}$ the variation of the loss:

$$
\begin{align*}
\mathcal{L}(\boldsymbol{\theta}_{t+1}) - \mathcal{L}(\boldsymbol{\theta}_t) \\
\approx \boldsymbol{\Delta}_2(\boldsymbol{\eta}) := \mathcal{L}(\boldsymbol{\theta}_t) - \boldsymbol{\eta}^T \bar{\mathbf{g}}_t + \frac{1}{2} \boldsymbol{\eta}^T \bar{\mathbf{H}}_t \boldsymbol{\eta}.
\end{align*}
$$

We have:

$$ \frac{\partial \boldsymbol{\Delta}_2}{\partial \boldsymbol{\eta}} = -\bar{\mathbf{g}}_t + \bar{\mathbf{H}}_t \boldsymbol{\eta}, $$

which is zero if, and only if:

$$ \bar{\mathbf{g}}_t = \bar{\mathbf{H}}_t \boldsymbol{\eta}. $$

If this linear system can be inverted, one may choose:

$$ \boldsymbol{\eta} = \boldsymbol{\eta}_t^* := \bar{\mathbf{H}}_t^{-1} \bar{\mathbf{g}}_t. $$

## C Link with Cauchy's steepest descent and Newton's method

**Cauchy's steepest descent.** Let us consider the trivial partition: $S = 1$, $\mathcal{I}_1 = \{1, \dots, P\}$. So, $\mathbf{I}_{S:P} = (1, \dots, 1)^T$. Therefore, the training step is:

$$
\begin{align*}
\boldsymbol{\theta}_{t+1} &:= \boldsymbol{\theta}_t - \mathbf{G}_t \mathbf{1}_S (\mathbf{1}_S^T \mathbf{G}_t \mathbf{H}_t \mathbf{G}_t \mathbf{1}_S)^{-1} \mathbf{1}_S^T \mathbf{G}_t \mathbf{g}_t \\
&= \boldsymbol{\theta}_t - \mathbf{g}_t \frac{\mathbf{g}_t^T \mathbf{g}_t}{\mathbf{g}_t^T \mathbf{H}_t \mathbf{g}_t},
\end{align*}
$$
since $G_t 1_S = g_t$. We recover Cauchy's steepest descent.

**Newton's method.** Since we aim to recover Newton's method, we assume that the Hessian $H_t$ is positive definite. Let us consider the discrete partition: $S = P$, $\mathcal{I}_s = \{s\}$. So, $I_{S:P} = I_P$, the identity matrix of $\mathbb{R}^{P \times P}$. Therefore, the training step is:

$$ \theta_{t+1} := \theta_t - G_t (G_t H_t G_t)^{-1} G_t g_t. $$

To perform the training step, we have to find $\mathbf{x} \in \mathbb{R}^P$ such that: $(G_t H_t G_t)^{-1} G_t g_t = \mathbf{x}$. That is, solve the linear system $G_t H_t G_t \mathbf{x} = G_t g_t$. In the case where all the coordinates of the gradient $g_t$ are nonzero, we can write:

$$ \mathbf{x} = G_t^{-1} H_t^{-1} G_t^{-1} G_t g_t = G_t^{-1} H_t^{-1} g_t, $$

so the training step becomes:

$$ \theta_{t+1} := \theta_t - G_t \mathbf{x} = \theta_t - H_t^{-1} g_t, $$

which corresponds to Newton's method.

## D Anisotropic Nesterov cubic regularization

Let $\mathbf{D}$ be a diagonal matrix whose diagonal coefficients are all strictly positive: $\mathbf{D} = \text{Diag}(d_1, \dots, d_S)$, with $d_i > 0$ for all $i$.

We want to minimize the function:

$$ T(\eta) := -\eta^T \bar{\mathbf{g}} + \frac{1}{2} \eta \bar{\mathbf{H}} \eta + \frac{\lambda_{\text{int}}}{6} \|\mathbf{D}\eta\|^3. $$

The function $T$ is strictly convex if, and only if, $\bar{\mathbf{H}}$ is positive definite. Moreover, $T$ is differentiable twice and has at least one global minimum $\eta_*$, so $\frac{\partial T}{\partial \eta}(\eta_*) = 0$. Therefore, we first look for the solutions of the equation $\frac{\partial T}{\partial \eta}(\eta) = 0$.

We have:

$$ \begin{aligned} \frac{\partial T}{\partial \eta}(\eta) &= -\bar{\mathbf{g}} + \bar{\mathbf{H}}\eta + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2 \eta \\ &= -\bar{\mathbf{g}} + \left(\bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2\right) \eta, \end{aligned} $$

which is equal to zero if, and only if:

$$ \bar{\mathbf{g}} = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2 \right) \eta. \quad (3) $$

Let $\eta' := \mathbf{D}\eta$. Eqn. (3) is then equivalent to:

$$ \begin{aligned} \bar{\mathbf{g}} &= \left( \bar{\mathbf{H}} \mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2} \|\eta'\| \mathbf{D} \right) \eta' \\ &= \frac{\lambda_{\text{int}}}{2} \mathbf{D} \left( \frac{2}{\lambda_{\text{int}}} \mathbf{D}^{-1} \bar{\mathbf{H}} \mathbf{D}^{-1} + \|\eta'\| \mathbf{I} \right) \eta' \end{aligned} $$

Let $\mathbf{K} := \frac{2}{\lambda_{\text{int}}} \mathbf{D}^{-1} \bar{\mathbf{H}} \mathbf{D}^{-1}$. We want to solve:

$$ \bar{\mathbf{g}} = \frac{\lambda_{\text{int}}}{2} \mathbf{D} (\mathbf{K} + \|\eta'\| \mathbf{I}) \eta' \quad (4) $$

Since $\mathbf{K}$ is positive definite if, and only if, $\bar{\mathbf{H}}$ is positive definite, we consider the following cases.

**Case 1: $\bar{\mathbf{H}}$ is positive definite.** In this case, Eqn. (4) is equivalent to:

$$ \eta' = \frac{2}{\lambda_{\text{int}}} (\mathbf{K} + \|\eta'\| \mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}}. $$

Now, let $r = \|\eta'\|$. We want to solve:

$$ r = \frac{2}{\lambda_{\text{int}}} \left\| (\mathbf{K} + r\mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}} \right\|. \quad (5) $$

Trivially: $\eta$ solution of (3) $\Rightarrow$ $\mathbf{D}\eta$ solution of (4) $\Rightarrow$ $\|\mathbf{D}\eta\|$ solution of (5). Reciprocally: $r$ solution of (5) $\Rightarrow$ $\eta' := (\bar{\mathbf{H}}\mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2} r\mathbf{D})^{-1} \bar{\mathbf{g}}$ solution of (4) $\Rightarrow$ $\mathbf{D}^{-1}\eta'$ solution of (3).

Therefore, in order to find the unique global minimum of $T$, it is sufficient to solve Eqn. (5). This is doable numerically.

**Case 2: $\bar{\mathbf{H}}$ is not positive definite.** We follow the procedure proposed in (Nesterov and Polyak, 2006, Section 5). Let $\lambda_{\min}$ be the minimum eigenvalue of $\mathbf{K}$. So, $\lambda_{\min} \le 0$. Following Nesterov and Polyak (2006), we look for the unique $\eta'$ belonging to $\mathcal{C} := \{\eta' \in \mathbb{R}^S : \|\eta'\| > |\lambda_{\min}|\}$, which is also the solution of maximum norm of Eqn. (4). Conditionally to $\eta' \in \mathcal{C}$, $(\mathbf{K} + \|\eta'\|\mathbf{I})$ is invertible. So we only need to solve:

$$ r > |\lambda_{\min}|: \quad r = \frac{2}{\lambda_{\text{int}}} \left\| (\mathbf{K} + r\mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}} \right\|, \quad (6) $$
since $G_t 1_S = g_t$. We recover Cauchy's steepest descent.

**Newton's method.** Since we aim to recover Newton's method, we assume that the Hessian $H_t$ is positive definite. Let us consider the discrete partition: $S = P$, $\mathcal{I}_s = \{s\}$. So, $\mathbf{I}_{S:P} = \mathbf{I}_P$, the identity matrix of $\mathbb{R}^{P \times P}$. Therefore, the training step is:

$$ \theta_{t+1} := \theta_t - G_t (G_t H_t G_t)^{-1} G_t g_t. $$

To perform the training step, we have to find $\mathbf{x} \in \mathbb{R}^P$ such that: $(G_t H_t G_t)^{-1} G_t g_t = \mathbf{x}$. That is, solve the linear system $G_t H_t G_t \mathbf{x} = G_t g_t$. In the case where all the coordinates of the gradient $g_t$ are nonzero, we can write:

$$ \mathbf{x} = G_t^{-1} H_t^{-1} G_t^{-1} G_t g_t = G_t^{-1} H_t^{-1} g_t, $$

so the training step becomes:

$$ \theta_{t+1} := \theta_t - G_t \mathbf{x} = \theta_t - H_t^{-1} g_t, $$

which corresponds to Newton's method.

# D Anisotropic Nesterov cubic regularization

Let $\mathbf{D}$ be a diagonal matrix whose diagonal coefficients are all strictly positive: $\mathbf{D} = \text{Diag}(d_1, \dots, d_S)$, with $d_i > 0$ for all $i$.

We want to minimize the function:

$$ T(\eta) := -\eta^T \bar{\mathbf{g}} + \frac{1}{2} \eta \bar{\mathbf{H}} \eta + \frac{\lambda_{\text{int}}}{6} \|\mathbf{D}\eta\|^3. $$

The function $T$ is strictly convex if, and only if, $\bar{\mathbf{H}}$ is positive definite. Moreover, $T$ is differentiable twice and has at least one global minimum $\eta_*$, so $\frac{\partial T}{\partial \eta}(\eta_*) = 0$. Therefore, we first look for the solutions of the equation $\frac{\partial T}{\partial \eta}(\eta) = 0$.

We have:

$$ \begin{aligned} \frac{\partial T}{\partial \eta}(\eta) &= -\bar{\mathbf{g}} + \bar{\mathbf{H}}\eta + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2 \eta \\ &= -\bar{\mathbf{g}} + \left(\bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2\right) \eta, \end{aligned} $$

which is equal to zero if, and only if:

$$ \bar{\mathbf{g}} = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta\| \mathbf{D}^2 \right) \eta. \quad (3) $$

Let $\eta' := \mathbf{D}\eta$. Eqn. (3) is then equivalent to:

$$ \begin{aligned} \bar{\mathbf{g}} &= \left( \bar{\mathbf{H}} \mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2} \|\eta'\| \mathbf{D} \right) \eta' \\ &= \frac{\lambda_{\text{int}}}{2} \mathbf{D} \left( \frac{2}{\lambda_{\text{int}}} \mathbf{D}^{-1} \bar{\mathbf{H}} \mathbf{D}^{-1} + \|\eta'\| \mathbf{I} \right) \eta' \end{aligned} $$

Let $\mathbf{K} := \frac{2}{\lambda_{\text{int}}} \mathbf{D}^{-1} \bar{\mathbf{H}} \mathbf{D}^{-1}$. We want to solve:

$$ \bar{\mathbf{g}} = \frac{\lambda_{\text{int}}}{2} \mathbf{D} (\mathbf{K} + \|\eta'\| \mathbf{I}) \eta' \quad (4) $$

Since $\mathbf{K}$ is positive definite if, and only if, $\bar{\mathbf{H}}$ is positive definite, we consider the following cases.

**Case 1: $\bar{\mathbf{H}}$ is positive definite.** In this case, Eqn. (4) is equivalent to:

$$ \eta' = \frac{2}{\lambda_{\text{int}}} (\mathbf{K} + \|\eta'\| \mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}}. $$

Now, let $r = \|\eta'\|$. We want to solve:

$$ r = \frac{2}{\lambda_{\text{int}}} \left\| (\mathbf{K} + r\mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}} \right\|. \quad (5) $$

Trivially: $\eta$ solution of (3) $\Rightarrow$ $\mathbf{D}\eta$ solution of (4) $\Rightarrow$ $\|\mathbf{D}\eta\|$ solution of (5). Reciprocally: $r$ solution of (5) $\Rightarrow$ $\eta' := (\bar{\mathbf{H}}\mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2} r\mathbf{D})^{-1} \bar{\mathbf{g}}$ solution of (4) $\Rightarrow$ $\mathbf{D}^{-1}\eta'$ solution of (3).

Therefore, in order to find the unique global minimum of $T$, it is sufficient to solve Eqn. (5). This is doable numerically.

**Case 2: $\bar{\mathbf{H}}$ is not positive definite.** We follow the procedure proposed in (Nesterov and Polyak, 2006, Section 5). Let $\lambda_{\min}$ be the minimum eigenvalue of $\mathbf{K}$. So, $\lambda_{\min} \le 0$. Following Nesterov and Polyak (2006), we look for the unique $\eta'$ belonging to $\mathcal{C} := \{\eta' \in \mathbb{R}^S : \|\eta'\| > |\lambda_{\min}|\}$, which is also the solution of maximum norm of Eqn. (4). Conditionally to $\eta' \in \mathcal{C}$, $(\mathbf{K} + \|\eta'\| \mathbf{I})$ is invertible. So we only need to solve:

$$ r > |\lambda_{\min}|: \quad r = \frac{2}{\lambda_{\text{int}}} \left\| (\mathbf{K} + r\mathbf{I})^{-1} \mathbf{D}^{-1} \bar{\mathbf{g}} \right\|, \quad (6) $$
which has exactly one solution $r_*. Then, we compute
$$ \eta_* := \mathbf{D}^{-1}(\mathbf{H}\mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2}r_*\mathbf{D})^{-1}\tilde{\mathbf{g}}. $$

## E Approximate invariance by subset-wise affine reparametrization

We consider a parameter $\tilde{\theta}$ such that:

$$ \theta = \varphi(\tilde{\theta}), $$

where $\varphi$ is an invertible map, affine on each subset of parameters. Therefore, its Jacobian is: $\mathbf{J} = \text{Diag}(\alpha_1, \dots, \alpha_p)$, where, for all $1 \le s \le S$ and $1 \le p_1, p_2 \le P$, we have:

$$ p_1, p_2 \in \mathcal{I}_s \Rightarrow \alpha_{p_1} = \alpha_{p_2} =: a_s. $$

Also, let $\tilde{\mathbf{J}} = \text{Diag}(a_1, \dots, a_S)$.

We want to compare the training trajectory of $\mathcal{L}(\theta)$ and $\mathcal{L}(\varphi(\tilde{\theta}))$ when using Method 4.1. For any quantity $\mathbf{x}$ computed with the parameterization $\theta$, we denote by $\tilde{\mathbf{x}}$ its counterpart computed with the parameterization $\tilde{\theta}$.

We compute $\tilde{\eta}_*$. Equation (2) gives:

$$ \tilde{\eta}_* = \left( \tilde{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\tilde{\mathbf{D}}\tilde{\eta}_*\| \tilde{\mathbf{D}}^2 \right)^{-1} \tilde{\mathbf{g}}. \quad (7) $$

Besides:

$$ \tilde{\mathbf{H}} := \mathbf{I}_{S:P} \tilde{\mathbf{U}} \tilde{\mathbf{H}} \tilde{\mathbf{U}} \mathbf{I}_{P:S} $$

$$ \tilde{\mathbf{g}} := \mathbf{I}_{S:P} \tilde{\mathbf{U}} \tilde{\mathbf{g}} $$

To go further, we need to do an assumption about the direction $\mathbf{u}$.

**Assumption E.1.** We assume that $\mathbf{U}_t$ is computed in such a way that $\tilde{\mathbf{U}}_t = \mathbf{JU}_t$ at every step.

This assumption holds typically when $\mathbf{u}_t$ is the gradient at time step $t$. It holds also when $\mathbf{u}_t$ is a linear combination of the past gradients:

$$ \mathbf{u}_1 := \mathbf{g}_1 $$

$$ \mathbf{u}_{t+1} := \mu \mathbf{u}_t + \mu' \mathbf{g}_{t+1}, $$

which includes the momentum.

To summarize, we have:

$$ \tilde{\mathbf{U}} = \mathbf{JU}, \quad \tilde{\mathbf{H}} = \mathbf{JHJ}, \quad \tilde{\mathbf{g}} = \mathbf{Jg}, $$

So:

$$ \tilde{\mathbf{H}} = \tilde{\mathbf{J}}^2 \mathbf{I}_{S:P} \mathbf{U} \mathbf{H} \mathbf{I}_{P:S} \tilde{\mathbf{J}}^2 = \tilde{\mathbf{J}}^2 \tilde{\mathbf{H}} \tilde{\mathbf{J}}^2, $$

$$ \tilde{\mathbf{g}} = \tilde{\mathbf{J}}^2 \mathbf{I}_{S:P} \mathbf{U} \mathbf{g} = \tilde{\mathbf{J}}^2 \tilde{\mathbf{g}}, $$

since $\mathbf{J}$ and $\mathbf{U}$ are diagonal. And:

$$ \text{since: } \mathbf{D}_{ii} = \left| \left( \mathbf{D}_{\theta}^{(3)}(\mathbf{u}) \right)_{iii} \right|^{1/3}, $$

$$ \text{then: } \tilde{\mathbf{D}}_{ii} = a_i^2 \mathbf{D}_{ii}, $$

$$ \text{thus: } \tilde{\mathbf{D}} = \tilde{\mathbf{J}}^2 \mathbf{D}. $$

Thus, Eqn. (7) becomes:

$$ \tilde{\eta}_* = \left( \tilde{\mathbf{J}}^2 \tilde{\mathbf{H}} \tilde{\mathbf{J}}^2 + \frac{\lambda_{\text{int}}}{2} \| \tilde{\mathbf{J}}^2 \tilde{\mathbf{D}} \tilde{\eta}_* \| \tilde{\mathbf{J}}^4 \mathbf{D}^2 \right)^{-1} \tilde{\mathbf{J}}^2 \tilde{\mathbf{g}}, $$

which can be rewritten (since $\tilde{\mathbf{J}}$ is invertible):

$$ \tilde{\mathbf{J}}^2 \tilde{\eta}_* = \left( \tilde{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \| \tilde{\mathbf{D}} \tilde{\mathbf{J}}^2 \tilde{\eta}_* \| \tilde{\mathbf{D}}^2 \right)^{-1} \tilde{\mathbf{g}}. $$

Therefore, $\tilde{\eta}_*$ is a solution of Eqn. (2) in the parameterization $\tilde{\theta}$ if, and only if, $\tilde{\mathbf{J}}^2 \tilde{\eta}_*$ is a solution in the parameterization $\theta$. Moreover, $\|\tilde{\mathbf{D}}\tilde{\eta}_*\| = \|\tilde{\mathbf{D}}\tilde{\mathbf{J}}^2\tilde{\eta}_*\|$, so $\tilde{\eta}_*$ is the solution of maximum norm $\|\tilde{\mathbf{D}}\tilde{\eta}_*\|$ of (2) with parameterization $\tilde{\theta}$ iff $\tilde{\mathbf{J}}^2\tilde{\eta}_*$ is a the solution of maximum norm $\|\tilde{\mathbf{D}}\tilde{\mathbf{J}}^2\tilde{\eta}_*\|$ of (2) with parameterization $\theta$.

Thus, $\eta_* = \tilde{\mathbf{J}}^2\tilde{\eta}_*$, and the update step in parameterization $\tilde{\theta}$ is:

$$ \begin{aligned} \tilde{\theta}_{t+1} &= \tilde{\theta}_t - \tilde{\mathbf{U}}_t \mathbf{I}_{P:S} \tilde{\eta}_* \\ &= \tilde{\theta}_t - \tilde{\mathbf{U}}_t \mathbf{I}_{P:S} \tilde{\mathbf{J}}^{-2} \eta_* \end{aligned} $$

which can be rewritten:

$$ \mathbf{J}^{-1} \theta_{t+1} = \mathbf{J}^{-1} \theta_t - \mathbf{U} \mathbf{J} \mathbf{I}_{P:S} \tilde{\mathbf{J}}^{-2} \eta_*, \quad (8) $$

since $\varphi$ is an affine function with factor $\mathbf{J}$. Finally, En. (8) boils down to:

$$ \theta_{t+1} = \theta_t - \mathbf{U} \mathbf{I}_{P:S} \eta_*, $$

which is exactly Method 4.1 in parameterization $\theta$.
which has exactly one solution $r_*. Then, we compute
$\eta_* := \mathbf{D}^{-1}(\mathbf{H}\mathbf{D}^{-1} + \frac{\lambda_{\text{int}}}{2}r_*\mathbf{D})^{-1}\tilde{\mathbf{g}}$.

E Approximate invariance by
subset-wise affine reparame-
terization

We consider a parameter $\tilde{\theta}$ such that:

$$
\theta = \varphi(\tilde{\theta}),
$$

where $\varphi$ is an invertible map, affine on each sub-
set of parameters. Therefore, its Jacobian is: $\mathbf{J} =$
Diag($\alpha_1, \dots, \alpha_p$), where, for all $1 \le s \le S$
$1 \le p_1, p_2 \le P$, we have:

$$
p_1, p_2 \in \mathcal{I}_s \Rightarrow \alpha_{p_1} = \alpha_{p_2} =: a_s.
$$

Also, let $\bar{\mathbf{J}} = \text{Diag}(a_1, \dots, a_S)$.

We want to compare the training trajectory of $\mathcal{L}(\theta)$ and $\mathcal{L}(\varphi(\tilde{\theta}))$ when using Method 4.1. For any quantity $x$ computed with the parameterization $\theta$, we denote by $\tilde{x}$ its counterpart computed with the parameterization $\tilde{\theta}$.

We compute $\tilde{\eta}_*$. Equation (2) gives:

$$
\tilde{\eta}_* = \left( \tilde{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\tilde{\mathbf{D}}\tilde{\eta}_*\| \tilde{\mathbf{D}}^2 \right)^{-1} \tilde{\mathbf{g}}. \quad (7)
$$

Besides:

$$
\tilde{\mathbf{H}} := \mathbf{I}_{S:P} \tilde{\mathbf{U}} \tilde{\mathbf{H}} \tilde{\mathbf{U}} \mathbf{I}_{P:S}
$$

$$
\tilde{\mathbf{g}} := \mathbf{I}_{S:P} \tilde{\mathbf{U}} \tilde{\mathbf{g}}
$$

To go further, we need to do an assumption about
the direction **u**.

**Assumption E.1.** We assume that $\mathbf{U}_t$ is computed
in such a way that $\tilde{\mathbf{U}}_t = \mathbf{JU}_t$ at every step.

This assumption holds typically when $\mathbf{u}_t$ is the gra-
dient at time step $t$. It holds also when $\mathbf{u}_t$ is a linear
combination of the past gradients:

$$
\mathbf{u}_1 := \mathbf{g}_1
$$

$$
\mathbf{u}_{t+1} := \mu \mathbf{u}_t + \mu' \mathbf{g}_{t+1},
$$

which includes the momentum.

To summarize, we have:

$$
\tilde{\mathbf{U}} = \mathbf{JU}, \quad \tilde{\mathbf{H}} = \mathbf{JHJ}, \quad \tilde{\mathbf{g}} = \mathbf{Jg},
$$

So:

$$
\tilde{\mathbf{H}} = \tilde{\mathbf{J}}^2 \mathbf{I}_{S:P} \mathbf{U} \mathbf{H} \mathbf{U}_{P:S} \tilde{\mathbf{J}}^2 = \tilde{\mathbf{J}}^2 \bar{\mathbf{H}} \tilde{\mathbf{J}}^2,
$$

$$
\tilde{\mathbf{g}} = \tilde{\mathbf{J}}^2 \mathbf{I}_{S:P} \mathbf{U} \mathbf{g} = \tilde{\mathbf{J}}^2 \bar{\mathbf{g}},
$$

since **J** and **U** are diagonal. And:

$$
\text{since: } \mathbf{D}_{ii} = \left| (\mathbf{D}_{\theta}^{(3)}(\mathbf{u}))_{iii} \right|^{1/3},
$$

$$
\text{then: } \tilde{\mathbf{D}}_{ii} = a_i^2 \mathbf{D}_{ii},
$$

thus: $\tilde{\mathbf{D}} = \tilde{\mathbf{J}}^2 \mathbf{D}$.

Thus, Eqn. (7) becomes:

$$
\tilde{\eta}_* = \left( \tilde{\mathbf{J}}^2 \bar{\mathbf{H}} \tilde{\mathbf{J}}^2 + \frac{\lambda_{\text{int}}}{2} \| \tilde{\mathbf{J}}^2 \tilde{\mathbf{D}} \tilde{\eta}_* \| \tilde{\mathbf{J}}^4 \mathbf{D}^2 \right)^{-1} \tilde{\mathbf{J}}^2 \bar{\mathbf{g}},
$$

which can be rewritten (since $\tilde{\mathbf{J}}$ is invertible):

$$
\tilde{\mathbf{J}}^2 \tilde{\eta}_* = \left( \tilde{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \| \tilde{\mathbf{D}} \tilde{\mathbf{J}}^2 \tilde{\eta}_* \| \tilde{\mathbf{D}}^2 \right)^{-1} \tilde{\mathbf{g}}.
$$

Therefore, $\tilde{\eta}_*$ is a solution of Eqn. (2) in the param-
eterization $\tilde{\theta}$ if, and only if, $\tilde{\mathbf{J}}^2\tilde{\eta}_*$ is a solution in the
parameterization $\theta$. Moreover, $\|\tilde{\mathbf{D}}\tilde{\eta}_*\| = \|\mathbf{D}\mathbf{J}^2\tilde{\eta}_*\|$,
so $\tilde{\eta}_*$ is the solution of maximum norm $\|\tilde{\mathbf{D}}\tilde{\eta}_*\|$ of (2)
with parameterization $\tilde{\theta}$ iff $\tilde{\mathbf{J}}^2\tilde{\eta}_*$ is a the solution of
maximum norm $\|\mathbf{D}\mathbf{J}^2\tilde{\eta}_*\|$ of (2) with parameteriza-
tion $\theta$.

Thus, $\eta_* = \tilde{\mathbf{J}}^2 \tilde{\eta}_*$, and the update step in parame-
terization $\tilde{\theta}$ is:

$$
\begin{align*}
\tilde{\theta}_{t+1} &= \tilde{\theta}_t - \tilde{U}_t \mathbf{I}_{P:S} \tilde{\eta}_* \\
&= \tilde{\theta}_t - \tilde{U}_t \mathbf{I}_{P:S} \tilde{\mathbf{J}}^{-2} \eta_*
\end{align*}
$$

which can be rewritten:

$$
\mathbf{J}^{-1} \theta_{t+1} = \mathbf{J}^{-1} \theta_t - \mathbf{U} \mathbf{J} \mathbf{I}_{P:S} \tilde{\mathbf{J}}^{-2} \eta_*, \quad (8)
$$

since $\varphi$ is an affine function with factor **J**. Finally,
En. (8) boils down to:

$$
\theta_{t+1} = \theta_t - \mathbf{U} \mathbf{I}_{P:S} \eta_*,
$$

which is exactly Method 4.1 in parameterization $\theta$.
## F Experimental details

**Practical implementation.** To implement the method proposed in Section 4, we propose Algorithm 1. The key function is compute\_lr($\lambda_{\text{int}}; \mathcal{L}, \theta, \tilde{Z}, \mathbf{u}$), which returns a solution $\eta_{*}$ of:

$$ \eta_{*} = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta_{*}\|\mathbf{D}^{2} \right)^{-1} \bar{\mathbf{g}}, $$

$$ \text{with: } \bar{\mathbf{H}} := \mathbf{I}_{S:P}\text{Diag}(\mathbf{u})\frac{\partial^2 \mathcal{L}}{\partial\theta^2}(\theta, \tilde{Z})\text{Diag}(\mathbf{u})\mathbf{I}_{P:S}, \\ \bar{\mathbf{g}} := \mathbf{I}_{S:P}\text{Diag}(\mathbf{u})\mathbf{g}_t, $$

$$ \mathbf{D} := \text{Diag} \left( \left( \left| \mathbf{D}_{\theta}^{(3)}(\mathbf{u}) \right|_{iii}^{1/3} \right)_{1 \le i \le S} \right). $$

“momentum($\mu$, $\mathbf{x}$, $\tilde{\mathbf{x}}$)” returns $\mathbf{x}$ if $\tilde{\mathbf{x}}$ is undefined, else $\mu\tilde{\mathbf{x}} + (1 - \mu)\mathbf{x}$. “schedule($\tau_{\text{sch}}$, $p_{\text{sch}}$, $f_{\text{sch}}$; ...)” corresponds to torch.optim.lr_scheduler.ReduceLROnPlateau called every $\tau_{\text{sch}}$ with patience $p_{\text{sch}}$ and factor $f_{\text{sch}}$, in order to reduce the damping $\lambda_t$ when the loss attains a plateau.<fnref n="7" /> The samplers $\mathcal{D}_g$ and $\mathcal{D}_{\text{newt}}$ are respectively used to compute the gradients $\mathbf{g}_t$ and $(\mathbf{H}, \mathbf{g})$ used in “compute\_lr”.

The hyperparameters are: the initial damping factor $\lambda_1$, the momentum $\mu_g$ on the gradients $\mathbf{g}_t$, the minibatch size $B$ to sample the $\tilde{Z}$ (used to compute $\bar{\mathbf{g}}$, $\mathbf{H}$ and $\mathbf{D}$), the number of steps $\tau$ between each call of compute\_lr, the momentum $\mu_{\eta}$ on the learning rates $\eta_t$, the internal damping $\lambda_{\text{int}}$, and the parameters of the scheduler $\tau_{\text{sch}}$, $p_{\text{sch}}$, $f_{\text{sch}}$.

**Explanation.** The “momentum” functions are used to deal with the stochastic part of the training process, since our method has not been designed to be robust against noise. The period $\tau$ is usually strictly greater than 1, in order to avoid calling “compute\_lr” at every step, which would be costly. The minibatch size $B$ should be large enough to reduce noise in the estimation of $\eta_{*}$. If we denote by $B_g$ the size of the minibatches in $\mathcal{D}_g$, then we recommend the following setup: $\tau = \frac{B}{B_g} = \frac{1}{1-\mu_g}$. That way, we ensure that the training data are sampled from $\mathcal{D}_g$ and $\mathcal{D}_{\text{newt}}$ at

**Algorithm 1** Complete implementation of the second-order optimization method described in Section 4. $\lambda_1$ and $\lambda_{\text{int}}$ are the only hyperparameter to be tuned across the experiments, the others are fixed.

Hyperpar.: $\lambda_1, \mu_g, B_g, B, \tau, \mu_{\eta}, \lambda_{\text{int}}, \tau_{\text{sch}}, p_{\text{sch}}, f_{\text{sch}}$
$\mathcal{D}_g \leftarrow$ sampler of minibatches of size $B_g$
$\mathcal{D}_{\text{newt}} \leftarrow$ sampler of minibatches of size $B$

for all $t \in [1, T]$ do

$Z_t := (X_t, Y_t) \sim \mathcal{D}_{\text{tr}}$ (sample minibatch)

$\mathcal{L}_t \leftarrow \mathcal{L}(\theta_t, Z_t)$ (forward pass)

$\mathbf{g}_t \leftarrow \frac{\partial \mathcal{L}}{\partial \theta}(\theta_t, Z_t)$ (backward pass)

$\tilde{\mathbf{g}}_t \leftarrow$ momentum($\mu_g; \mathbf{g}_t, \tilde{\mathbf{g}}_{t-1}$)

if $t \% \tau == 0$ then

sample $\tilde{Z}_t \sim \mathcal{D}_{\text{newt}}$

$\eta_t \leftarrow$ compute\_lr($\lambda_{\text{int}}; \mathcal{L}, \theta_t, \tilde{Z}_t, \tilde{\mathbf{g}}_t$)

$\tilde{\eta}_t \leftarrow$ momentum($\mu_{\eta}; (\eta_t)_+, \tilde{\eta}_{t-1}$)

end if

$\theta_{t+1} \leftarrow \theta_t - \lambda_t \text{Diag}(\tilde{\mathbf{g}}_t) \mathbf{I}_{P:S} \tilde{\eta}_t$ (training step)

$\lambda_{t+1} \leftarrow$ schedule($\tau_{\text{sch}}, p_{\text{sch}}, f_{\text{sch}}; t, \mathcal{L}_t, \lambda_t$)

end for

the same rate, and that $\tilde{\mathbf{g}}_t$ memorizes the preceding gradients $\mathbf{g}_t$ for $\tau$ steps. Besides, we have to take the positive part $(\eta_t)_+$ of $\eta_t$ in order to avoid negative learning rates.

**Experimental setup.** We provide in Table 2 the hyperparameters fixed for all the experiments. In Table 3, we report the results of the grid-search for the hyperparameters of the 3 tested optimization methods.

<table><caption>Table 2: Hyperparameters fixed in all the series of experiments. $N_e$ is the number of training steps per epoch.</caption><thead><tr><th>μ<sub>g</sub></th><th>B<sub>g</sub></th><th>B</th><th>τ</th><th>μ<sub>η</sub></th><th>τ<sub>sch</sub></th><th>p<sub>sch</sub></th><th>f<sub>sch</sub></th></tr></thead><tbody><tr><td>0.9</td><td>10<sup>2</sup></td><td>10<sup>3</sup></td><td>10</td><td>0.5</td><td>N<sub>e</sub></td><td>5</td><td>0.5</td></tr></tbody></table>

<footnote n="7">⁷See `torch.optim.lr_scheduler.ReduceLROnPlateau`.</footnote>
## F Experimental details

**Practical implementation.** To implement the method proposed in Section 4, we propose Algorithm 1. The key function is `compute_lr`($\lambda_{\text{int}}$; $\mathcal{L}, \theta, \tilde{Z}, \mathbf{u}$), which returns a solution $\eta_{*}$ of:

$$ \eta_{*} = \left( \bar{\mathbf{H}} + \frac{\lambda_{\text{int}}}{2} \|\mathbf{D}\eta_{*}\|\mathbf{D}^{2} \right)^{-1} \bar{\mathbf{g}}, $$

$$ \text{with: } \bar{\mathbf{H}} := \mathbf{I}_{S:P}\text{Diag}(\mathbf{u})\frac{\partial^2 \mathcal{L}}{\partial\theta^2}(\theta, \tilde{Z})\text{Diag}(\mathbf{u})\mathbf{I}_{P:S}, $$

$$ \bar{\mathbf{g}} := \mathbf{I}_{S:P}\text{Diag}(\mathbf{u})\mathbf{g}_t, $$

$$ \mathbf{D} := \text{Diag} \left( \left( \left| \mathbf{D}_{\theta}^{(3)}(\mathbf{u}) \right|_{iii}^{1/3} \right)_{1 \le i \le S} \right). $$

“momentum($\mu$, $\mathbf{x}$, $\tilde{\mathbf{x}}$)” returns $\mathbf{x}$ if $\tilde{\mathbf{x}}$ is undefined, else $\mu\tilde{\mathbf{x}} + (1 - \mu)\mathbf{x}$. “schedule($\tau_{\text{sch}}$, $p_{\text{sch}}$, $f_{\text{sch}}$; ...)” corresponds to `torch.optim.lr_scheduler ReduceLROnPlateau` called every $\tau_{\text{sch}}$ with patience $p_{\text{sch}}$ and factor $f_{\text{sch}}$, in order to reduce the damping $\lambda_t$ when the loss attains a plateau.<fnref n="7" /> The samplers $\mathcal{D}_g$ and $\mathcal{D}_{\text{newt}}$ are respectively used to compute the gradients $\mathbf{g}_t$ and $(\mathbf{H}, \mathbf{g})$ used in “compute_lr”.

The hyperparameters are: the initial damping factor $\lambda_1$, the momentum $\mu_g$ on the gradients $\mathbf{g}_t$, the minibatch size $B$ to sample the $\tilde{Z}$ (used to compute $\bar{\mathbf{g}}$, $\mathbf{H}$ and $\mathbf{D}$), the number of steps $\tau$ between each call of `compute_lr`, the momentum $\mu_\eta$ on the learning rates $\eta_t$, the internal damping $\lambda_{\text{int}}$, and the parameters of the scheduler $\tau_{\text{sch}}$, $p_{\text{sch}}$, $f_{\text{sch}}$.

**Explanation.** The “momentum” functions are used to deal with the stochastic part of the training process, since our method has not been designed to be robust against noise. The period $\tau$ is usually strictly greater than 1, in order to avoid calling “compute_lr” at every step, which would be costly. The minibatch size $B$ should be large enough to reduce noise in the estimation of $\eta_{*}$. If we denote by $B_g$ the size of the minibatches in $\mathcal{D}_g$, then we recommend the following setup: $\tau = \frac{B}{B_g} = \frac{1}{1-\mu_g}$. That way, we ensure that the training data are sampled from $\mathcal{D}_g$ and $\mathcal{D}_{\text{newt}}$ at

**Algorithm 1** Complete implementation of the second-order optimization method described in Section 4. $\lambda_1$ and $\lambda_{\text{int}}$ are the only hyperparameter to be tuned across the experiments, the others are fixed.

Hyperpar.: $\lambda_1, \mu_g, B_g, B, \tau, \mu_\eta, \lambda_{\text{int}}, \tau_{\text{sch}}, p_{\text{sch}}, f_{\text{sch}}$
$\mathcal{D}_g \leftarrow$ sampler of minibatches of size $B_g$
$\mathcal{D}_{\text{newt}} \leftarrow$ sampler of minibatches of size $B$
**for all** $t \in [1, T]$ **do**
$Z_t := (X_t, Y_t) \sim \mathcal{D}_{\text{tr}}$ (sample minibatch)
$\mathcal{L}_t \leftarrow \mathcal{L}(\theta_t, Z_t)$ (forward pass)
$\mathbf{g}_t \leftarrow \frac{\partial \mathcal{L}}{\partial \theta}(\theta_t, Z_t)$ (backward pass)
$\tilde{\mathbf{g}}_t \leftarrow$ momentum($\mu_g; \mathbf{g}_t, \tilde{\mathbf{g}}_{t-1}$)
**if** $t \% \tau == 0$ **then**
sample $\tilde{Z}_t \sim \mathcal{D}_{\text{newt}}$
$\eta_t \leftarrow$ compute_lr($\lambda_{\text{int}}; \mathcal{L}, \theta_t, \tilde{Z}_t, \tilde{\mathbf{g}}_t$)
$\tilde{\eta}_t \leftarrow$ momentum($\mu_\eta; (\eta_t)_+, \tilde{\eta}_{t-1}$)
**end if**
$\theta_{t+1} \leftarrow \theta_t - \lambda_t \text{Diag}(\tilde{\mathbf{g}}_t) \mathbf{I}_{P:S} \tilde{\eta}_t$ (training step)
$\lambda_{t+1} \leftarrow$ schedule($\tau_{\text{sch}}, p_{\text{sch}}, f_{\text{sch}}; t, \mathcal{L}_t, \lambda_t$)
**end for**

the same rate, and that $\tilde{\mathbf{g}}_t$ memorizes the preceding gradients $\mathbf{g}_t$ for $\tau$ steps. Besides, we have to take the positive part $(\eta_t)_+$ of $\eta_t$ in order to avoid negative learning rates.

**Experimental setup.** We provide in Table 2 the hyperparameters fixed for all the experiments. In Table 3, we report the results of the grid-search for the hyperparameters of the 3 tested optimization methods.

<table><caption>Table 2: Hyperparameters fixed in all the series of experiments. $N_e$ is the number of training steps per epoch.</caption><thead><tr><th>μ<sub>g</sub></th><th>B<sub>g</sub></th><th>B</th><th>τ</th><th>μ<sub>η</sub></th><th>τ<sub>sch</sub></th><th>p<sub>sch</sub></th><th>f<sub>sch</sub></th></tr></thead><tbody><tr><td>0.9</td><td>10<sup>2</sup></td><td>10<sup>3</sup></td><td>10</td><td>0.5</td><td>N<sub>e</sub></td><td>5</td><td>0.5</td></tr></tbody></table>

<footnote n="7">⁷See `torch.optim.lr_scheduler ReduceLROnPlateau`.</footnote>
<table><caption>Table 3: Hyperparams tuned for each series of experiments. $\eta$: learning rate, $\lambda_1$: initial damping factor.</caption><thead><tr><th></th><th>MLP</th><th>LeNet</th><th>VGG-11'</th><th>BigMLP</th></tr></thead><tbody><tr><td>Adam: $\eta$</td><td>3 &middot; 10<sup>-4</sup></td><td>3 &middot; 10<sup>-4</sup></td><td>10<sup>-5</sup></td><td>10<sup>-5</sup></td></tr><tr><td>KFAC: $\eta$</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>3 &middot; 10<sup>-4</sup></td><td>10<sup>-5</sup></td></tr><tr><td>KFAC: $\lambda_1$</td><td>10<sup>-2</sup></td><td>3 &middot; 10<sup>-2</sup></td><td>3 &middot; 10<sup>-2</sup></td><td>10<sup>-2</sup></td></tr><tr><td>Ours: $\lambda_1$</td><td>10<sup>-1</sup></td><td>3 &middot; 10<sup>-1</sup></td><td>3 &middot; 10<sup>-1</sup></td><td>10<sup>-1</sup></td></tr><tr><td>Ours: $\lambda_{int}$</td><td>10</td><td>3</td><td>3</td><td>10</td></tr></tbody></table>

## G Very deep multilayer perceptron

**Grouping the layers.** In addition to the neural networks considered in Section 5, we have also tested “VBigMLP”, a very deep multilayer perceptron with 100 layers of size 1024 trained on CIFAR-10. Instead of considering $S = 2L = 200$ groups of parameters, we split the sequence of layers of VBigMLP into 5 chunks. Then, each chunk is divided into 2 parts, one containing the weight tensors, and the other the bias vectors. Finally, we have $S = 10$ subsets of parameters, grouped by role (weight/bias) and by position inside the network.

**Experimental results.** We show in Figure 5 the matrices $\mathbf{H}$ and $\mathbf{H}^{-1}$ at different stages of training. At initialization, even if the neural network is very deep, we observe that all the chunks of the network interact together, even the first one with the last one. However, after several training steps, the long-range interactions seem to disappear. Incidentally, the matrices become tridiagonal, which ties in with the block-tridiagonal approximation of the inverse of the Hessian done by Martens and Grosse (2015).

In Figure 6, we observe the evolution of the learning rates $\eta_*$ computed according to (2). First, there are all increasing during training. Second, the biases in the last layers of the network seem to need larger learning rates than biases in the first layers. Third, the learning rate computed for the weights of the first chunk of layers is smaller than the others.

Finally, the training curves in Figure 7 indicate

<figure><img src="image_5.png" /><figcaption>Figure 5: VBigMLP trained by SGD on CIFAR-10. Submatrices of $\mathbf{H}$ (first row) and $\mathbf{H}^{-1}$ (second row), at initialization and before the 100th epoch.</figcaption></figure>

that our method can be used to train very deep networks, but quite slowly. In this setup, it is far from being competitive with Adam. Besides, we did not manage to tune the learning rate and the damping of K-FAC to make it work in this setup.
<table><caption>Table 3: Hyperparams tuned for each series of experiments. $\eta$: learning rate, $\lambda_1$: initial damping factor.</caption><thead><tr><th></th><th>MLP</th><th>LeNet</th><th>VGG-11'</th><th>BigMLP</th></tr></thead><tbody><tr><td>Adam: $\eta$</td><td>3 &middot; 10<sup>-4</sup></td><td>3 &middot; 10<sup>-4</sup></td><td>10<sup>-5</sup></td><td>10<sup>-5</sup></td></tr><tr><td>KFAC: $\eta$</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>3 &middot; 10<sup>-4</sup></td><td>10<sup>-5</sup></td></tr><tr><td>KFAC: $\lambda_1$</td><td>10<sup>-2</sup></td><td>3 &middot; 10<sup>-2</sup></td><td>3 &middot; 10<sup>-2</sup></td><td>10<sup>-2</sup></td></tr><tr><td>Ours: $\lambda_1$</td><td>10<sup>-1</sup></td><td>3 &middot; 10<sup>-1</sup></td><td>3 &middot; 10<sup>-1</sup></td><td>10<sup>-1</sup></td></tr><tr><td>Ours: $\lambda_{int}$</td><td>10</td><td>3</td><td>3</td><td>10</td></tr></tbody></table>

## G Very deep multilayer perceptron

**Grouping the layers.** In addition to the neural networks considered in Section 5, we have also tested “VBigMLP”, a very deep multilayer perceptron with 100 layers of size 1024 trained on CIFAR-10. Instead of considering $S = 2L = 200$ groups of parameters, we split the sequence of layers of VBigMLP into 5 chunks. Then, each chunk is divided into 2 parts, one containing the weight tensors, and the other the bias vectors. Finally, we have $S = 10$ subsets of parameters, grouped by role (weight/bias) and by position inside the network.

**Experimental results.** We show in Figure 5 the matrices $\mathbf{H}$ and $\mathbf{H}^{-1}$ at different stages of training. At initialization, even if the neural network is very deep, we observe that all the chunks of the network interact together, even the first one with the last one. However, after several training steps, the long-range interactions seem to disappear. Incidentally, the matrices become tridiagonal, which ties in with the block-tridiagonal approximation of the inverse of the Hessian done by Martens and Grosse (2015).

In Figure 6, we observe the evolution of the learning rates $\eta_*$ computed according to (2). First, there are all increasing during training. Second, the biases in the last layers of the network seem to need larger learning rates than biases in the first layers. Third, the learning rate computed for the weights of the first chunk of layers is smaller than the others.

Finally, the training curves in Figure 7 indicate

<figure><img src="image_7.png" /><figcaption>Figure 5: VBigMLP trained by SGD on CIFAR-10. Submatrices of $\mathbf{H}$ (first row) and $\mathbf{H}^{-1}$ (second row), at initialization and before the 100th epoch.</figcaption></figure>

that our method can be used to train very deep networks, but quite slowly. In this setup, it is far from being competitive with Adam. Besides, we did not manage to tune the learning rate and the damping of K-FAC to make it work in this setup.
<figure><img src="image_6.png" /><figcaption>Figure 6: VBigMLP trained by SGD on CIFAR-10. Legend: solid lines: weights; dotted lines: biases. Learning rates $\eta_*$ computed according to (2), specific to each subset of parameters. For each epoch $k \in \{25, 75, 125, 175, 225, 275, 325, 375\}$, the reported value has been averaged over the epochs $[k-25, k+24]$ to remove the noise.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Training curves with VBigMLP, CIFAR-10. The training of Adam has been stopped at 200 epochs, since no significant improvement was observed in the last 50 epochs.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: VBigMLP trained by SGD on CIFAR-10. Legend: solid lines: weights; dotted lines: biases. Learning rates $\eta_*$ computed according to (2), specific to each subset of parameters. For each epoch $k \in \{25, 75, 125, 175, 225, 275, 325, 375\}$, the reported value has been averaged over the epochs $[k-25, k+24]$ to remove the noise.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Training curves with VBigMLP, CIFAR-10. The training of Adam has been stopped at 200 epochs, since no significant improvement was observed in the last 50 epochs.</figcaption></figure>