††This work was performed under the auspices of the U.S. Department of Energy
by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344
(LLNL-JRNL-827231) and was supported by the LLNL-LDRD program under Project
No. 19-ERD-019.

# Multilevel-in-Layer Training for Deep Neural Network Regression

Colin Ponce
Center for Applied Scientific Computing
Lawrence Livermore National Laboratory
<EMAIL_ADDRESS>&Ruipeng Li
Center for Applied Scientific Computing
Lawrence Livermore National Laboratory
<EMAIL_ADDRESS>&Christina Mao
Global Security Computing Applications Division
Lawrence Livermore National Laboratory
<EMAIL_ADDRESS>&Panayot Vassilevski
Fariborz Maseeh Department of Mathematics + Statistics
Portland State University
<EMAIL_ADDRESS>
###### Abstract

A common challenge in regression is that for many problems, the degrees of
freedom required for a high-quality solution also allows for overfitting.
Regularization is a class of strategies that seek to restrict the range of
possible solutions so as to discourage overfitting while still enabling good
solutions, and different regularization strategies impose different types of
restrictions. In this paper, we present a multilevel regularization strategy
that constructs and trains a hierarchy of neural networks, each of which has
layers that are wider versions of the previous network’s layers. We draw
intuition and techniques from the field of Algebraic Multigrid (AMG),
traditionally used for solving linear and nonlinear systems of equations, and
specifically adapt the Full Approximation Scheme (FAS) for nonlinear systems
of equations to the problem of deep learning. Training through V-cycles then
encourage the neural networks to build a hierarchical understanding of the
problem. We refer to this approach as _multilevel-in-width_ to distinguish
from prior multilevel works which hierarchically alter the depth of neural
networks. The resulting approach is a highly flexible framework that can be
applied to a variety of layer types, which we demonstrate with both fully-
connected and convolutional layers. We experimentally show with PDE regression
problems that our multilevel training approach is an effective regularizer,
improving the generalize performance of the neural networks studied.

_Keywords_ Full Approximation Scheme $\cdot$ FAS $\cdot$ Multigrid $\cdot$
Multilevel Methods $\cdot$ DNN Training $\cdot$ Neural Networks

##  1 Introduction

A feedforward deep neural network (DNN) consists of a set of layers, each of
which applies an affine mapping to its input vectors followed by a nonlinear
_activation function_ acting componentwise on the transferred vectors. That
affine mapping is controlled by a set of tunable parameters, the selection of
which determines the overall action of the DNN. The task of training a deep
neural network (DNN) is that of selecting these network parameters in order to
enable the DNN to approximately capture a desired pattern. This is done by
solving, or approximately solving, an optimization problem with an objective
function that maps a choice of DNN parameters to a quantitative measure of the
DNN’s performance on the problem of interest; this objective function is
commonly called a _loss function_ [11].

Training is typically accomplished through stochastic gradient descent (SGD)
or its variants, which are cost-efficient and well-supported by existing
software [1, 20]. A common challenge in the design of neural networks and
their training is that, often, providing sufficient degrees of freedom in the
network to obtain a high-quality solution also allows for overfitting, which
harms generalization accuracy. _Regularization_ is a class of strategies that
seeks to restrict the freedom of the training method, therefore encouraging
simpler solutions. Ideally, this regularization discourages overfitting while
still enabling high-quality solutions. A variety of regularization strategies
exist, each of which encourages a different kind of simplicity, such as L2
regularization, L1 regularization, or dropout regulariztaion [11].

In this paper, we present a novel _multilevel_ regularization strategy that
constructs and trains a hierarchy of neural networks. In this approach, all
networks have the same number and type of layers, but each layer is a wider
version of next network’s layers, so networks get ”thinner” as we move through
the hierarchy. We base this approach on the techniques of Algebraic Multigrid
(AMG), a class of methods traditionally used for solving linear and nonlinear
systems of equations, and we leverage teh Full Approximation Scheme (FAS), a
multilevel framework traditionally used for nonlinear systems of equations. We
use the classic V-cycle in training, and develop a method to pass problem and
parameter information, including momentum information, between the levels of
the hierarchy during training. The resulting training method encourages the
neural networks to learn the problem in a hierarchical manner. This approach
is highly flexible, allowing each layer to coarsen (or not) independently of
each other layer, and can be applied to a variety of layer types. We
specifically demonstrate this for fully-connected and convolutional layers.
Our approach is not a replacement for traditional, one-level learning methods
such as SGD; instead, it uses that traditional optimizer as an inner iterator
within the multilevel framework. We demonstrate the generalization benefits of
our approach with a set of PDE-based regression problems.

We refer to our approach as ”multilevel-in-width” to distinguish it from prior
multilevel works that change the number of layers in a neural network, which
we think of as ”multilevel-in-depth.” Other authors, taking very different
approaches from that presented here, have also observed regularization
benefits from their multilevel methods [27, 7]. The present paper contributes
to a small but growing body of evidence that suggests that multilevel neural
network training methods have the potential to act as a beneficial regularizer
over a significant range of problems, network architectures, and multilevel
approaches, though this requires further investigation to fully understand.

Our open-source software, which we call MTNN (short for _Multilevel Training
of Neural Networks_ , pronounced “mountain”) is built on PyTorch and is free
to download at https://github.com/LLNL/MTNN. It is designed for modularity and
extensibility to facilitate the study of new problems and algorithmic choices.

In Section 2, we discuss prior work in AMG, current methods for training
neural networks, and prior work on building notions of hierarchy into DNNs. We
then give a summary of the FAS algorithm from an optimization point of view in
Section 3. Next, in Section 4, we formulate a two-level FAS in the setting of
DNN training. Our specific choices of the operators (restriction and
prolongation) are given in Section 5. Section 6 presents our two-level full-
approximation scheme, which can be recursively extended to more levels.
Section 7 presents the extension of our method to convolutional neural
networks. Finally, in Section 8 we demonstrate our techniques on a set of PDE
regression problems, and in Section 9 we conclude.

##  2 Background and Related Work

###  2.1 Multilevel Methods

Multilevel methods are a class of iterative methods long used in scientific
computing as a scalable solver for linear and non-linear systems of equations
$g(\mathbf{x})=\mathbf{f}$. Multilevel methods is a more general class of
algorithms emerging out of _multigrid methods_ , where this method was first
applied to partial differential problems on grid. Traditional multigrid
algorithms construct a hierarchy of discretizations of a geometric problem
domain in which each discretization is a coarsened variant of the previous.
Algebraic multigrid (AMG) also constructs a hierarchy of problems, but does in
a fully-algebraic manner, independent of any underlying geometry.

When applied to (deterministic) systems of equations, multilevel methods
accelerate convergence towards a solution through _smoothing_ , also known as
relaxation methods, and _coarse-grid corrections_ , down-sampling residual
error on the coarse grid back to the fine grid. Passing solution information
between fine and coarse grids is known as _restriction_ and _prolongation_.

A single multilevel iteration is known as a _cycle_. Perhaps the most common
of these is the V-cycle, the linear two-level version of which proceeds as
follows:

  1. 1.

Smooth the solution iterate on the fine level problem $g$.

  2. 2.

Restrict problem residual to the coarse level problem $g_{c}$ with some
restriction operator $R$.

  3. 3.

Solve the problem on the coarse level to obtain a coarse-grid correction.

  4. 4.

Prolong the coarse-grid correction to the fine level and add it to the
solution iterate with some prolongation operator $P$.

  5. 5.

Smooth the solution iterate on the fine level problem $g$.

This V-cycle is a framework in which one must make algorithmic choices on how
to accomplish individual steps such as smoothing, restriction, and
prolongation. For further reference, a more detailed introduction to
multilevel and multigrid methods can be provided by Briggs et al. [5],
Hackbusch [13] and Wesseling [28].

A classic problem for which linear multilevel methods are extremely effective
is the Poisson problem, $-\nabla\leavevmode\nobreak\ \cdot\leavevmode\nobreak\
(\kappa\leavevmode\nobreak\ \nabla\leavevmode\nobreak\ u)\leavevmode\nobreak\
=\leavevmode\nobreak\ f$, defined over some problem domain, which models
diffusion processes. As a result, data from Poisson problems are often
amenable to the type of hierarchical representations discussed in the
Introduction.

####  2.1.1 Full Approximation Scheme

When multilevel methods are applied to nonlinear systems of equations, the
V-cycle is similar, but passing the residual to the coarse level may no longer
be appropriate. A generalization appropriate for nonlinear problems is the
_Full Approximation Scheme_ (FAS) [4]. In FAS, one still has a fine problem,
$g$, a coarse problem, $g_{c}$, and restriction and prolongation operators,
$R$ and $P$.

However, instead of restricting the residual as in Step 2 above, one instead
restricts the solution iterate, and alters the right-hand side of the coarse
problem to capture residual information. The FAS coarse-grid problem, then, is

${\widehat{g}_{c}}(\mathbf{x}_{c})=R\mathbf{b}+(g_{c}(R\mathbf{x})-Rg(\mathbf{x})).$
(1)

The above modification to the right-hand side is sometimes called a _$\tau$
-correction_; one can check that in the linear case it reduces back to the
coarse-grid residual equation.

An extension of FAS to discretized optimization problems, “MG/OPT”, was
introduced by Nash, who also demonstrated its convergence properties [21]. The
MG/OPT algorithm extends FAS by performing a line search on applying the
coarse-grid correction (Step 4 in the above V-cycle).

###  2.2 Deep Neural Networks (DNN)

A neural network (NN) consists of $n_{L}$ layers: an input layer, $n_{L}-2$
hidden layers, and an output layer. An input vector is denoted by
$\mathbf{y}^{in}$. The hidden layers consists of weight matrices $W_{k}$ and
bias vectors $\mathbf{b}_{k}$, $k=1,...,n_{L}$, and an activation function
$\varphi$ such as the well-known ReLU function. The output of the neural
network calculation is denoted by $\mathbf{y}^{out}$. There is no clear
definition for the number of hidden layers required to consider a neural
network a “deep” neural network, though often any network with more than 2
hidden layers is considered a DNN.

Forward propagation, also known as a _forward pass_ , is the procedure whereby
the input (feature) vector is passed through the neural network layer by
layer, applying an activation function after each layer to produce
intermediate variables and finally an output vector. We discuss the notation
and operation of a forward pass in greater detail in Section 4.

####  2.2.1 Training Methods for Neural Networks

An important difference between neural network training and traditional
optimization is that the true objective function of interest, performance over
the entire population of possible examples, is inaccessible; the performance
over a training sample is merely an approximation to this true function of
interest. This motivates a number of important design choices in neural
network training methods. Training is typically done via a stochastic first-
order method, such as stochastic gradient descent, AdaGrad [8], RMSProp, Adam
[15], or other variations [19, 18]; this both accelerates computation and
reduces the chances of getting stuck in a local minimum that is an artifact is
the particular training set. Momentum maintains an exponentially decaying sum
of past gradients, which helps to accelerate training and jump over local
minima by allowing the ”big ideas” to compound [25]. Finally, training all the
way to convergence on the training set would likely result in severe
overfitting, even with regularization, and so early stopping is almost
universal [11].

These methods are all fundamentally variations on gradient descent. Second-
order methods such as Newton-Raphson have been standard in optimization for
some time [3], and in fact is often much more effective than gradient descent-
style methods. However, second order methods require the computation or
approximation of the Hessian matrix of the objective function, requiring the
computation of $O(n^{2})$ Hessian elements for $n$ parameters. In the case of
DNNs, even moderately sized networks contain millions of parameters, and so a
Hessian matrix is very slow to compute (and sometimes does not fit in memory);
consequently, second-order methods have not become mainstream for DNNs.
Approximate second-order methods for DNNs with linear storage requirements,
however, continues to be an active area of research and are beginning to
become more competitive [26].

####  2.2.2 Hierarchy and Multilevel Methods in Neural Network Architectures

Although SGD and its variations do not naturally capture hierarchy, notions of
hierarchy are quite common in neural network architectures. Convolution neural
networks (CNNs) for image processing are often described as building up
concepts in a hierarchical manner, capturing simple ideas such as edges and
shading in the first layer, and capturing more complex image elements in later
layers. Beyond CNNs, the autoencoder is one of the most popular architectural
designs today [11]. There are many variations and applications of autoencoders
but they all fundamentally work by hierarchically compressing input data into
a low-dimension representation before hierarchically decompressing it back to
its original size.

Another network architecture that deserves special mention is U-net [24]. This
network, through its layers, progressively makes an images smaller and smaller
(less height and width) but deeper and deeper (more “color” channels), before
then passing the data through layers that make it smaller and shallower again,
eventually producing a processed image of the original size. The architecture
design is highly analogous to the V-Cycle that is ubiquitous in multilevel
methods [5], and works because it processes images in a hierarchical fashion.

In recent years, a small but growing body of literature has explored
multilevel methods in the training itself. In [6], the authors develop a
multilevel-in-width Levenberg-Marquardt method for training a neural network
with one hidden layer to find a discretization-free approximate solution of a
PDE. Most other works currently in this space focus on multilevel-in-depth
methods for deep residual networks (ResNets), exploiting the dynamical
perspective that views a ResNet as taking forward Euler steps on an ODE:
Günther et. al. [12] uses multilevel parallel-in-time methods to achieve
parallelism in the forward and backward propagation; Kirby et. al. [16]
extends this approach to GPUs, and Cyr et. al. [7] uses this framework to
develop a cascadic initialization strategy with regularization benefits.
Gaedke-Merzhäuser et. al., [10] and Planta et. al. [27] also study ResNets,
but develop variants of MG/OPT to solve the problem in contrast to the
parallel-in-time approach.

Both [27] and [7] show that, even though multilevel methods are often
motivated by speed, the networks that result from multilevel training methods
can have better performance than those trained without them. We observe this
in our present work as well.

##  3 Application of Full Approximation Scheme (FAS) to an optimization
problem

In this section, we develop the framework for applying the Full Approximation
Scheme (FAS) to optimization problems and present some results that will guide
the development of FAS-based methods for deep learning.

Our optimization problem of interest is

$g(\mathbf{x})\mapsto\min,$ (2)

for a given nonlinear functional $g:\;{\mathbb{R}}^{n}\mapsto{\mathbb{R}}$. We
assume that the gradient of $g$, $\nabla g$ is computationally available, and
we will be using (in the presentation) the Hessian of $g$, denoted $H_{g}$,
which is a matrix of its second partial derivatives. Using a Taylor expansion,
we have the following approximate expansion

$g(\mathbf{x}+\mathbf{q})\approx g(\mathbf{x})+\mathbf{q}^{T}\nabla
g(\mathbf{x})+\frac{1}{2}\mathbf{q}^{T}H_{g}(\mathbf{x})\mathbf{q}.$

Now, to define a two-level FAS, we need a coarse version of $g$,
$g_{c}:\;{\mathbb{R}}^{n_{c}}\mapsto{\mathbb{R}}$, where $n_{c}\ll n$. We
introduce three mappings:

  * •

A coarse approximation to a fine quantity,
$\mathbf{x}_{c}={\boldsymbol{\Pi}}\mathbf{x}$ for a mapping
${\boldsymbol{\Pi}}:\;{\mathbb{R}}^{n}\mapsto{\mathbb{R}}^{n_{c}}$.

  * •

A fine-to-coarse restriction mapping
$\mathbf{R}:\;{\mathbb{R}}^{n}\mapsto{\mathbb{R}}^{n_{c}}$.

  * •

A coarse-to-fine interpolation mapping
${\mathbf{P}}:\;{\mathbb{R}}^{n_{c}}\mapsto{\mathbb{R}}^{n}$.

To be precise, given a fine vector $\mathbf{x}$, the coarse vector
$\mathbf{x}_{c}={\boldsymbol{\Pi}}\mathbf{x}$ is such that
${\mathbf{P}}\mathbf{x}_{c}={\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}$
approximates (in some norm) the vector $\mathbf{x}$.

The coarse FAS functional ${\widehat{g}}_{c}$ is defined in terms of the
coarse functional $g_{c}$, and in TL FAS is chosen so that ${\widehat{g}}_{c}$
approximates $g$. Specifically, let $\mathbf{x}$ be a current approximation to
the solution of the original problem. Then, the FAS method is based on the
solution of the auxiliary coarse problem

${\widehat{g}_{c}}(\mathbf{x}_{c})\equiv
g_{c}(\mathbf{x}_{c})-\mathbf{x}^{T}_{c}\left(\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})-\mathbf{R}\nabla
g(\mathbf{x})\right)\mapsto\min.$ (3)

Once $\mathbf{x}_{c}$ is computed, the next approximation to (2) is

$\mathbf{x}:=\mathbf{x}+{\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}).$
(4)

The term $\nabla g_{c}({\boldsymbol{\Pi}}\mathbf{x})-\mathbf{R}\nabla
g(\mathbf{x})$ is sometimes refereed to as the _$\tau$ -correction_. It
replaces $\nabla g_{c}({\boldsymbol{\Pi}}\mathbf{x})$, the gradient of the
coarse functional at the coarse approxmation point, with $\mathbf{R}\nabla
g(\mathbf{x})$, the restriction of the fine gradient. When computed exactly,
this makes ${\widehat{g}_{c}}(\mathbf{x}_{c})$ match, up to first order, a
restricted version of $g(\mathbf{x})$ (though as we will see in Section 5.6 we
do not compute this exactly).

###  3.1 The FAS coarse problem solution

Let $\mathbf{x}_{c}$ solve the coarse problem (3), so that
$\nabla{\widehat{g}}(\mathbf{x}_{c})=0$. Hence,

$\nabla g_{c}(\mathbf{x}_{c})-\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})+\mathbf{R}\nabla g(\mathbf{x})=0.$

Using Taylor expansion of $\nabla g_{c}$ about ${\boldsymbol{\Pi}}\mathbf{x}$,
gives

$\nabla g_{c}(\mathbf{x}_{c})=\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})+\sum\limits_{|{\boldsymbol{\alpha}}|\geq
1}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})^{{\boldsymbol{\alpha}}}\partial^{{\boldsymbol{\alpha}}}\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x}).$

Above (and in what follows), we use multi-index notation, whereby for any
${\boldsymbol{\alpha}}=(\alpha_{i})^{d}_{i=1}\in{\mathbb{R}}^{d}$ and any
vector $\mathbf{y}=(y_{i})^{d}_{i=1}\in{\mathbb{R}}^{d}$,
$\mathbf{y}^{{\boldsymbol{\alpha}}}$ is shorthand for $\prod\limits_{i:\;1\leq
i\leq d}y^{\alpha_{i}}_{i}$. Also, for a scalar function $f=f(\mathbf{x})$ of
several variables $\mathbf{x}=(x_{i})^{d}_{i=1}\in{\mathbb{R}}^{d}$, we denote
$\partial_{i}f=\frac{\partial f}{\partial x_{i}}$, and by
$\partial^{{\boldsymbol{\alpha}}}f(\mathbf{x})$ the repeated partial
derivatives of corresponding order $\alpha_{i}$ with respect to the $i$th
variable $x_{i}$ of $f$.

Combining the last two equations above, we obtain

$\mathbf{R}\nabla g(\mathbf{x})+\sum\limits_{|{\boldsymbol{\alpha}}|\geq
1}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})^{{\boldsymbol{\alpha}}}\partial^{{\boldsymbol{\alpha}}}\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})=0.$ (5)

###  3.2 The FAS-correction

Once we have computed $\mathbf{x}_{c}$, we update the new fine approximation
$\mathbf{x}_{new}=\mathbf{x}+{\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})$.
Our goal is to see how close the updated $\mathbf{x}_{new}$ is to the exact
solution, or rather how close $\nabla g(\mathbf{x}_{new})$ is to zero. For
that purpose, we expand $\nabla g(\mathbf{x}_{new})$ about $\mathbf{x}$ and
multiply by $\mathbf{R}$ on the left which gives

$\mathbf{R}\nabla g(\mathbf{x}_{new})=\mathbf{R}\nabla
g(\mathbf{x})+\sum\limits_{|{\boldsymbol{\alpha}}|\geq
1}({\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}))^{{\boldsymbol{\alpha}}}\mathbf{R}\partial^{{\boldsymbol{\alpha}}}\nabla
g(\mathbf{x}).$ (6)

If we want to make the first order terms of (6) and (5) equal, we obtain

$\sum\limits_{|{\boldsymbol{\alpha}}|=1}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})^{{\boldsymbol{\alpha}}}\partial^{{\boldsymbol{\alpha}}}\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})=\sum\limits_{|{\boldsymbol{\alpha}}|=1}({\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}))^{{\boldsymbol{\alpha}}}\mathbf{R}\partial^{{\boldsymbol{\alpha}}}\nabla
g(\mathbf{x}),$

which one can rewrite as

$\sum\limits^{n_{c}}_{i=1}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})_{i}\partial_{i}\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})=\sum\limits^{n}_{i=1}({\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}))_{i}\mathbf{R}\partial_{i}\nabla
g(\mathbf{x}).$ (7)

Introducing the Hessian matrices (note that they are symmetric)

$H_{g}(\mathbf{y})=(\partial_{i}\partial_{j}g(\mathbf{y}))^{n}_{i,j=1}\text{
and
}H_{g_{c}}(\mathbf{y}_{c})=(\partial_{i}\partial_{j}g_{c}(\mathbf{y}_{c}))^{n_{c}}_{i,j=1},$

identity (7) can be written more concisely as

$H_{g_{c}}({\boldsymbol{\Pi}}\mathbf{x})(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})=\mathbf{R}H_{g}(\mathbf{x}){\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}).$

In other words, we have proved the following result.

######  Theorem 3.1.

Let $\mathbf{R},{\mathbf{P}},{\boldsymbol{\Pi}}$ and $g$, $g_{c}$ satisfy the
identity

$\mathbf{R}H_{g}(\mathbf{x}){\mathbf{P}}=H_{g_{c}}({\boldsymbol{\Pi}}\mathbf{x}).$
(8)

Then the updated approximation $\mathbf{x}_{new}$ by the FAS algorithm
satisfies the approximate first order optimality condition

$\mathbf{R}\nabla g(\mathbf{x}_{new})=0.$

∎

To satisfy condition (8), since $H_{g_{c}}$ is a symmetric matrix, we require

$\mathbf{R}={\mathbf{P}}^{T}.$ (9)

For the application of our main interest (training of DNNs), we do not form
the mappings ${\mathbf{P}},\mathbf{R},{\boldsymbol{\Pi}}$ as matrices, but
will instead derive methods for computing the action of these operators. Also,
in general we do not choose $\mathbf{R}={\boldsymbol{\Pi}}$. Instead, we
choose ${\boldsymbol{\Pi}}$ and ${\mathbf{P}}$ such that the product
$\mathbf{Q}={\mathbf{P}}{\boldsymbol{\Pi}}$ admits an approximation property
(that we motivate in Section 3.3) and such that

${\boldsymbol{\Pi}}{\mathbf{P}}=I,$ (10)

which implies that $\mathbf{Q}$ is a projection ($\mathbf{Q}^{2}=\mathbf{Q}$).

######  Remark 3.1.

The above theorem shows that the gradient of the functional $g$ is orthogonal
(up to first order terms) to the coarse space
$\text{Range}(\mathbf{R}^{T})=\text{Range}({\mathbf{P}})$. Hence, what remains
to be done, is to reduce the error in a space complementary to the coarse
space, which is customarily referred to as _smoothing_. In other words, we are
looking at $\mathbf{x}^{\perp}$ such that

$\nabla g(\mathbf{x}_{new}+\mathbf{x}^{\perp})=0.$

The two-level FAS algorithm alternates between a smoothing step and a coarse
level solution step as described in detail in Section 6.

###  3.3 Motivating guidelines for the design of the coarse functional $g_{c}$

One choice is to define

$g_{c}(\mathbf{x}_{c})=g({\mathbf{P}}\mathbf{x}_{c}),$ (11)

which admits excellent approximation but is non-optimal in terms of cost. Let
us examine the choice (11) to motivate an approximation. Assume that
${\mathbf{P}}$ and ${\boldsymbol{\Pi}}$ are such that
${\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}$ approximates $\mathbf{x}$, i.e.,
$\|\mathbf{x}-{\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}\|$ is “small” in some
norm $\|.\|$. Let ${\mathbf{P}}$, represented as a matrix, admit the form
${\mathbf{P}}=(p_{i,i_{c}})$, $1\leq i\leq n$, $1\leq i_{c}\leq n_{c}$.

Under the coarse functional choice (11), we have

$\partial_{j_{c}}g_{c}(\mathbf{x}_{c})=\sum\limits^{n}_{j=1}\partial_{j}g({\mathbf{P}}\mathbf{x}_{c})p_{j,j_{c}},$

and similarly

$\partial_{i_{c}}\partial_{j_{c}}g_{c}(\mathbf{x}_{c})=\sum\limits^{n}_{i,j=1}\partial_{i}\partial_{j}g({\mathbf{P}}\mathbf{x}_{c})p_{j,j_{c}}p_{i,i_{c}}.$

In other words, we have

$H_{g_{c}}(\mathbf{x}_{c})={\mathbf{P}}^{T}H_{g}({\mathbf{P}}\mathbf{x}_{c}){\mathbf{P}}.$

Now, suppose $\mathbf{x}_{c}$ is derived from the restriction
$\mathbf{x}_{c}={\boldsymbol{\Pi}}\mathbf{x}$. Then, we have

$H_{g_{c}}({\boldsymbol{\Pi}}\mathbf{x})={\mathbf{P}}^{T}H_{g}({\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}){\mathbf{P}}={\mathbf{P}}^{T}H_{g}(\mathbf{x}){\mathbf{P}}+{\mathbf{P}}^{T}\left(H_{g}({\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x})-H_{g}(\mathbf{x})\right){\mathbf{P}}.$

Because $\|\mathbf{x}-{\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}\|$ is small,
the second term in this formula is small, so

$H_{g_{c}}({\boldsymbol{\Pi}}\mathbf{x})\approx{\mathbf{P}}^{T}H_{g}(\mathbf{x}){\mathbf{P}},$

so that condition (8) is approximately satisfied.

The exact choice of $g_{c}$ is problem-dependent, and we design such a coarse
functional for neural networks below, but this motivates some guidelines:

  1. (A)

It is desirable that $g_{c}$ be chosen so that
$H_{g_{c}}(\mathbf{x}_{c})\approx{\mathbf{P}}^{T}H_{g}({\mathbf{P}}\mathbf{x}_{c}){\mathbf{P}}$

  2. (B)

It is desirable that $\mathbf{R}={\mathbf{P}}^{T}$.

  3. (C)

It is desirable that ${\mathbf{P}}$ and ${\boldsymbol{\Pi}}$ be chosen so that
$\|{\mathbf{P}}^{T}\left(H_{g}({\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x})-H_{g}(\mathbf{x})\right){\mathbf{P}}\|$
is small. As the Hessian matrix is a continuous function of $\mathbf{x}$, this
guideline can be simplified to a desire that
$\|\mathbf{x}-{\mathbf{P}}{\boldsymbol{\Pi}}\mathbf{x}\|_{2}$ be small.

These properties are desirable because, when they are true, Theorem 3.1
approximately holds, ensuring an approximation for first-order consistency.

##  4 Application of Full Approximation Scheme to Deep Neural Network training

In the previous section, we derived some useful properties of two-level FAS
for optimization along with a set of advice to help guide the development of
${\mathbf{P}}$ and ${\boldsymbol{\Pi}}$ operators. In this section, we apply
those results to the problem of training deep neural networks via FAS.

A feedforward DNN (deep neural network) is defined as follows:

  * •

A set of $n_{L}$ layers, where the $k$th layer is associated with a vector
variable $\mathbf{y}_{k}$.

  * •

Two consecutive layers are connected by an activation function $\varphi$ which
acts on each component of the vectors involved.

  * •

Also, two consecutive layers, $k$ and $k+1$, are connected by a transition
coefficient matrix $W_{k}$, and a shift vector $\mathbf{b}_{k}$.

  * •

To compute the parameters $\\{W_{k}\\}^{n_{L}}_{k=0}$ and
$\\{\mathbf{b}_{k}\\}^{n_{L}}_{k=0}$, the DNN makes use of two sets of given
quantities, $\\{\mathbf{y}^{in}_{s}\\}$ \- input, and
$\\{\mathbf{y}^{out}_{s}\\}$\- the respective output.

The DNN is _trained_ by computing the weights $\\{W_{k}\\}$ and
$\\{\mathbf{b}_{k}\\}$, such that for the set of given inputs
$\\{\mathbf{y}^{in}_{s}\\}$, it produces an accurate approximations to the
respective given set of outputs $\\{\mathbf{y}^{out}_{s}\\}$, i.e., for each
$s$, $DNN(\mathbf{y}^{in}_{s})\approx\mathbf{y}^{out}_{s}$. The DNN parameters
$\\{W_{k}\\}$ and $\\{\mathbf{b}_{k}\\}$, are computed by solving a suitable
optimization problem.

For any given set of parameters $W=\\{W_{k}\\}$ and
$\mathbf{b}=\\{\mathbf{b}_{k}\\}$, let $\mathbf{x}$ be a vector quantity that
unrolls and concatenates them, that is, let

$\mathbf{x}=\left[W_{0}[:],W_{1}[:],\dots,\;W_{n_{L}}[:],\mathbf{b}_{0}[:],\mathbf{b}_{1}[:],\dots,\mathbf{b}_{n_{L}}[:]\right]$
(12)

be the unrolled vector of all learnable parameters in the DNN. Given
$\mathbf{x}$, the DNN acts on any input vector data $\mathbf{y}^{in}$, i.e.,
computes
$\mathbf{y}^{in}\mapsto\mathbf{y}^{out}:=DNN(\mathbf{x};\;\mathbf{y}^{in})$,
by evaluating the recurrence relation

$\begin{array}[]{rl}\mathbf{y}_{1}&=\varphi\left(W_{0}\mathbf{y}^{in}+\mathbf{b}_{0}\right),\\\
\mathbf{y}_{k+1}&=\varphi\left(W_{k}\mathbf{y}_{k}+\mathbf{b}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\mathbf{y}^{out}&=\varphi\left(W_{n_{L}}\mathbf{y}_{n_{L}}+\mathbf{b}_{n_{L}}\right).\end{array}$
(13)

######  Remark 4.1.

The above definition extends to convolutional neural networks (CNNs) as well,
using Toeplitz structured $W_{k}$ matrices and $b_{k}$ vectors. We discuss
this and the application of FAS to CNNs in Section 7.

Training the network selects the weight matrices $W=\\{W_{k}\\}$ and bias
vectors $\mathbf{b}=\\{\mathbf{b}_{k}\\}$ by approximately solving the
optimization problem

$g(\mathbf{x})=g(W,\mathbf{b})=\sum_{s=1}^{N}L_{s}(DNN(\mathbf{x};\;\mathbf{y}^{in}_{s})),$
(14)

where $L_{s}$ is some loss function associated with the $s$th training sample.

One example of a _loss functional_ is the mean-squared error

$g(\mathbf{x})=\frac{1}{N}\sum\limits^{N}_{s=1}\|\mathbf{y}^{out}_{s}-DNN(\mathbf{x};\;\mathbf{y}^{in}_{s})\|^{2}.$
(15)

In the present setting, to train the DNN means to find a set of parameters
$\mathbf{x}=(W,\;\mathbf{b})$ that approximately minimize the above functional
$g(\mathbf{x})$. However, as discussed in Section 2.2.1, one does not
generally wish to exactly minimize this functional due to overfitting
concerns. Rather, one attempts to find a set of parameters that approximately
minimize $g(\mathbf{x})$, and achieve a small but nonzero gradient.

The popularity of the DNNs is due to the fact that the gradients of the loss
functional $\nabla g(\mathbf{x})$ are readily computable by a
_backpropagation_ , and a number of popular libraries such as TensorFlow [20]
and PyTorch[1] readily compute them. Moreover, in the optimization algorithms
employed in the DNN training we may use only a portion of the sum in (15) (or
(14)) by employing randomness. That is, for a _minibatch_ of $m$ randomly
selected indices $\mathbf{s}=(s_{1},s_{2},\;\dots,\;s_{m})$, we compute the
gradient of

$g_{\mathbf{s}}(\mathbf{x})=\frac{1}{m}\sum\limits^{m}_{j=1}L_{s_{j}}(DNN(\mathbf{x};\;\mathbf{y}^{in}_{s_{j}})).$
(16)

This partial (stochastic) gradient, $\nabla g_{\mathbf{s}}(\mathbf{x})$, is
used to update the iterates $\mathbf{x}:=\mathbf{x}-\tau\nabla
g_{\mathbf{s}}(\mathbf{x})$ for a chosen step length $\tau\in(0,1]$, in a
gradient descent algorithm to minimize $g(\mathbf{x})$.

##  5 Restriction, Interpolation and $\tau$-correction operators for the DNN
training

In this section, we motivate and define our restriction and interpolation
operators, as well as the $\tau$-correction for the case of training DNNs, in
a manner consistent with the guidelines of Section 3.3.

First, we clarify notation. Consider a DNN with the vector quantity
$\mathbf{x}$ containing its learnable parameters as in Equation (12). To
define a two-level FAS algorithm we need an interpolation mapping
${\mathbf{P}}$ and a fine-to-coarse restriction mapping ${\boldsymbol{\Pi}}$
such that $\mathbf{Q}={\mathbf{P}}{\boldsymbol{\Pi}}$ is a projection. We will
choose ${\mathbf{P}}$ and ${\boldsymbol{\Pi}}$ such that
${\boldsymbol{\Pi}}{\mathbf{P}}=I$.

As the parameters of a DNN are defined layerwise, we will also decompose the
${\mathbf{P}}$ and ${\boldsymbol{\Pi}}$ operators in this manner. Consider the
sets $P=\\{P_{k}\\}^{n_{L}+1}_{k=0}$ and $\pi=\\{\pi_{k}\\}^{n_{L}+1}_{k=0}$,
where $P_{k}$ and $\pi_{k}$ act on vectors associated with the $k$th layer;
that is, given a data vector $\mathbf{y}_{k}$ at layer $k$,
$\pi_{k}\mathbf{y}_{k}$ converts $\mathbf{y}_{k}$ to a coarsened
$\mathbf{y}_{k}^{c}$. We set $\pi_{in}=\pi_{0}$, $P_{in}=P_{0}$ and
$\pi_{out}=\pi_{n_{L}+1},P_{out}=P_{n_{L}+1}$. We also assume that for each
$k$, $\pi_{k}$ and $P_{k}$ make a projection, that is,

$\pi_{k}P_{k}=I.$ (17)

Given a set of coarse parameters $\\{W^{c}_{k}\\}$ and
$\\{\mathbf{b}^{c}_{k}\\}$, neural networks in fact provide a natural way to
produce a coarse functional $g_{c}$ that satisfies Guideline (A) of Section
3.3. We simply define a coarse DNN as a coarse version of the same recursion
as in (13):

$\begin{array}[]{rl}\mathbf{y}^{c}_{1}&=\varphi\left(W^{c}_{0}\mathbf{y}^{in}_{c}+\mathbf{b}^{c}_{0}\right),\\\
\mathbf{y}^{c}_{k+1}&=\varphi\left(W^{c}_{k}\mathbf{y}^{c}_{k}+\mathbf{b}^{c}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\mathbf{y}^{out}_{c}&=\varphi\left(W^{c}_{n_{L}}\mathbf{y}_{n_{L}}+\mathbf{b}^{c}_{n_{L}}\right).\end{array}$
(18)

That is, we have
$\mathbf{y}^{out}_{c}=DNN_{c}(\mathbf{x}_{c};\;\mathbf{y}^{in}_{c})$, where
$\mathbf{x}_{c}$ combines all parameters $\\{W^{c}_{k}\\}$ and
$\\{\mathbf{b}^{c}_{k}\\}$ as in Equation (12).

The training objective function is a summation over a coarse set of training
samples:

$g_{c}(\mathbf{x}_{c})=g_{c}(W^{c},b^{c})=\sum_{s=1}^{N}L^{c}_{s}(DNN_{c}(\mathbf{x}_{c};\;\mathbf{y}^{in}_{c,s})).$
(19)

If we do not have coarse counterparts of the training samples
$\mathbf{y}^{in}_{c}$ and $\mathbf{y}^{out}_{c}$, we use the only given (fine
level) set of samples, i.e., we let
$\mathbf{y}^{in}_{c,s}=\mathbf{y}^{in}_{s}$ and
$\mathbf{y}^{out}_{c,s}=\mathbf{y}^{out}_{s}$. This implies that we let
$\pi_{0}=I,P_{0}=I$ and $\pi_{n_{L}+1}=I,P_{n_{L}+1}=I$. Also, in that case we
have $L^{c}_{s}=L_{s}$. Because $g_{c}$ is designed to approximate $g$ here,
Guideline (A) is satsified.

In what follows, we construct a coarse DNN from a fine DNN parameters; and
vice versa, given a coarse DNN, we embed it in a finer DNN.

###  5.1 Fine-to-coarse DNN restriction

Let us assume that $\|\mathbf{y}_{k}-P_{k}\pi_{k}\mathbf{y}_{k}\|$ is small in
the sense discussed in Section 3.3. Then the following is an approximation to
the standard DNN recurrence (13):

$\begin{array}[]{rl}\pi_{1}\mathbf{y}_{1}&=\pi_{1}\varphi\left(W_{0}P_{0}\pi_{in}\mathbf{y}^{in}+\mathbf{b}_{0}\right),\\\
\pi_{k+1}\mathbf{y}_{k+1}&=\pi_{k+1}\varphi\left(W_{k}P_{k}\pi_{k}\mathbf{y}_{k}+\mathbf{b}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\pi_{out}\mathbf{y}^{out}&=\pi_{out}\varphi\left(W_{n_{L}}P_{n_{L}}\pi_{n_{L}}\mathbf{y}_{n_{L}}+\mathbf{b}_{n_{L}}\right).\end{array}$
(20)

By letting $\mathbf{y}^{c}_{k}=\pi_{k}\mathbf{y}_{k}$, we have the following
recurrence involving coarse vectors $\mathbf{y}^{c}_{k}$,

$\begin{array}[]{rl}\mathbf{y}^{c}_{1}&=\pi_{1}\varphi\left(W_{0}P_{0}\mathbf{y}^{in}_{c}+\mathbf{b}_{0}\right),\\\
\mathbf{y}^{c}_{k+1}&=\pi_{k+1}\varphi\left(W_{k}P_{k}\mathbf{y}^{c}_{k}+\mathbf{b}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\mathbf{y}^{out}_{c}&=\pi_{n_{L}+1}\varphi\left(W_{n_{L}}P_{n_{L}}\mathbf{y}^{c}_{n_{L}}+\mathbf{b}_{n_{L}}\right).\end{array}$
(21)

Finally, we assume the approximation

$\varphi(\pi\mathbf{z})\approx\pi\varphi(\mathbf{z}).$ (22)

We discuss in Section 5.4 when this is valid. Given this, recurrence (21) can
be approximated as

$\begin{array}[]{rl}\mathbf{y}^{c}_{1}&=\varphi\left(\pi_{1}W_{0}P_{0}\mathbf{y}^{in}_{c}+\pi_{1}\mathbf{b}_{0}\right),\\\
\mathbf{y}^{c}_{k+1}&=\varphi\left(\pi_{k+1}W_{k}P_{k}\mathbf{y}^{c}_{k}+\pi_{k+1}\mathbf{b}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\mathbf{y}^{out}_{c}&=\varphi\left(\pi_{n_{L}+1}W_{n_{L}}P_{n_{L}}\mathbf{y}^{c}_{n_{L}}+\pi_{n_{L}}\mathbf{b}_{n_{L}}\right).\end{array}$

That is, letting

$W^{c}_{k}=\pi_{k+1}W_{k}P_{k}\text{ and
}\mathbf{b}^{c}_{k}=\pi_{k+1}\mathbf{b}_{k},$ (23)

we end up with the recurrence that implements the actions of a coarse DNN,
$\mathbf{y}^{in}_{c}\mapsto\mathbf{y}^{out}_{c}=DNN_{c}(\mathbf{x}_{c};\;\mathbf{y}^{out}_{c})$:

$\begin{array}[]{rl}\mathbf{y}^{c}_{1}&=\varphi\left(W^{c}_{0}\mathbf{y}^{in}_{c}+\mathbf{b}^{c}_{0}\right),\\\
\mathbf{y}^{c}_{k+1}&=\varphi\left(W^{c}_{k}\mathbf{y}^{c}_{k}+\mathbf{b}^{c}_{k}\right),\;\text{
for }0\leq k\leq n_{L}-1,\\\
\mathbf{y}^{out}_{c}&=\varphi\left(W^{c}_{n_{L}}\mathbf{y}^{c}_{n_{L}}+\mathbf{b}^{c}_{n_{L}}\right).\end{array}$
(24)

This motivates our fine-to-coarse mapping, ${\boldsymbol{\Pi}}$, which we
define as

$\begin{array}[]{rl}\mathbf{x}=&[(W_{k}),\;(\mathbf{b}_{k})],\\\
{\boldsymbol{\Pi}}\mathbf{x}=&[(\pi_{k+1}W_{k}P_{k}),\;(\pi_{k+1}\mathbf{b}_{k})].\end{array}$
(25)

Note that, despite the layer-wise decomposition, this is still a linear
mapping of $\mathbf{x}$, so there exists a matrix ${\boldsymbol{\Pi}}$ whose
action is this.

###  5.2 Coarse-to-fine DNN prolongation

Consider the coarse DNN recurrence (24). We assume the approximation

$P\varphi(\mathbf{z}_{c})\approx\varphi(P\mathbf{z}_{c}),$ (26)

discussed in Section 5.4, and using the property $\pi_{k}P_{k}=I$, we have

$\begin{array}[]{rll}P_{1}\mathbf{y}^{c}_{1}&=P_{1}\varphi\left(W^{c}_{0}\pi_{0}(P_{0}\mathbf{y}^{in}_{c})+\mathbf{b}^{c}_{0}\right)&\approx\varphi\left(P_{1}W^{c}_{0}\pi_{0}(P_{0}\mathbf{y}^{in}_{c})+P_{1}\mathbf{b}^{c}_{0}\right),\\\
P_{k+1}\mathbf{y}^{c}_{k+1}&=P_{k+1}\varphi\left(W^{c}_{k}\pi_{k}(P_{k}\mathbf{y}^{c}_{k})+\mathbf{b}^{c}_{k}\right)&\approx\varphi\left(P_{k+1}W^{c}_{k}\pi_{k}(P_{k}\mathbf{y}^{c}_{k})+P_{k+1}\mathbf{b}^{c}_{k}\right),\;k\leq
n_{L}-1,\\\
P_{out}\mathbf{y}^{out}_{c}&=P_{out}\varphi\left(W^{c}_{n_{L}}\pi_{n_{L}}(P_{n_{L}}\mathbf{y}^{c}_{n_{L}})+\mathbf{b}^{c}_{n_{L}}\right)&\approx\varphi\left(P_{n_{L}+1}W^{c}_{n_{L}}\pi_{n_{L}}(P_{n_{L}}\mathbf{y}^{c}_{n_{L}})+P_{n_{L}+1}\mathbf{b}^{c}_{n_{L}}\right).\end{array}$

The latter formulas motivates us to define the interpolation mapping
${\mathbf{P}}$ as follows. Let
$\mathbf{x}_{c}=[(W^{c}_{k}),\;(\mathbf{b}^{c}_{k})]$, then

${\mathbf{P}}\mathbf{x}_{c}=[(P_{k+1}W^{c}_{k}\pi_{k}),\;(P_{k+1}\mathbf{b}^{c}_{k})].$
(27)

Note that this does not require that the $W_{k}$ matrices be square. It only
requires that the input dimension to layer $k$ match the output dimension of
layer $k-1$.

###  5.3 The restriction mapping $\mathbf{R}={\mathbf{P}}^{T}$

The definition of $\mathbf{R}$ is already given by (9). The layerwise
definition is

$\mathbf{R}\mathbf{x}=[(P^{T}_{k+1}W_{k}\pi^{T}_{k}),\;(P^{T}_{k}\mathbf{b}_{k})].$
(28)

The following result, satisfying Guideline (B) of Section 3.3, holds.

######  Lemma 5.1.

The restriction mapping $\mathbf{R}$ defined in (28) and the interpolation
mapping ${\mathbf{P}}$ defined in (27) are transpose to each other, i.e.,

$\mathbf{R}={\mathbf{P}}^{T}.$

###### Proof.

We start with the inner product of $\mathbf{x}=[(W_{k}),\;(\mathbf{b}_{k})]$
and
${\mathbf{P}}\mathbf{x}_{c}={\mathbf{P}}[(W^{c}_{k}),\;(\mathbf{b}^{c}_{k})]$,

$(\mathbf{x}_{c})^{T}({\mathbf{P}}^{T}\mathbf{x})=\mathbf{x}^{T}({\mathbf{P}}\mathbf{x}_{c})=\sum\limits_{i}x_{i}({\mathbf{P}}\mathbf{x}_{c})_{i}.$

With $\mathbf{b}_{k}=(b_{k,i})$, $\mathbf{b}^{c}_{k}=(b^{c}_{k,i_{c}})$,
$W_{k}=\left((W_{k})_{r,s}\right)$ and
$W^{c}_{k}=\left((W^{c}_{k})_{r_{c},s_{c}}\right)$, based on the Definition
(27), we have the representation

$\mathbf{x}^{T}({\mathbf{P}}\mathbf{x}_{c})=\sum\limits_{i}x_{i}({\mathbf{P}}\mathbf{x}_{c})_{i}=\sum\limits_{k}\sum\limits_{i}b_{k,i}(P_{k}\mathbf{b}^{c}_{k})_{i}+\sum\limits_{k}\sum\limits_{i=(r,s)}(W_{k})_{r,s}(P_{k+1}W^{c}_{k}\pi_{k})_{r,s}.$
(29)

The first term on the right gives

$\sum\limits_{k}\sum\limits_{i}b_{k,i}(P_{k}\mathbf{b}^{c}_{k})_{i}=\sum\limits_{k}(\mathbf{b}^{c}_{k})^{T}P^{T}_{k}\mathbf{b}_{k}=\sum\limits_{k}(\mathbf{b}^{c}_{k})^{T}\mathbf{b}^{c}_{k}.$
(30)

For the second term, by rearranging terms and using the definition (28) of
$\mathbf{R}\mathbf{x}$, we have

$\begin{array}[]{rl}\sum\limits_{k}\sum\limits_{i=(r,s)}(W_{k})_{r,s}(P_{k+1}W^{c}_{k}\pi_{k})_{r,s}&=\sum\limits_{k}\sum\limits_{i=(r,s)}\sum\limits_{i_{c}=(r_{c},s_{c})}(W_{k})_{r,s}(P_{k+1})_{r,r_{c}}(W^{c}_{k})_{r_{c},s_{c}}(\pi_{k})_{s_{c},s}\\\
&=\sum\limits_{k}\sum\limits_{i_{c}=(r_{c},s_{c})}(W^{c}_{k})_{i_{c}}\sum\limits_{r,s}(P_{k+1})_{r,r_{c}}(W_{k})_{r,s}(\pi_{k})_{s_{c},s}\\\
&=\sum\limits_{k}\sum\limits_{i_{c}=(r_{c},s_{c})}(W^{c}_{k})_{i_{c}}(P^{T}_{k+1}W_{k}\pi^{T}_{k})_{i_{c}}\\\
&=\sum\limits_{k}\sum\limits_{i_{c}}(W^{c}_{k})_{i_{c}}(W^{c}_{k})_{i_{c}}.\end{array}$
(31)

Thus substituting (30) and (31) in (29), using formulas (28) defining the
components of $\mathbf{R}\mathbf{x}$, we arrive at

$(\mathbf{x}_{c})^{T}({\mathbf{P}}^{T}\mathbf{x})=\mathbf{x}^{T}({\mathbf{P}}\mathbf{x}_{c})=(\mathbf{x}_{c})^{T}\mathbf{R}\mathbf{x},$

which proves the desired result. ∎

######  Remark 5.1.

We note that the linear mappings ${\boldsymbol{\Pi}}$, ${\mathbf{P}}$, as well
as $\mathbf{R}={\mathbf{P}}^{T}$ are not explicitly formed as matrices.
However, from their definitions it is clear that their actions are readily
available via the actions of the layer matrices $P_{k}$, $\pi_{k}$ and their
transposes. That is, in an implementation, only the layer matrices $\pi_{k}$
and $P_{k}$ need to be constructed and stored explicitly as (sparse) matrices.

######  Lemma 5.2.

With the definitions (25) of ${\boldsymbol{\Pi}}$ and (27) of ${\mathbf{P}}$,
assuming (17) for the layer matrices $\pi_{k}$ and $P_{k}$, the following main
property holds

${\boldsymbol{\Pi}}{\mathbf{P}}=I,$

which implies that $\mathbf{Q}={\mathbf{P}}{\boldsymbol{\Pi}}$ is a projection
($\mathbf{Q}^{2}=\mathbf{Q}$).

###### Proof.

Given $\mathbf{x}_{c}=[(W^{c}_{k}),\;(\mathbf{b}^{c}_{k})]$, define
$\mathbf{x}=[(W_{k}),\;(\mathbf{b}_{k})]={\mathbf{P}}\mathbf{x}_{c}$, i.e.,
$W_{k}=P_{k+1}W^{c}_{k}\pi_{k}$ and
$\mathbf{b}_{k}=P_{k+1}\mathbf{b}^{c}_{k}$. Then for
${\boldsymbol{\Pi}}\mathbf{x}=[(\pi_{k+1}W_{k}P_{k}),\;(\pi_{k+1}\mathbf{b}_{k})]$
since $\pi_{j}P_{j}=I$, we have

$\begin{array}[]{rl}\pi_{k+1}W_{k}P_{k}&=\pi_{k+1}P_{k+1}W^{c}_{k}\pi_{k}P_{k}=W^{c}_{k},\\\
\pi_{k+1}\mathbf{b}_{k}&=\pi_{k+1}P_{k+1}\mathbf{b}^{c}_{k}=\mathbf{b}^{c}_{k}.\end{array}$

That is,
${\boldsymbol{\Pi}}\mathbf{x}={\boldsymbol{\Pi}}{\mathbf{P}}\mathbf{x}_{c}=\mathbf{x}_{c}$,
which completes the proof. ∎

###  5.4 A simple choice for $\pi_{k}$ and $P_{k}$

Finally, we must choose ${\boldsymbol{\Pi}}$ and ${\mathbf{P}}$ to satisfy
Guideline (C) of Section 3.3. Intuitively, if two neurons have nearly the same
weights and biases, then we could average them together to construct a coarse
neuron that approximates the two fine neurons. Therefore, we choose $\pi_{k}$
to perform pairwise averaging and $P_{k}$ to perform piecewise constant
interpolation over selected pairwise aggregates. The choice of the pairwise
aggregates is done through a Heavy Edge Matching (HEM) method [14]. The
mappings are simply written as follows

$P_{k}=\begin{bmatrix}1\\\ 1\\\ &1\\\ &1\\\
&&\ddots\end{bmatrix}\in\mathbb{R}^{n_{k}\times
n^{c}_{k}},\quad\pi_{k}=\begin{bmatrix}1/2&1/2\\\ &&1/2&1/2\\\
&&&&\ddots\end{bmatrix}\in\mathbb{R}^{n^{c}_{k}\times n_{k}}.$ (32)

When the activation function $\varphi$ is the commonly-used ReLU function,

$\varphi(x)=\left\\{\begin{array}[]{rl}x,&x\geq 0\\\
0,&x<0\end{array}\right.,$ (33)

these choices for $P_{k}$ and $\pi_{k}$ justify the approximations (22) and
(26). In particular, (22) is valid because the HEM (discussed in Section 5.5
and shown in Algorithm 1) matches together neurons whose parameter vectors
point in approximately the same direction, and as a result, tend to produce
positive and negative dot products with input vectors at the same time. The
nonlinearity of the ReLU depends entirely on that positive/negative property
of its input, and so (22) will be exactly true except when dot products with
paired neurons produce different signs. Similarly, (26) is in fact an
equality, by design, due to the piecewise-constant prolongation. Other
activation functions are likely effective here as well; in general, if we
approximate $\varphi(\mathbf{z})$ near $\mathbf{z}$ linearly with a _matrix_
$\widehat{\varphi}$, then (26) holds when
$\pi\widehat{\varphi}\approx\widehat{\varphi}\pi$. In particular, other
piecewise-linear activation functions such as the ”leaky ReLU” satisfy this.

A variant to (32) that is often found better is to scale $P_{k}$ by the row
norms of $W_{k}$, i.e.,

$\tilde{P}_{k}=D_{k}P_{k},\quad\tilde{\pi}_{k}=(P_{k}^{T}D_{k}P_{k})^{-1}P_{k}^{T},$
(34)

where $D_{k}$ is the diagonal matrix that contains the row norms of $W_{k}$,
so that the rows that correspond to the paired neurons in an aggregate are
invariant with the projector $Q$ if they are linear dependent, i.e., they have
the maximum similarity score in the HEM algorithm. Apparently, (32) is a
special case where $D_{k}=I$. Furthermore, $P_{k}^{T}D_{k}P_{k}$ is a diagonal
matrix due to the structures of $P_{k}$ and $D_{k}$.

######  Remark 5.2.

Although we write Equations (32) and (34) as matrices, these are very sparse
matrices and the _actions_ of these matrices can be efficiently computed
without ever directly forming $\tilde{P}_{k}$ and $\tilde{\pi}_{k}$. In our
software, we have implemented matrix-free versions for fast, practical
computation.

###  5.5 HEM algorithm for coarsening

Standard HEM algorithms have been used as coarsening algorithms for
aggregation-based AMG, see, e.g., [23]. Finding the optimal HEM of a graph is
in general difficult, but greedy algorithms can usually work well. A sketch of
greedy HEM is given in Algorithm 1, where input $S_{k}$ is the strength-of-
connection matrix, which defines the similarities between neurons. Cosine
angles of the rows of $W_{k}$ are used for $S_{k}$. In Algorithm 1, neurons
are visited in a given order and at each time if neuron $i$ has not been
matched, neuron $j$ that has the heaviest edge between $i$ and $j$ among all
the remaining unmatched neurons is selected. Edges that have weights smaller
than a predefined threshold $\theta$ are disregarded in the matching. The
results are kept in $\pi$ and $\zeta$, where $\pi(i)$ keeps the matched neuron
of $i$, and $\zeta(i)$ stores the index of the aggregate of $i$ and $\pi(i)$.
Clearly, on exit of the algorithm, we have $\pi(\pi(i))=i$ (if $i=\pi(i)$, $i$
is a singleton), and $\zeta(i)=\zeta(\pi(i))$. Roughly speaking, a coarsening
factor close to 2 is typically obtained. The HEM method can also be
parallelized for faster processing on GPUs via parallel matching methods such
as [22].

Algorithm 1 A greedy HEM algorithm

1:Input: Strength-of-connection matrix $S_{k}$ and strength threshold $\theta$

2:Output: matching $\pi$, aggregation $\zeta$, and the number of aggregates
$n_{k}^{c}$

3:for $i=1,2,\ldots,n_{k}$ do

4: $\pi(i):=0$

5:end for

6:$n_{k}^{c}:=0$

7:for $i=1,2,\ldots,n_{k}$ do

8: if $\pi(i)=0$ then

9: if $\exists j\neq i:\pi(j)=0\land S_{ij}>\theta\land\forall k:S_{ij}\geq
S_{ik},\pi(k)=0$ then

10: $\pi(i):=j,\;\pi(j):=i$

11: $n_{k}^{c}:=n_{k}^{c}+1,\;\zeta(i):=n_{k}^{c},\;\zeta(j):=n_{k}^{c}$

12: else$\triangleright$ match $i$ to itself

13: $\pi(i):=i,\;n_{k}^{c}:=n_{k}^{c}+1,\;\zeta(i):=n_{k}^{c}$

14: end if

15: end if

16:end for

###  5.6 Constructing a tau-correction

As discussed above, the tau correction in optimization is $\tau=\nabla
g_{c}({\boldsymbol{\Pi}}\mathbf{x})-\mathbf{R}\nabla g(x)$, where
$\mathbf{R}={\mathbf{P}}^{T}$, and ${\boldsymbol{\Pi}}$ are the two linear
fine-to-coarse mappings,(28)-(25), defined in our DNN setting. Note that the
operation $\mathbf{R}\nabla g(\mathbf{x})$ is well-defined, as $\nabla
g(\mathbf{x})$ has the same structure as $\mathbf{x}$.

One wrinkle is that the calculation of a gradient requires a summation over
all training samples in the training set. Often, training sets have a large
number of samples in them, so computing the true gradient would be
computationally unwieldy. We approximate this by simply calculating a
stochastic tau-correction, in which we average over a small subset
$\mathbf{s}=(s_{1},\;s_{2},\;\dots,\;s_{m})$ of the total training set, and
then scale appropriately:

$\begin{array}[]{rl}\nabla_{\mathbf{x}}g_{\mathbf{s}}(\mathbf{x})&\displaystyle=\sum_{j=1}^{m}\nabla_{\mathbf{x}}L_{s_{j}}(DNN(\mathbf{x},\;y^{in}_{s_{j}})\\\
\nabla_{\mathbf{x}_{c}}g^{c}_{\mathbf{s}}(\mathbf{x}_{c})&\displaystyle=\sum_{j=1}^{m}\nabla_{\mathbf{x}_{c}}L^{c}_{s_{j}}(DNN_{c}(\mathbf{x}_{c},y^{in}_{c,s_{j}})\\\
\mathbf{\tau}&\displaystyle=\frac{N}{m}\left(\left(\nabla_{\mathbf{x}_{c}}g^{c}_{\mathbf{s}}\right)({\boldsymbol{\Pi}}\mathbf{x})-\mathbf{R}\nabla_{\mathbf{x}}g_{\mathbf{s}}(\mathbf{x})\right),\end{array}$
(35)

where $m$ is the size of a tau batch, $s_{j}$ is the index of the $j$’th
sample in the tau batch, $\nabla_{\mathbf{y}}$ stands for a gradient with
respect to variable $\mathbf{y}$, and $N$ is the number of minibatches in the
dataset as a whole.

Due to this stochastic approximation, we in general will never construct a
coarse problem such that $\nabla\widehat{g}(\mathbf{x}_{c})=0$ exactly. This
is consistent with how neural networks are trained: Specifically, because
traditional training is performed in a stochastic manner over minibatches,
exact gradients are never computed, and one never achieves an $\mathbf{x}$
such that $\nabla g(\mathbf{x})=0$ exactly, or even iteratively converges to
such a point. In fact, this behavior is often considered desirable, see
Section 2.2.1.

##  6 The DNN FAS Algorithm

The original objective function $g$ is defined in Equation (14), and the
auxiliary objective function in Equation (19). Given a tau correction, we can
construct a tau-corrected auxiliary problem like in Equation (3).

We recall that the FAS step cannot be our only method for improving
$\mathbf{x}$. As discussed earlier, Remark 3.1, we need a correction in a
space complementary to the coarse space $\text{Range}({\mathbf{P}})$, a step
customarily referred to as _smoothing_. That is because $\mathbf{x}_{c}$
exists in a reduced-dimension space and so cannot represent all the possible
changes we might make.

There exist a wide variety of methods one might choose for the smoothing step.
In general, it should be cost-efficient and effectively target the complement
to $\text{Range}({\mathbf{P}})$. The neural network community has already
provided a number of efficient first-order training methods which can be
leveraged for smoothing. We choose here to use stochastic gradient descent
(SGD). Our reason for using this is that SGD in this context is quite
analogous to the commonly-used Jacobi smoother in traditional multigrid
methods. Indeed, when applying the Full Approximation Scheme to the
deterministic problem of minimizing quadratic forms, (deterministic) gradient
descent is exactly a Richardson iteration, very similar to the Jacobi
iteration commonly used as a multigrid smoother [5]. A two-level FAS V-cycle,
then, consists of the following steps:

  1. 1.

Apply one or more iterations of SGD to the original problem
$g(\mathbf{x})=g(W,\mathbf{b})$.

  2. 2.

Calculate the stochastic tau correction $\mathbf{\tau}$.

  3. 3.

Apply the restriction operation to calculate the elements of the coarse neural
network $W^{c},\mathbf{b}^{c}$.

  4. 4.

Apply one or more iterations of SGD to the corrected auxiliary problem
${\widehat{g}}(\mathbf{x}_{c})={\widehat{g}}(W^{c},\mathbf{b}^{c})$ (see (3)).

  5. 5.

Apply the prolongation operation
${\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})$ and then the FAS
update
$\mathbf{x}\leftarrow\mathbf{x}_{new}=\mathbf{x}+{\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})$.

  6. 6.

Apply one or more iterations of SGD to the original problem
$g(\mathbf{x})=g(W,b)$.

FAS V-cycles of greater than two levels are achieved through recursion at Step
4.

###  6.1 A Few Notes on Algorithmic Tuning

####  6.1.1 Minibatch Selection

We perform minibatch selection in a manner intended to closely match that of
traditional one-level training. We first shuffle the entire training dataset,
and then proceed by cycling through the dataset. When computing the stochastic
tau correction, we draw the next $m$ minibatches in shuffled order, compute
$\mathbf{\tau}$ over those $m$ minibatches as in Equation (35), and then use
those same $m$ minibatches for SGD at the coarse level. A new set of $m$
minibatches are drawn when restricting to the next-coarser level. Because SGD
at the coarse level is performed using the same minibatches as were used for
the tau-correction, one might view this as a _stochastic V-cycle_. Because
multigrid V-cycles typically only apply a small number of smoothing steps at
each level, we tend to keep $m$ fairly small as well, often between 2 and 8.

####  6.1.2 Momentum

Momentum [25] is an important component in training methods, as it helps to
smooth out variations and move parameters more consistently in directions that
matter. In fact, SGD with momentum (and many traditional training methods) use
gradient information to update momentum vectors, and only directly update
parameter vectors through momentum. In our setting, important information
about desirable directions of movement is learned at the coarse levels, and we
should not throw this important information away upon prolongation. Therefore,
we modify this approach slightly to restrict and prolong _momentum_ as well as
network parameters. As there is a momentum parameter for each network
parameter, a momentum vector $\mathbf{m}$ has the same structure as a
parameter vector $\mathbf{x}$.

Upon restriction, we also restrict the momentum vector using the same
restriction operator $\mathbf{m}_{c}={\boldsymbol{\Pi}}\mathbf{m}$, and upon
prolongation, we perform the same coarse-grid correction

$\mathbf{m}=\mathbf{m}+{\mathbf{P}}(\mathbf{m}_{c}-{\boldsymbol{\Pi}}\mathbf{m}).$
(36)

We show the two-level FAS algorithm, with multilevel momentum, in Algorithm 2.

####  6.1.3 HEM Re-Matching

The HEM matching in Algorithm 1 depends on the weights, which change during
the training process. Therefore, matchings tend to change somewhat throughout
the training process. We find that it is generally beneficial to periodically
re-compute matchings, but that the training is not particularly sensitive to
the rate of re-matching. In our tests, we choose to re-match at a rate that
keeps the computational cost of matching minimal compared to the overall
training cost, such as every 10 or 50 V-Cycles.

Algorithm 2 Two-Level Full Approximation Scheme for Deep Neural Networks

$\mathbf{x},g$

while not converged do

$\mathbf{x}\leftarrow\text{smooth}(\mathbf{x},g)$

Construct ${\boldsymbol{\Pi}},{\mathbf{P}}$ via Heavy Edge Matching

$\mathbf{x}_{c}\leftarrow{\boldsymbol{\Pi}}\mathbf{x},\quad\mathbf{m}_{c}\leftarrow{\boldsymbol{\Pi}}\mathbf{m}$
$\triangleright$ (25)

$\mathbf{\tau}\leftarrow\frac{N}{m}\left(\left(\nabla_{\mathbf{x}_{c}}g^{c}_{\mathbf{s}}\right)({\boldsymbol{\Pi}}\mathbf{x})-\mathbf{R}\nabla_{\mathbf{x}}g_{\mathbf{s}}(\mathbf{x})\right)$
$\triangleright$ (35)

$\mathbf{x}_{c}\leftarrow\text{smooth}(\mathbf{x}_{c},{\widehat{g}_{c}}(\mathbf{x}_{c}))$,
where
${\widehat{g}_{c}}(\mathbf{x}_{c})=g_{c}(\mathbf{x}_{c})-\mathbf{x}_{c}^{T}\mathbf{\tau}$
$\triangleright$ (3)

$\mathbf{x}\leftarrow\mathbf{x}+{\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x}),\quad\mathbf{m}\leftarrow\mathbf{m}+{\mathbf{P}}(\mathbf{m}_{c}-{\boldsymbol{\Pi}}\mathbf{m})$
$\triangleright$ (4)

$\mathbf{x}\leftarrow\text{smooth}(\mathbf{x},g)$

end while

####  6.1.4 Maintaining Stability

In some cases, this multilevel hierarchies developed here can cause
instabilities, and in severe cases, divergence of the method. Here we discuss
a number of hyperparameters that can be used in the method to maintain
stability when this happens:

  * •

The learning rate can be decreased at coarser levels of the hierarchy. We
often choose to decrease the learning rate by a fixed factor $\eta$ at each
hierarchy level. This can prevent the coarse parameters from moving outside
the area in which the coarse problem approximates the fine problem
effectively. However, the learning rate does not impact the stored momentum
vector [11], so that must be stabilized separately.

  * •

The coarse-grid parameter correction and coarse-grid momentum correction in
Algorithm 2 and in Equations (4), (36) can be scaled down with parameters
$\alpha_{p},\alpha_{m}\in(0,1]$ as follows:

$\displaystyle\mathbf{x}\leftarrow$
$\displaystyle\mathbf{x}+\alpha_{p}{\mathbf{P}}(\mathbf{x}_{c}-{\boldsymbol{\Pi}}\mathbf{x})$
$\displaystyle\mathbf{m}\leftarrow$
$\displaystyle\mathbf{m}+\alpha_{m}{\mathbf{P}}(\mathbf{m}_{c}-{\boldsymbol{\Pi}}\mathbf{m}).$

This is akin to the $\alpha$ parameter used in the line search in MG/OPT [21];
however, we use constants here rather than the line search of MG/OPT due to
the added computational cost.

  * •

The $\tau$-correction in Algorithm 2 and in Equation (3) can be scaled down by
$\gamma\in[0,1]$, so that

${\widehat{g}_{c}}(\mathbf{x}_{c})=g_{c}(\mathbf{x}_{c})-\gamma\mathbf{x}_{c}^{T}\mathbf{\tau}.$

Because our $\tau$-correction is stochastic, even when we are near the
optimum, any given set of minibatches could potentially result in a much
larger $\mathbf{\tau}$, which could then cause an inappropriate coarse-grid
correction.

We note that in MG/OPT [21], it is also suggested one could add a quadratic
penalty to the coarse problem to prevent too much movements; however, MG/OPT
expects one to fully solve the coarse problem, whereas we only take a few SGD
steps. In our experiments, tuning the weight of that quadratic penalty caused
its own stability challenges, so we do not use that strategy.

When making these adjustments in Section 8, we typically choose
$\eta\in[\sqrt{2},2\sqrt{2}],\alpha_{p}=1,\alpha_{m}=0.2,\gamma=0.125$. These
were chosen by trial and error.

##  7 Extension to Convolutional Neural Networks

In this section we describe how the above method can be extended from fully-
connected networks to convolutional neural networks (CNNs). A CNN typically
consists of one or more convolutional layers followed by one or more fully-
connected layers. They are often effective in learning patterns over array-
structured data. In a CNN, the object analogous to a fully-connected neuron is
the _channel_ ; therefore, our goal in this section will be to develop methods
for coarsening over the channels of a CNN. The resulting convolution operators
will have the same width and stride, but the number of channels in each layer
will be restricted.

Convolutional operators are often represented as higher-order tensors, but
this tensor acts on the inputs in a linear fashion, and so the action can be
represented as matrix-vector multiplication. As a result, convolutional layers
are mathematically quite similar to fully-connected layers, and the above
method extends in a straightforward manner. That is, the action of a
convolutional layer can still be viewed in the form

$\mathbf{y}_{k+1}=\varphi(W_{k}\mathbf{y}_{k}+\mathbf{b}_{k}),$

where $W_{k}$ and $\mathbf{b}_{k}$ have a particular structure.

To apply our restriction and prolongation operators, we require an
understanding of the structure of the matrix $W$ associated with a
convolutional layer. This depends on the dimension of the input array, the
number of input channels, and the number of output channels. We build up our
understanding progressively as follows:

  * •

Consider input $\mathbf{y}$ to a convolutional layer that is a 1-dimensional
array and has a single channel; for example, audio data has this structure.
For simplicity, suppose we have a convolutional layer with a single
convolutional channel, which we represent as a vector
$\mathbf{c}=(c)_{i=-\ell}^{\ell}$ and a scalar $b$. Then the action of the
convolution is

$\sum_{j=-\ell}^{\ell}y(i+\ell-j)\cdot c(j)+b.$ (37)

We can represent this as a matrix-vector multiplication $W\mathbf{y}$ where
$W$ is the banded Toeplitz matrix

$W=\begin{bmatrix}c_{-\ell}&c_{-\ell+1}&\cdots&c_{k-1}&c_{k}&0&\cdots&0&0&0\\\
0&c_{-\ell}&c_{-\ell+1}&\cdots&c_{-1}&c_{k}&0&\cdots&0&0\\\
0&0&c_{-\ell}&c_{-\ell+1}&\cdots&c_{-1}&c_{k}&0&\cdots&0\\\
\vdots&&&&&&&&&\vdots\end{bmatrix}.$ (38)

  * •

More complex convolutional layers extend this in a natural way. One-
dimensional, single-channel arrays with multiple convolutional channels
results in a matrix

$W=[W^{1};\cdots;W^{p}]^{T}$

that is a column of banded Toeplitz blocks.

  * •

Two-dimensional, single-channel arrays with a single convolutional channel
produces a weight matrix $W$ that has banded block-Toeplitz structure. See
Figure 1.

Figure 1: The banded block-Toeplitz structure of a convolutional weight matrix
with a single channel applied to a 2D input with a single channel.

  * •

Adding additional output channels creates a block-column matrix in which each
block has this structure. Thus, two-dimensional, single-channel arrays with
multiple convolutional channels produce a column of matrix blocks each of the
form of Figure 1. We might say this then has _block banded block Toeplitz_
structure.

  * •

Finally, multi-dimensional input arrays with multiple inputs channels and
multiple output channels results in a weight matrix $W$ of the form

$W=\begin{bmatrix}W^{1,1}&W^{1,2}&\cdots&W^{1,\beta}\\\
\cdots&\vdots&&\cdots\\\
W_{\alpha,1}&W_{\alpha,2}&\cdots&W^{\alpha,\beta}\end{bmatrix}$ (39)

in which each submatrix $W^{i,j}$ is a matrix of the form of Figure 1.

A key observation here is that every block row corresponds to an output
channel and every block column corresponds to an input channel. As each
channel (block row) here is analogous to a neuron (individual row) in a fully-
connected layer, we see that, analogous to Section 5.4, we can coarsen two
channels together by simply averaging their channels together.

Therefore, we see that applying Algorithm 2 to CNNs proceeds in exactly the
same manner. Restrictions occurs through averaging channels as in $\pi_{k}$ of
Equation (32), and prolongation occurs through piecewise constant
interpolation of channels as in $P_{k}$ of Equation (32).

##  8 Experiments

In this section we demonstrate the use of our multilevel FAS trainer on a set
of PDE-based regression problems to compare this to SGD. These experiments
were performed on the Lassen compute cluster at Lawrence Livermore National
Laboratory. For each experiment, we used a single IBM Power9 nodes with 44 CPU
cores at 3.5 Ghz and 256 GB of memory and a single NVIDIA V100 Voltage
graphics card. What we find with these results is that in some cases the use
of FAS can improve $L_{2}$ loss, and in all cases it achieves $L_{2}$ loss at
least as good as SGD. In addition, in all cases tested, it _also_ improves
upon $L_{\infty}$ loss even though the training was against $L_{2}$. This
implies that FAS-trained neural networks has better worst-case performance
than those trained with traditional SGD, suggesting better generalization
power.

Note that in these experiments, to achieve a more faithful comparison in
computational work, we borrow a strategy from the multigrid literature and
report _work units_ rather than epochs. A single work unit, as we define it,
is the computational cost to compute a forward and backpropagation over one
minibatch. We do this because coarser levels of the multilevel hierarchy are
less expensive to compute over; as a rule of thumb, half as many neurons in
each level implies weight matrices with 1/4 as many elements, and so
computation at a coarse level requires roughly 1/4 as much work as the next-
finer level. For example, if we perform a V-Cycle over a 2-level hierarchy
with 2 minibatch presmoothing iterations, 2 minibatch coarse grid iterations,
and then 2 minibatch postsmoothing iterations, that costs 4.5 work units. In
addition, computing the tau-correction requires both a fine-level and a
coarse-level forward and backpropagation over each minibatch that will be used
at the coarse level; in our example, that would add an additional 2.5 work
units, for a total of 7 work units for the V-Cycle. In general, if $s$ is the
cost of a forward/backward pass, then work units can be computed as

$\displaystyle\mbox{(smoothing cost)}+\mbox{(tau correction work for current
level)}+\mbox{(tau correction work for next level)}$ $\displaystyle=$
$\displaystyle\left(\sum_{i=0}^{k-2}(1/4)^{i}\cdot 2s+(1/4)^{k-1}\cdot
s\right)+\left(\sum_{i=1}^{k-2}(1/4)^{i}\cdot 2s+(1/4)^{k-1}\cdot
s\right)+\left(\sum_{i=0}^{k-3}(1/4)^{i}\cdot 2s+(1/4)^{k-1}\cdot s\right)$
$\displaystyle=$ $\displaystyle 6s-10s/4^{k-1}.$ (40)

We see that, for any hierarchy depth, Work Units is bounded from above by
$6s$.

As discussed in Section 2.2.1, deep learning, despite being formalized in the
same way as a traditional optimization problem, has different goals from
traditional optimization. Most significantly, the actual objective function of
interest, the neural network’s performance over the entire population of
_possible_ data, is inaccessible, and the training performance is only an
approximation to this. For this reason, the goal is not to achieve asymptotic
convergence to a fixed minimum point, but rather to achieve good performance
on a validation set.

###  8.1 Poisson Problem

We consider solving

$-\nabla\cdot(\kappa\nabla u)=f\quad\mbox{in}\quad\Omega=[0,1]^{2}$ (41)

on unit square $\Omega=[0,1]^{2}$ with a $32\times 32$ mesh and $u=0$ on the
boundary, where

$\displaystyle\kappa(x,y)$
$\displaystyle=1.1+\cos[k_{x}\pi(x^{\prime}+a_{x})]\cos[k_{y}\pi(y^{\prime}+a_{y})]$
$\displaystyle x^{\prime}$
$\displaystyle=\cos\alpha(x-0.5)-\sin\alpha(y-0.5)+0.5$ $\displaystyle
y^{\prime}$ $\displaystyle=\sin\alpha(x-0.5)+\cos\alpha(y-0.5)+0.5$
$\displaystyle f$ $\displaystyle=32e^{-4[(x-0.25)^{2}+(y-0.25)^{2}]}$

The data set has $10,000$ samples: frequencies $k_{x},k_{y}\in(0.5,4)$, phase
shifts $a_{x},a_{y}\in(0,0.5)$, and rotation angle $\alpha\in(0,\pi/2)$,
generated using FEniCS [2] and DOLFIN [17]. The goal is to train NNs to learn
the solution $u$ from $\kappa$, $f$, and the mesh point coordinates $(x,y)$.
So, the dimension of the input is $3\times 32\times 32$ and the dimension of
the output is $32\times 32$. Figure 2 shows an example of an exact solution
$u$ and the corresponding solution learned from NN.

(a) $\kappa$

(b) exact $u$

(c) learned $u$

Figure 2: Diffusion coefficient and the solution

We tested two neural network architectures:

  * •

A fully-connected network with two hidden layers of 400 neurons each, with
ReLU activation functions.

  * •

A convolutional network with two convolutional layers of 200 channels each,
kernel width 7 and stride 2, following by two fully-connected layers of 400
neurons each, with ReLU activation functions.

(a) training $L_{2}$ (b) training $L_{\infty}$

(c) validation $L_{2}$ (d) validation $L_{\infty}$

Figure 3: (a) Smoothed training and validation loss ($L_{2}$ and $L_{\infty}$)
for one-level SGD training and multilevel FAS trainings (with 2-4 levels) for
the Poisson problem using fully-connected NNs.

(a) training $L_{2}$ (b) training $L_{\infty}$

(c) validation $L_{2}$ (d) validation $L_{\infty}$

Figure 4: (a) Smoothed training and validation loss ($L_{2}$ and $L_{\infty}$)
for one-level SGD training and multilevel FAS trainings (with 2-4 levels) for
the Poisson problem using CNNs.

We trained these neural networks using traditional (one-level) SGD and
multilevel FAS with SGD as the smoother (from 2 levels and up to 4 levels). In
these cases, we used a learning rate of 0.01, momentum of 0.9, weight decay of
0, and minibatches of size 200. For the multilevel FAS, we gradually decrease
the learning rates for lower levels and recompute the neuron matching for
constructing the $P$ and $\Pi$ operators every 50 V-cycles. We perform 4 SGD
minibatch iterations at each smoothing step. The results of validation and
training losses are shown in Figure 3 for the fully-connected network and
Figure 4 for the convolutional network. We also show the best losses achieved
for each network with each method in Table 1 up to the first $10,000$ work
units.

|  | Validation Loss | Training Loss
---|---|---|---
|  | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$
Fully-Connected | One-Level SGD | 2.94e-3 | 0.626 | 2.84e-3 | 0.814
Two-Level FAS | 1.61e-3 | 0.496 | 1.59e-3 | 0.664
Convolutional | One-Level SGD | 5.79e-3 | 0.772 | 5.81e-3 | 1.022
Two-Level FAS | 5.38e-3 | 0.756 | 5.37e-3 | 0.991
Table 1: Best $L_{2}$ and $L_{\infty}$ losses achieved for fully-connected and
convolutional networks using both SGD and two-level FAS. We can see that the
the fully-connected network outperforms convolutional networks, and the two-
level FAS consistently outperforms SGD.

For the fully-connected NNs, we see that the multilevel FAS training methods
have significantly faster convergence rates with respect to work units than
the 1-level SGD in terms of both the $L_{2}$ and $L_{\infty}$ errors. With the
3- and 4-level methods, the speed of the convergence can be further improved
slightly over the 2-level method. Similar behavior can also be observed with
CNNs, where multilvel FAS methods yield superior convergence than SGD. We note
here that the 4-level method appears to struggle with stability at the
beginning of the training process. So, in this experiment, we did not further
increase the number of levels. Improving the stability of the proposed
approach with higher number of levels is left for our future work.

###  8.2 Darcy Flow Problem

Figure 5: A spatially-correlated permeability field used as input to a Darcy
problem.

The test data correspond to the finite element solution of the Darcy equation.
More specifically, for a single phase porous media flow, we compute the flux
$\mathbf{q}$ and the pressure $p$ by solving the following mixed system posed
on the unit square $D$: Given two (Raviart-Thomas) finite element spaces,
$\mathbf{V}_{h}\subset\mathbf{H}(\operatorname{div},\;D)$ and $W_{h}\subset
L_{2}(D)$, find $\mathbf{q}\in\mathbf{V}_{h}$ and $p\in W_{h}$, such that

$\begin{array}[]{lll}(k^{-1}\mathbf{q},\;\mathbf{s})-(\operatorname{div}\mathbf{s},\;p)&=(\mathbf{f},\;\mathbf{s}),&\text{
for all }\mathbf{s}\in{\mathbf{V}}_{h},\\\
(\operatorname{div}\mathbf{q},\;w)&=0,&\text{ for all }w\in W_{h},\end{array}$
(42)

subject to Dirichlet boundary conditions $p=p_{D}$ on $\Gamma_{D}$ enforced by
the right-hand side $\mathbf{f}$, and Neumann boundary condition
$\mathbf{q}\cdot\mathbf{n}=0$ on $\Gamma_{N}$, where $\Gamma_{D}$ and
$\Gamma_{N}$ are nonintersecting and cover the boundary $\partial D$.

For a given mesh ${\mathcal{T}}_{h}$, we have the following data:

  * •

The input is $k=\exp(u)$, where $u$ is a realization of log-normal
permeability field with piecewise constant values over the elements of
${\mathcal{T}}_{h}$.

  * •

The output data is either the solution $p\in W_{h}$, also a piecewise constant
quantity over the elements of ${\mathcal{T}}_{h}$, or a specific quantity of
interest, QoI. In our case, we choose as QoI,

$Q=\frac{1}{|\Gamma_{out}|}\;\int\limits_{\Gamma_{out}}\mathbf{q}\cdot\mathbf{n}\;dS,$
(43)

representing the average flux across the outflow boundary $\Gamma_{out}$,
where $\mathbf{n}$ is outer unit vector normal to $\Gamma_{out}\subset\partial
D$ [9].

This type of data $(u,\;Q)$ is used in Monte Carlo simulations to estimate the
expected value of the quantity of interest. See Figure 5 for an example of a
permeability field $u$.

(a)  (b)

Figure 6: (a) Smoothed $L_{2}$ loss and (b) smoothed $L_{\infty}$ loss for
one-level training (i.e. SGD) and two-level FAS for the Darcy problem using a
fully-connected network.

####  8.2.1 2-Level Without Tuning

(a)  (b)

Figure 7: (a) Smoothed $L_{2}$ loss and (b) smoothed $L_{\infty}$ loss for
one-level training (i.e. SGD) and two-level FAS for the Darcy problem using a
CNN.

We show the validation losses achieved for two types of networks:

  * •

A fully-connected network with 3 layers, each with 1024 neurons.

  * •

A convolutional neural network with 3 convolutional layers. The channels,
kernel width, and stride of these layers are (100, 200, 400), (11, 7, 3), and
(2, 1, 1) respectively. Following these are two fully-connected layers each
with 2048 neurons.

We trained these two networks using both traditional (one-level) SGD and our
two-level FAS. In these cases we used a learning rate of 0.01, momentum of
0.9, weight decay of 1e-6, with minibatches of size 200. For two-level FAS, we
recompute matchings with which to construct the $P$ and $\Pi$ operators every
10 V-cycles, and perform 2 SGD minibatch iterations at each smoothing step.

|  | Validation Loss | Training Loss
---|---|---|---
|  | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$
Fully-Connected | One-Level SGD | 1.49e-3 | 0.707 | 2.37e-5 | 0.0546
Two-Level FAS | 1.39e-3 | 0.472 | 4.75e-5 | 0.0540
Convolutional | One-Level SGD | 6.61e-4 | 0.361 | 2.92e-5 | 0.0469
Two-Level FAS | 4.84e-4 | 0.312 | 7.47e-5 | 0.0696
Table 2: Best $L_{2}$ and $L_{\infty}$ losses achieved for fully-connected and
convolutional networks using both SGD and two-level FAS. We can see that the
CNN consistenly outperforms the fully-connected network, and the two-level FAS
consistenly outperforms SGD.

The results can be seen in Figure 6 for the fully-connected network and Figure
7 for the convolutional network using both traditional SGD and our two-level
FAS. As stochastic training methods tend to produce high performance variation
from step to step, raw plots of performance over many steps can be messy; to
improve the visualization these plots are smoothed using a moving arithmetic
mean with a window size of 33, chosen as the smallest window size we felt
produced informative visuals. We also show the best losses achieved for each
network with each method in Table 2. We can see from these results that, as
with the Poisson problem, two-level FAS produces significant improvements over
SGD.

###  8.3 Deeper Hierarchies with Stability Tuning

Here we investigate the use of deeper hierarchies. We find that the 3-level
hierarchy tends to suffer from stability problems, while the 4-level and
5-level hierarchies are much more stable. Understanding this stability
phenomena is the subject of future work.

In order to maintain stability, we make the following hyperparameter choices:

  * •

The momentum coarse-grid correction scaling parameter is set to
$\alpha_{m}=0.2$.

  * •

The tau correction scaling parameter is set to $\gamma=0.125$.

  * •

The learning rate on each hierarchy level is divided by a factor of
$\eta=2\sqrt{2}$ for the first 3 levels, after which the learning rate is kept
the same as the third level.

|
---|---
|
|
|

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 8: $L_{2}$ loss (left) and $L_{\infty}$ validation loss (right) for
fully-connected networks trained with 2, 3, 4, and 5 hierarchy levels. Each
plot shows 1-level loss (the same in each plot), multilevel loss, and first
auxiliary loss.

The results can be seen in Figure 8, where we show validation error from both
the original neural network and its first-coarsened neural network. In each
plot in this Figure, we also show the 1-level SGD results as a point of common
comparison. We find that not only does the fine neural network perform well,
but that the first auxiliary network (the first-coarsened network) often has
strong performance as well. In fact, what one can observe in Figure 8 is that
the original neural network typically performs best in the L2 error, while the
first auxiliary network performs best in the $L_{\infty}$ error. We believe
this is because the first-coarsened network has roughly half as many neurons
and therefore has less opportunity for overfitting. We have found that this
strong auxiliary performance only occurs when the tau scaling parameter
$\gamma$ is sufficiently small. We also find that the 3-level hierarchy
appears unusually unstable; this is not fully understood, and deserves further
investigation.

We also show the best losses seen, for both validation and training loss
across different hierarchy depths with both fully-connected and convolutional
networks, in Table 3. These results demonstrate the regularization impact of
multilevel training. When dealing with training loss, either the one-level
method or a shallower hierarchy produces the best losses. But when considering
validation loss, which is the stronger metric of generalization performance of
the network, the deeper hierarchies produce significantly stronger
performance. This is especially true when considering $L_{\infty}$ loss; this
is consistent with the notion of regularization as a mechanism to prevent
inappropriate corner-case behavior.

| Full-Connected Network | Convolutional Network
---|---|---
| Validation Loss | Training Loss | Validation Loss | Training Loss
Levels | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$
1 | 1.49e-3 | 0.677 | 2.20e-5 | 0.0544 | 6.61e-4 | 0.361 | 2.92e-5 | 0.0469
2 | 1.21e-3 | 0.434 | 4.66e-6 | 0.0165 | 6.06e-4 | 0.371 | 5.38e-5 | 0.0867
2aux | 1.35e-3 | 0.492 | 4.25e-4 | 0.107 | 6.32e-4 | 0.392 | 2.53e-4 | 0.0954
3 | 2.15e-3 | 0.573 | 1.31e-3 | 0.747 | 1.05e-3 | 0.361 | 6.84e-4 | 0.686
3aux | 3.16e-3 | 0.642 | 2.28e-3 | 0.656 | 1.34e-3 | 0.512 | 1.03e-3 | 0.395
4 | 1.20e-3 | 0.484 | 6.71e-6 | 0.0195 | 6.65e-4 | 0.293 | 6.69e-5 | 0.0768
4aux | 1.26e-3 | 0.345 | 4.24e-4 | 0.110 | 6.13e-4 | 0.238 | 2.56e-4 | 0.0991
5 | 1.42e-3 | 0.456 | 1.16e-5 | 0.0303 | 7.28e-4 | 0.378 | 8.58e-5 | 0.101
5aux | 1.46e-3 | 0.336 | 5.05e-4 | 0.124 | 6.03e-4 | 0.398 | 2.40e-4 | 0.0818
Table 3: Best $L_{2}$ and $L_{\infty}$ losses achieved for fully-connected and
convolutional networks using both SGD and FAS with stability control. Level
$k$aux refers to the first auxiliary network in the $k$-level hierarchy.

###  8.4 Experiments with Learning Rates

|
---|---

(a)  (b)

Figure 9: $L_{2}$ validation loss (left) and $L_{\infty}$ validation loss
(right) for varying hierarchical learning $\eta$ parameters. A smoothing
window of size 65 was used

There are many possible hyperparameter and architecture choices one could
make, and we cannot present a full exploration of all hyperparameter choices
here. We instead present experiments in which we vary how the learning rate
changes with hierarchy for the 4-level hierarchy. In the previous section we
divided the learning rate by a factor of $2\sqrt{2}$ as the hierarchy deepens
for the first three levels; here we also consider division constants of 1,
$\sqrt{2}$, 2, and 4. These results can be see in Figure 9. Best losses seen
across the training runs, of both validation and training losses, can be see
in Table 4.

| Full-Connected Network
---|---
| Validation Loss | Training Loss
$\eta$ | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$
1 | 1.21e-3 | 0.412 | 1.57e-5 | 0.0328
1, Aux1 | 1.40e-3 | 0.347 | 3.00e-4 | 0.0871
$\sqrt{2}$ | 1.08e-3 | 0.364 | 9.01e-6 | 0.0218
$\sqrt{2}$, Aux1 | 1.26e-3 | 0.297 | 3.40e-4 | 0.0921
2 | 1.09e-3 | 0.602 | 7.10e-6 | 0.0196
2, Aux1 | 1.24e-3 | 0.345 | 3.83e-4 | 0.106
$2\sqrt{2}$ | 1.20e-3 | 0.484 | 6.71e-6 | 0.0195
$2\sqrt{2}$, Aux1 | 1.26e-3 | 0.345 | 4.24e-4 | 0.110
4 | 1.29e-3 | 0.396 | 1.21e-5 | 0.0260
4, Aux1 | 1.31e-3 | 0.345 | 4.45e-4 | 0.116
Table 4: Best $L_{2}$ and $L_{\infty}$ losses achieved for fully-connected
networks trained using 4-level FAS with varying learning rate decay $\eta$
parameters. Learning rate decays for the first 3 levels of the hierarchy, and
then remains the same for the fourth.

Like previous experiments, we see from these results that the first auxiliary
network tends to produce stronger $L_{\infty}$ performance, which we
hypothesize to be due to the reduced degrees of freedom having less
opportunity for overfitting. Of these options, $\eta=\sqrt{2}$ seems to
perform best: its $L_{2}$ loss initially drops as quickly as $\eta=1$ but
eventually outperforms it, and its auxiliary network achieves the best
$L_{\infty}$ loss.

##  9 Conclusion

In this paper we have presented a multilevel-in-width method for training
neural networks by developing a variation of the Full Approximation Scheme as
originally used in Algebraic Multigrid. We have shown that DNNs trained in
this way for regression problems can achieve better generalization power by
improving on multiple validation performance metrics. We also observe that, in
some cases, the first auxiliary network achieves better worst-case performance
than the original network, likely due to having less opportunity for
overfitting.

This paper is perhaps most interesting in the context of [27] and [7], which
also showed regularization benefits to multilevel training even though all
have taken quite different approaches. Future investigations should include
work that seeks to more-fully understand the settings in which multilevel
training confers regularization benefits, and to develop principled approaches
for making specific algorithmic choices in those settings.

Our open-source PyTorch-based software, MTNN, to implement these methods, is
written in a modular fashion to facilitate the study of different algorithmic
choices and is available on Github at https://github.com/LLNL/MTNN. We have
also provided open-source data and drivers that allow users to easily recreate
our experiments. Readers are encouraged to download and use our software.

## Acknowledgments

We would like to thank Tom Benson and Ulrike Yang for their insight and
suggestions throughout the research and writing of this paper.

## References

  * [1] et. al. Adam Paszke.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [2] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells.  Archive of numerical software: The FEniCS project version 1.5.  Archive of Numerical Software, 3(100), 2015.
  * [3] Mordecai Avriel.  Nonlinear Programming: Analysis and Methods.  Dover Publishing, 2003.
  * [4] Achi Brandt.  Multi-level adaptive solutions to boundary-value problems.  Mathematics of computation, 31(138):333–390, 1977.
  * [5] William L. Briggs, Van Emden Hensen, and Steve F. McCormick.  A Multigrid Tutorial, 2nd Ed. Society for Industrial and Applied Mathematics, 2000.
  * [6] H. Calandra, Serge Gratton, Elisa Riccietti, and Xavier Vasseur.  On a multilevel levenberg-marquardt method for the training of artificial neural networks and its application to the solution of partial differential equations.  Optimization Methods and Software, pages 1–26, 6 2020.
  * [7] Eric C. Cyr, Stefanie Günther, and Jacob B. Schroder.  Multilevel initialization for layer-parallel deep neural network training, 2019.
  * [8] J. Duchi, E. Hazan, and Y. Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12:2121–2159, 2011.
  * [9] Hillary R. Fairbanks, Umberto Villa, and Panayot S. Vassilevski.  Multilevel hierarchical decomposition of finite element white noise with application to multilevel Markov chain Monte Carlo.  SIAM Journal on Scientific Computing, pages S293–S316, 2020.
  * [10] Lisa Gaedke-Merzhäuser, Alena Kopaničáková, and Rolf Krause.  Multilevel minimzation for deep residual networks.  In ESAIM: Proceedings and Surveys, volume 71, pages 131–144, 2021\.
  * [11] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [12] Stefanie Günther, Lars Ruthotto, Jacob Schroder, Eric Cyr, and Nicolas Gauger.  Layer-parallel training of deep residual neural networks.  SIAM Journal on the Mathematics of Data Science, 2(1):1–23, 2020\.
  * [13] Wolfgang Hackbusch.  Multi-grid methods and applications, volume 4.  Springer Science & Business Media, 2013.
  * [14] George Karypis and Vipin Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM Journal on Scientific Computing, 20(1):359–392, 1998.
  * [15] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [16] Andrew Kirby, Siddharth Samsi, Michael Jones, Albert Reuther, Jeremy Kepner, and Vijay Gadepally.  Layer-parallel training with gpu concurrency of deep residual neural networks via nonlinear multigrid.  In 2020 IEEE High Performance Extreme Computing Conference (HPEC), pages 1–7, 2020.
  * [17] Anders Logg and Garth N. Wells.  Dolfin: Automated finite element computing.  ACM Trans. Math. Softw., 37(2), April 2010.
  * [18] James Lucas, Shengyang Sun, Richard Zemel, and Roger Grosse.  Aggregated momentum: Stability through passive damping.  In International Conference on Learning Representations, 2019.
  * [19] Jerry Ma and Denis Yarats.  Quasi-hyperbolic momentum and Adam for deep learning.  In International Conference on Learning Representations, 2019.
  * [20] et. al. Martín Abadi.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [21] Stephen G. Nash.  A multigrid approach to discretized optimization problems.  Optimization Methods and Software, 14(1-2):99–116, 2000.
  * [22] M. Naumov, M. Arsaev, P. Castonguay, J. Cohen, J. Demouth, J. Eaton, S. Layton, N. Markovskiy, I. Reguly, N. Sakharnykh, V. Sellappan, and R. Strzodka.  Amgx: A library for gpu accelerated algebraic multigrid and preconditioned iterative methods.  SIAM Journal on Scientific Computing, 37(5):S602–S626, 2015.
  * [23] Yvan Notay.  Aggregation-based algebraic multigrid for convection-diffusion equations.  SIAM Journal on Scientific Computing, 34(4):A2288–A2316, 2012.
  * [24] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, pages 234–241, Cham, 2015. Springer International Publishing.
  * [25] Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton.  On the importance of initialization and momentum in deep learning.  In Proceedings of the 30th International Conference on Machine Learning, page III–1139–III–1147. JMLR.org, 2013.
  * [26] Hong Hui Tan and King Hann Lim.  Review of second-order optimization techniques in artificial neural networks backpropagation.  IOP Conference Series: Materials Science and Engineering, 495:012003, June 2019.
  * [27] Cyrill von Planta, Alena Kopanicáková, and Rolf Krause.  Training of deep residual networks with stochastic MG/OPT.  CoRR, abs/2108.04052, 2021.
  * [28] Pieter Wesseling.  Introduction to multigrid methods.  Technical report, 1995.