#  Bilevel Methods for Image Reconstruction

Crockett  Caroline  Fessler  Jeffrey A  Department of EECS, University of
Michigan, Ann Arbor, Michigan, USA<EMAIL_ADDRESS>
###### Abstract

This review discusses methods for learning parameters for image reconstruction
problems using bilevel formulations. Image reconstruction typically involves
optimizing a cost function to recover a vector of unknown variables that
agrees with collected measurements and prior assumptions. State-of-the-art
image reconstruction methods learn these prior assumptions from training data
using various machine learning techniques, such as bilevel methods.

One can view the bilevel problem as formalizing hyperparameter optimization,
as bridging machine learning and cost function based optimization methods, or
as a method to learn variables best suited to a specific task. More formally,
bilevel problems attempt to minimize an upper-level loss function, where
variables in the upper-level loss function are themselves minimizers of a
lower-level cost function.

This review contains a running example problem of learning tuning parameters
and the coefficients for sparsifying filters used in a regularizer. Such
filters generalize the popular total variation regularization method, and
learned filters are closely related to convolutional neural networks
approaches that are rapidly gaining in popularity. Here, the lower-level
problem is to reconstruct an image using a regularizer with learned
sparsifying filters; the corresponding upper-level optimization problem
involves a measure of reconstructed image quality based on training data.

This review discusses multiple perspectives to motivate the use of bilevel
methods and to make them more easily accessible to different audiences. We
then turn to ways to optimize the bilevel problem, providing pros and cons of
the variety of proposed approaches. Finally we overview bilevel applications
in image reconstruction.

Caroline Crockett
University of Michigan
<EMAIL_ADDRESS>and Jeffrey A. Fessler
University of Michigan
<EMAIL_ADDRESS>copyrightowner=A. Heezemans and M. Casey, volume = xx, issue
= xx, pubyear = 2021, isbn = xxx-x-xxxxx-xxx-x, eisbn = xxx-x-xxxxx-xxx-x, doi
= 10.1561/XXXXXXXXX, firstpage = 1, lastpage = 18 bilevel,litreview.bib
bilevel,litreview,jeff.bib

###### Contents

  1. 1 Introduction
    1. 1.1 Notation
    2. 1.2 Bilevel Method: Set-up and Running Example
    3. 1.3 Conclusion
  2. 2 Background: Image Reconstruction
    1. 2.1 Image Reconstruction
    2. 2.2 Sparsity-Based Regularizers
      1. 2.2.1 Synthesis Regularizers
      2. 2.2.2 Analysis Regularizers
      3. 2.2.3 Comparing Analysis and Synthesis Approaches
    3. 2.3 Brief History of Analysis Regularizer Learning
    4. 2.4 Summary
  3. 3 Background: Hyperparameter Optimization
    1. 3.1 Image Quality Metrics
      1. 3.1.1 Full Reference IQA
      2. 3.1.2 No Reference IQA
    2. 3.2 Parameter Search Strategies
      1. 3.2.1 Model-free
      2. 3.2.2 Model-based
    3. 3.3 Summary
  4. 4 Gradient Based Bilevel Methodology: The Groundwork
    1. 4.1 Set-up
    2. 4.2 Minimizer Approach
      1. 4.2.1 Implicit Function Theorem Perspective
      2. 4.2.2 KKT Conditions
      3. 4.2.3 Summary of Minimizer Approach
      4. 4.2.4 Computational Costs
    3. 4.3 Translation to a Single Level
    4. 4.4 Unrolled Approaches
      1. 4.4.1 Computational and Memory Complexity
    5. 4.5 Summary
  5. 5 Gradient-Based Bilevel Optimization Methods
    1. 5.1 Double-Loop Algorithms
      1. 5.1.1 Complexity Analysis
    2. 5.2 Single-Loop Algorithms
      1. 5.2.1 Complexity Analysis
    3. 5.3 Summary of Methods
  6. 6 Survey of Applications
    1. 6.1 Lower-level Cost Function Design
    2. 6.2 Upper-Level Loss Function Design
    3. 6.3 Conclusion
  7. 7 Connections and Future Directions
    1. 7.1 Connection: Learnable Optimization Algorithms
    2. 7.2 Connection: Plug-and-play Priors
    3. 7.3 Connection: Single-Level Parameter Learning
    4. 7.4 Future Directions
    5. 7.5 Summary of Advantages and Disadvantages
  8. A Additional Running Example Results
    1. A.1 Derivatives for Convolutional Filters
    2. A.2 Evaluating Conditions for the Running Example
      1. A.2.1 Upper-level Loss Conditions
      2. A.2.2 Lower-level Cost Conditions

##  Chapter 1 Introduction

Methods for image recovery aim to estimate a good-quality image from noisy,
incomplete, or indirect measurements. Such methods are also known as
computational imaging. For example, image denoising and image deconvolution
attempt to recover a clean image from a noisy and/or blurry input image, and
image inpainting tries to complete missing measurements from an image. Medical
image reconstruction aims to recover images that humans can interpret from the
indirect measurements recorded by a system like an MRI or CT scanner. Such
image reconstruction applications are a type of inverse problem [engl:96].

New methods for image reconstruction attempt to lower complexity, decrease
data requirements, or improve image quality for a given input data quality.
For example, in CT one goal is to provide doctors with information to help
their patients while reducing radiation exposure [mccollough:17:ldc]. To
achieve these lower radiation doses, the CT system must collect data with
lower beam intensity or fewer views. Similarly, in MRI collecting fewer
k-space samples can reduce scan times. Such “undersampling” leads to an under-
determined problem, with fewer knowns (measurements from a scanner) than
unknowns (pixels in the reconstructed image), requiring advanced image
reconstruction methods.

Existing reconstruction methods make different assumptions about the
characteristics of the images being recovered. Historically, the assumptions
are based on easily observed (or assumed) characteristics of the desired
output image, such as a tendency to have smooth regions with few edges or to
have some form of sparsity [eldar:12:cs]. More recent machine learning
approaches use training data to discover image characteristics. These
learning-based methods often outperform traditional methods, and are gaining
popularity in part because of increased availability of training data and
computational resources [wang:16:apo, hammernik:2020:machinelearningimage].

There are many design decisions in learning-based reconstruction methods. How
many parameters should be learned? What makes a set of parameters “good?” How
can one learn these good parameters? Using a bilevel methodology is one
systematic way to address these questions.

Bilevel methods are so named because they involve two “levels” of
optimization: an upper-level loss function that defines a goal or measure of
goodness (equivalently, badness) for the learnable parameters and a lower-
level cost function that uses the learnable parameters, typically as part of a
regularizer. The main benefits of bilevel methods are learning task-based
hyperparameters in a principled approach and connecting machine learning
techniques with image reconstruction methods that are defined in terms of
optimizing a cost function, often called model-based image reconstruction
methods. Conversely, the main challenge with bilevel methods is the
computational complexity. However, like with neural networks, that complexity
is highest during the training process, whereas deployment has lower
complexity because it requires only the lower-level problem.

The methods in this review are broadly applicable to bilevel problems, but we
focus on formulations and applications where the lower-level problem is an
image reconstruction cost function that uses regularization based on analysis
sparsity. The application of bilevel methods to image reconstruction problems
is relatively new [hammernik:2020:machinelearningimage], but there are a
growing number of promising research efforts in this direction. We hope this
review serves as a primer and unifying treatment for readers who may already
be familiar with image reconstruction problems and traditional regularization
approaches but who have not yet delved into bilevel methods.

This review lies at the intersection of a specific machine learning method,
bilevel, and a specific application, filter learning for image reconstruction.
For a broad overview of machine learning in image reconstruction, see
[hammernik:2020:machinelearningimage]. For an overview of image reconstruction
methods, including classical, variational, and learning-based methods, see
[mccann:2019:biomedicalimagereconstruction]. Finally, for a historical
overview of bilevel optimization and perspective on its use in a wide variety
of fields, see [dempe:2020:bileveloptimizationadvances]. Within the image
recovery field, bilevel methods have also been used, e.g., in learning
synthesis dictionaries [mairal:2012:taskdrivendictionarylearning,
chambolle:2021:learningconsistentdiscretizations].

The structure of this review is as follows. The remainder of the introduction
defines our notation and presents a running example bilevel problem. Section 2
provides background information on the lower-level image reconstruction cost
function and analysis regularizers. Section 3 provides background information
on the upper-level loss function, specifically loss function design and
hyperparameter optimization strategies. These background sections are not
meant to be exhaustive overviews of these broad topics. Section 4 presents
building blocks for optimizing a bilevel problem. Section 5 uses these
building blocks to discuss optimization methods for the upper-level loss
function. Section 6 discusses previous applications of the bilevel method in
image recovery problems, including signal denoising, image inpainting, and
medical image reconstruction. It also overviews bilevel formulations for blind
learning and learning space-varying tuning parameters. Finally, Section 7
offers summarizing commentary on the benefits and drawbacks of bilevel methods
for computational imaging and our thoughts for future directions for the
field.

###  1.1 Notation

This review considers continuous-valued, discrete space signals. Some papers,
e.g., [calatroni:2017:bilevelapproacheslearning,
delosreyes:2017:bilevelparameterlearning], analyze signals in function space,
arguing that the goal of high resolution imagery is to approximate a
continuous space reality and that analysis in continuous domain can yield
insights and optimization algorithms that are resolution independent. However,
the majority of bilevel methods are motivated and described in discrete space.
The review does not include discrete-valued settings, such as image
segmentation [ochs:2015:bileveloptimizationnonsmooth]; those problems often
require different techniques to optimize the lower-level cost function.

The literature is inconsistent in how it refers to variables in machine
learning problems. For consistency within this document, we define the
following terms:

  * •

Hyperparameters: Any adjustable parameters that are part of a model. Tuning
parameters and model parameters are both sub-types of hyperparameters. This
document uses $\bm{\upgamma}$ to denote a vector of hyperparameters.

  * •

Tuning parameters: Scalar parameters that weight terms in a cost function to
determine the relative importance of each term. This review uses $\beta$ to
denote individual tuning parameters.

  * •

Model parameters: Parameters, generally in vector or matrix form, that are
used in the structure of an cost or loss function, typically as part of the
regularization term. In the running example in the next section, the model
parameters are typically filter coefficients, denoted $\bm{c}$.

We write vectors as column vectors and use bold to denote matrices (uppercase
letters) and vectors (lowercase letters). Subscripts index vector elements, so
$x_{i}$ is the $i$th element in $\bm{x}$. For functions that are applied
element-wise to vectors, we use notation following the Julia programming
language [bezanson:17:jaf], where $f.({\bm{x}})$ denotes the function $f$
applied element wise to its argument:

${\bm{x}}\in\mathbb{F}^{N}\implies f.({\bm{x}})=\begin{bmatrix}f(x_{1})\\\
\vdots\\\ f(x_{N})\end{bmatrix}\in\mathbb{F}^{N}.$

The field $\mathbb{F}$ can be either $\mathbb{R}$ or $\mathbb{C}$, depending
on the application.

Variable |  Dim |  Description
---|---|---
${\bm{x}}^{\text{true}}_{j}$ |  $N$ |  One of $J$ clean, noiseless training signals. Often used in a supervised training set-up.
$\bm{A}$ |  $M\times N$ |  Forward operator for the system of interest.
${\bm{y}}_{j}$ |  $M$ |  During the bilevel learning process, this refers to simulated measurements, where ${\bm{y}}_{j}={\bm{A}}{\bm{x}}^{\text{true}}_{j}+{\bm{n}}_{j}$. Once $\bm{\upgamma}$ is learned, $\bm{y}$ refers to collected measurements.
$\bm{n}$ |  $N$ |  A noise realization.
${\bm{\hat{x}}}$ |  $N$ |  A reconstructed image.
$\bm{\upgamma}$ |  $R$ |  The vector of parameters to learn using bilevel methods. This often includes ${\bm{c}}_{k}$ and/or $\beta_{k}$.
${\bm{c}}_{k}$ |  $S$ |  One of $K$ convolutional filters. A 2D filter might be $\sqrt{S}\times\sqrt{S}$.
$\tilde{{\bm{c}}}_{k}$ |  $S$ |  Conjugate mirror reversal of filter ${\bm{c}}_{k}$
${\bm{C}}_{k}$ |  $N\times N$ |  The convolution matrix such that ${\bm{C}}_{k}{\bm{x}}={\bm{c}}_{k}\circledast{\bm{x}}$ and ${\bm{C}}_{k}^{\prime}{\bm{x}}=\tilde{{\bm{c}}}_{k}\circledast{\bm{x}}$.
$\beta_{k}$ |  $\mathbb{R}$ |  The tuning parameter associated with ${\bm{c}}_{k}$.
$\beta_{0}$ |  $\mathbb{R}$ |  An overall regularization (tuning) parameter, appearing as $e^{\beta_{0}}$ the running example.
$\bm{z}$ |  Varies |  A sparse vector, often from ${\bm{C}}_{k}{\bm{x}}$.
$\epsilon$ |  $\mathbb{R}_{+}$ |  Parameter used to define $\phi$. Typically determines the amount of corner-rounding.
$t$ |  $0,\ldots,T$ |  Iteration counter for lower-level optimization iterates, e.g., ${\bm{x}}^{(t)}$ is the estimate of the lower-level optimization variable $\bm{x}$ at the $t$th iteration.
$u$ |  $0,1,\ldots$ |  Iteration counter for the upper-level optimization iterates, e.g., ${\bm{\upgamma}}^{(u)}$.
Table 1.1:  Overview of frequently used symbols in the review.

Convolution between a vector, $\bm{x}$, and a filter, $\bm{c}$, is denoted as
${\bm{c}}\circledast{\bm{x}}$. This review assumes all convolutions use
circular boundary conditions. Thus, the same convolution is equivalent to
multiplication with a square, circulant matrix:

${\bm{c}}\circledast{\bm{x}}={\bm{C}}{\bm{x}}.$

The conjugate mirror reversal of $\bm{c}$ is denoted as $\tilde{{\bm{c}}}$ and
its application is equivalent to multiplying with the adjoint of $\bm{C}$:

$\tilde{{\bm{c}}}\circledast{\bm{x}}={\bm{C}}^{\prime}{\bm{x}}.$

Finally, for partial derivatives, we use the notation that

$\displaystyle\nabla_{\bm{x}}f({\bm{x}},{\bm{y}})$
$\displaystyle=\frac{\partial
f({\bm{x}},{\bm{y}})}{\partial{\bm{x}}}\in\mathbb{F}^{N},$
$\displaystyle\nabla_{{\bm{x}}{\bm{y}}}f({\bm{x}},{\bm{y}})$
$\displaystyle=\left[\frac{\partial^{2}f({\bm{x}},{\bm{y}})}{\partial
x_{i}\partial y_{j}}\right]\in\mathbb{F}^{N\times M},$ (1.1)
$\displaystyle\nabla_{{\bm{x}}{\bm{y}}}f({\bm{\hat{x}}},{\bm{y}})$
$\displaystyle=\nabla_{{\bm{x}}{\bm{y}}}f({\bm{x}},{\bm{y}})\bigg{\rvert}_{{\bm{x}}={\bm{\hat{x}}}}.$

Tables 1.1 and 1.2 summarize our frequently used notation for variables and
functions.

Function |  Description
---|---
$\ell({\bm{\upgamma}})\mapsto\mathbb{R}$ or $\ell({\bm{\upgamma}},{\bm{x}})\mapsto\mathbb{R}$ |  Upper-level loss function used as a fitness measure of $\bm{\upgamma}$. Although $\ell$ is a function of $\bm{\upgamma}$, it is often helpful to write it with two inputs, where typically ${\bm{x}}={\bm{\hat{x}}}$.
$\Upphi({\bm{x}}\,;{\bm{\upgamma}})\mapsto\mathbb{R}$ |  Lower-level cost function used to reconstruct an image.
$R({\bm{x}})\mapsto\mathbb{R}$ |  Regularization function. Incorporates prior information about likely image characteristics.
$d({\bm{x}},{\bm{y}})\mapsto\mathbb{R}$ |  Data-fit term.
$\phi(z)\mapsto\mathbb{R}$ |  Sparsity promoting function, e.g., 0-norm, 1-norm, or corner-rounded 1-norm. Typically used in $R$.
Table 1.2:  Overview of frequently used functions in the review.

###  1.2 Bilevel Method: Set-up and Running Example

This section introduces a generic bilevel problem as well as a specific
bilevel problem that serves as a running example throughout the review. Later
sections discuss many of the ideas presented here more thoroughly. Our hope is
that an early introduction to the formal problem motivates readers and that
this section acts as a quick-reference guide to our notation.

This review considers the image reconstruction problem where the goal is to
form an estimate ${\bm{\hat{x}}}\in\mathbb{F}^{N}$ of a (vectorized) latent
image, given a set of measurements ${\bm{y}}\in\mathbb{F}^{M}$. For denoising
and inpainting problems, $N=M$, but the two dimensions may differ
significantly in more general image reconstruction problems. The forward
operator, ${\bm{A}}\in\mathbb{F}^{M\times N}$ models the physics of the system
such that one would expect ${\bm{y}}={\bm{A}}{\bm{x}}$ in an ideal (noiseless)
system. We focus on linear imaging systems here, but the concepts generalize
readily to nonlinear forward models. When known (in a supervised training
setting), we denote the true, underlying signal as
${\bm{x}}^{\text{true}}\in\mathbb{F}^{N}$.

We focus on model-based image reconstruction methods where the goal is to
estimate $\bm{x}$ from $\bm{y}$ by solving an optimization problem of the form

${\bm{\hat{x}}}={\bm{\hat{x}}}({\bm{\upgamma}})=\operatorname*{argmin}_{{\bm{x}}\in\mathbb{F}^{N}}\Upphi({\bm{x}}\,;{\bm{\upgamma}},{\bm{y}}).$
(1.2)

To simplify notation, we drop $\bm{y}$ from the list of $\Upphi$ arguments
except where needed for clarity. The quality of the estimate ${\bm{\hat{x}}}$
can depend greatly on the choice of the hyperparameters $\bm{\upgamma}$.
Historically there have been numerous approaches pursued for choosing
$\bm{\upgamma}$, such as cross validation [stone:78:cva], generalized cross
validation [golub:79:gcv], the discrepancy principle [phillips:62:atf] and
Bayesian methods [saquib:98:mpe], among others.

Bilevel methods provide a framework for choosing hyperparameters. Without
further ado, a bilevel problem for learning hyperparameters $\bm{\upgamma}$
has the following “double minimization” form:

$\displaystyle\hat{{\bm{\upgamma}}}=\operatorname*{argmin}_{{\bm{\upgamma}}\in\mathbb{F}^{R}}$
$\displaystyle\underbrace{\ell({\bm{\upgamma}}\,;\,{\bm{\hat{x}}}({\bm{\upgamma}}))}_{\ell({\bm{\upgamma}})}\text{
where }$ (UL)
$\displaystyle{\bm{\hat{x}}}({\bm{\upgamma}})=\operatorname*{argmin}_{{\bm{x}}\in\mathbb{F}^{N}}\Upphi({\bm{x}}\,;{\bm{\upgamma}}).$
(LL)

Fig. 1.1 depicts a generic bilevel problem for image reconstruction. The
upper-level (UL), or outer, loss function,
$\ell:\mathbb{R}^{R}\times\mathbb{F}^{N}\mapsto\mathbb{R}$, quantifies how
(not) good is a vector $\bm{\upgamma}$ of learnable parameters; it depends on
the solution to the lower-level (LL) cost function, $\Upphi$, which depends on
$\bm{\upgamma}$. We will also write the upper-level loss function with a
single parameter as $\ell({\bm{\upgamma}})$. We write the lower-level cost as
an optimization problem with “argmin” and thus implicitly assume that $\Upphi$
has unique minimizer with respect to $\bm{x}$. (See Section 4 for more
discussion of this point).

Figure 1.1:  Depiction of a typical bilevel problem for image reconstruction,
illustrated using XCAT phantom from [segars:10:4xp]. The upper box represents
the training process, with the upper-level loss and lower-level cost function.
During training, one minimizes the upper-level loss with respect to a vector
of parameters, $\bm{\upgamma}$, that are used in the image reconstruction
task. Once learned, $\hat{{\bm{\upgamma}}}$ is typically deployed in the same
image reconstruction task, shown in the lower box.

Bilevel methods typically use training data. Specifically, one often assumes
that a given a set of $J$ good quality images
${\bm{x}}^{\text{true}}_{1},\ldots,{\bm{x}}^{\text{true}}_{J}\in\mathbb{F}^{N}$
are representative of the images of interest in a given application. (For
simplicity of notation we assume the training images have the same size, but
they can have different sizes in practice.) We generate corresponding
simulated measurements for each training image using the imaging system model:

${\bm{y}}_{j}={\bm{A}}{\bm{x}}^{\text{true}}_{j}+{\bm{n}}_{j},\quad
j=1,\ldots,J,$ (1.3)

where ${\bm{n}}_{j}\in\mathbb{F}^{M}$ denotes an appropriate random noise
realization111 A more general system model allows the noise to depend on the
data and system model, i.e., ${\bm{n}}_{j}({\bm{A}},{\bm{x}}_{j})$. This
generality is needed for applications with certain noise distributions such as
Poisson noise. . In (1.3), we add one noise realization to each of the $J$
images; in practice one could add multiple noise realizations to each
${\bm{x}}^{\text{true}}_{j}$ to augment the training data. We then use the
training pairs $({\bm{x}}^{\text{true}}_{j},{\bm{y}}_{j})$ to learn a good
value of $\bm{\upgamma}$. After those parameters are learned, subsequent test
images are reconstructed using (1.2) with the learned hyperparameters,
$\hat{{\bm{\upgamma}}}$.

An alternative to (UL) is the following stochastic formulation of bilevel
learning:

$\displaystyle\hat{{\bm{\upgamma}}}=\operatorname*{argmin}_{{\bm{\upgamma}}\in\mathbb{F}^{R}}$
$\displaystyle\underbrace{\mathbb{E}\left[\ell({\bm{\upgamma}})\right]}_{\approx\frac{1}{J}\sum_{j=1}^{J}\ell({\bm{\upgamma}}\,;\,{\bm{\hat{x}}}_{j}({\bm{\upgamma}}))}\text{
where }$ (1.4)
$\displaystyle{\bm{\hat{x}}}_{j}({\bm{\upgamma}})=\operatorname*{argmin}_{{\bm{x}}\in\mathbb{F}^{N}}\Upphi({\bm{x}}\,;{\bm{\upgamma}},{\bm{y}}_{j}).$
(1.5)

The expectation, taken with respect to the training data and noise
distributions, is typically approximated as a sum over the training data.

To offer a concrete example, this review will frequently refer to the
following running example (Ex), a filter learning bilevel problem:

$\displaystyle\hat{{\bm{\upgamma}}}=$
$\displaystyle\operatorname*{argmin}_{{\bm{\upgamma}}\in\mathbb{F}^{R}}\dfrac{1}{2}\|{\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}\|^{2}_{2},\text{
where }$ $\displaystyle{\bm{\hat{x}}}({\bm{\upgamma}})$
$\displaystyle=\operatorname*{argmin}_{{\bm{x}}\in\mathbb{F}^{N}}\dfrac{1}{2}\left\lVert{\bm{A}}{\bm{x}}-{\bm{y}}\right\rVert^{2}_{2}+e^{\beta_{0}}\sum_{k=1}^{K}e^{\beta_{k}}{\bm{1}}^{\prime}\phi.({\bm{c}}_{k}\circledast{\bm{x}};\epsilon),$
(Ex)

where ${\bm{\upgamma}}\in\mathbb{F}^{R}$ contains all variables that we wish
to learn: the filter coefficients ${\bm{c}}_{k}\in\mathbb{F}^{S}$ and tuning
parameters $\beta_{k}\in\mathbb{R}$ for all $k\in[1,K]$. We include an
auxiliary tuning parameter, $\beta_{0}\in\mathbb{R}$, for easier comparison to
other models. Fig. 1.2 depicts the running example. In (Ex), we drop the sum
over $J$ training images for simplicity; the methods easily extend to multiple
training signals. In general, the filters may be of different lengths with
minimal impact on the methods presented in this review. However, for ease of
notation, we will consider ${\bm{c}}_{k}$ to be of length $S$ for all $k$,
e.g., a 2D filter might be $\sqrt{S}\times\sqrt{S}$.

Figure 1.2: Example bilevel problem in (Ex). The vector of learnable
hyperparameters, $\bm{\upgamma}$, includes the tuning parameters, $\beta_{k}$,
and the filter coefficients, ${\bm{c}}_{k}$, shown as example filters.
Although this review will generally consider learning filters of a single
size, the figure depicts how the framework easily extends to 2d filters of
different sizes.

The function $\phi$ in (Ex) is a sparsity-promoting function. If we were to
choose $\phi(z)=|z|$, then the regularizer would involve 1-norm terms of the
type common in compressed sensing formulations:

${\bm{1}}^{\prime}\phi.({\bm{c}}_{k}\circledast{\bm{x}})=\left\lVert{\bm{c}}_{k}\circledast{\bm{x}}\right\rVert_{1}.$

However, to satisfy differentiability assumptions (see Section 4), this review
will often consider $\phi$ to denote the following “corner rounded” 1-norm
having the shape of a hyperbola with the corresponding first and second
derivative:

$\displaystyle\phi(z)$ $\displaystyle=\sqrt{z^{2}+\epsilon^{2}}$ (CR1N)
$\displaystyle{{\dot{\phi}}}(z)$
$\displaystyle=\frac{z}{\sqrt{z^{2}+\epsilon^{2}}}\in[0,1)$
$\displaystyle{\ddot{\phi}}(z)$
$\displaystyle=\frac{\epsilon^{2}}{\left(z^{2}+\epsilon^{2}\right)^{3/2}}\in(0,\frac{1}{\epsilon}],$

where $\epsilon$ is a small, relative to the expected range of $z$, parameter
that controls the amount of corner rounding. (Here, we use a dot over the
function rather than $\nabla$ to indicate a derivative because $\phi$ is a
scalar function.)

###  1.3 Conclusion

Bilevel methods for selecting hyperparameters offer many benefits. Previous
papers motivate them as a principled way to approach hyperparameter
optimization [holler:2018:bilevelapproachparameter,
dempe:2020:bileveloptimizationadvances], as a task-based approach to learning
[peyre:2011:learninganalysissparsity,
haber:2003:learningregularizationfunctionals,
delosreyes:2017:bilevelparameterlearning], and/or as a way to combine the
data-driven improvements from learning methods with the theoretical guarantees
and explainability provided by cost function-based approaches
[chen:2020:learnabledescentalgorithm,
calatroni:2017:bilevelapproacheslearning, kobler:2020:totaldeepvariation]. A
corresponding drawback of bilevel methods are their computational cost; see
Sections 4 and 5 for further discussion.

The task-based nature of bilevel methods is a particularly important
advantage; Section 7.3 exemplifies why by comparing the bilevel problem to
single-level, non-task-based approaches for learning sparsifying filters.
Task-based refers to the hyperparameters being learned based on how well they
work in the lower-level cost function–the image reconstruction task in our
running example. The learned hyperparameters also adapt to the training
dataset and noise characteristics. The task-based nature yields other
benefits, such as making constraints or regularizers on the hyperparameters
generally unnecessary; Section 6.2 presents some exceptions and
[dempe:2020:bileveloptimizationadvances] discusses further bilevel methods for
applications with constraints.

There are three main elements to a bilevel approach. First, the lower-level
cost function in a bilevel problem defines a goal, such as image
reconstruction, including what hyperparameters can be learned, such as filters
for a sparsifying regularizer. Section 2 provides background on this element
specifically for image reconstruction tasks, such as the one in (Ex). Section
6.1 reviews example tasks used in bilevel methods.

Second, the upper-level loss function determines how the hyperparameters
should be learned. While the mean square error (MSE) loss function in the
running example is a common choice, Section 3 discusses other loss functions
based on supervised and unsupervised image quality metrics. Section 6.2 then
reviews example loss functions used in bilevel methods.

While not apparent in the written optimization problem, the third main element
for a bilevel problem is the optimization approach. Section 3.2 briefly
discusses various hyperparameter optimization strategies, then Sections 4 and
5 present multiple gradient-based bilevel optimization strategies. Throughout
the review, we refer to the running example to show how the bilevel
optimization strategies apply.

##  Chapter 2 Background: Image Reconstruction

This review focuses on bilevel problems having image reconstruction as the
lower-level problem. Image reconstruction involves undoing any transformations
inherent in an imaging system, e.g., a camera or CT scanner, and removing
measurement noise, e.g., thermal and shot noise, to realize an image that
captures an underlying object of interest, e.g., a patient’s anatomy. Fig. 2.1
shows an example image reconstruction pipeline for CT data. The following
sections formally define image reconstruction, discuss why regularization is
important, and overview common approaches to regularization.

Figure 2.1: Example image reconstruction pipe-line, illustrated using XCAT
phantom from [segars:10:4xp]. Here $\mathcal{A}$ denotes the actual physical
mapping of the imaging system and $\bm{A}$ denotes the numerical system matrix
used for reconstruction.

###  2.1 Image Reconstruction

Although the true object is in continuous space, image reconstruction is
almost always performed on sampled, discretized signals [lewitt:03:oom].
Without going into detail of the discretization process, we define
${\bm{x}}^{\text{true}}\in\mathbb{F}^{N}$ as the “true,” discrete signal. The
goal of image reconstruction is to recover an estimate
${\bm{\hat{x}}}\approx{\bm{x}}^{\text{true}}$ given corrupted measurements
${\bm{y}}\in\mathbb{F}^{N}$. Although we define the signal as a one-
dimensional vector for notational convenience, the mathematics generalize to
arbitrary dimensions.

To find ${\bm{\hat{x}}}$, image reconstruction involves minimizing a cost
function, $\Upphi({\bm{x}}\,;{\bm{\upgamma}})$, with two terms:

$\displaystyle{\bm{\hat{x}}}=\operatorname*{argmin}_{{\bm{x}}\in\mathbb{F}^{N}}\underbrace{\overbrace{d({\bm{x}}\,;{\bm{y}})}^{\text{Data-
fit}}+\;\;\;\beta\overbrace{R({\bm{x}}\,;{\bm{\upgamma}})}^{\text{Regularizer}}}_{\Upphi({\bm{x}}\,;{\bm{\upgamma}})}$
(2.1)

The first term, $d({\bm{x}}\,;{\bm{y}})$, is a data-fit term that captures the
physics of the ideal (noiseless) system using the matrix
${\bm{A}}\in\mathbb{F}^{M\times N}$; that matrix models the physical system
such that we expect an observation, $\bm{y}$, to be
${\bm{y}}\approx{\bm{A}}{\bm{x}}$. The most common data-fit term penalizes the
square Euclidean norm of the “measurement error,”
$d({\bm{x}}\,;{\bm{y}})=\left\lVert{\bm{A}}{\bm{\hat{x}}}-{\bm{y}}\right\rVert^{2}_{2}$.
This intuitive data-fit term can be derived from a maximum likelihood
perspective, assuming a white Gaussian noise distribution [elad:07:avs].

Many image reconstruction problems have linear system models. In image
denoising problems, one takes ${\bm{A}}=\bm{I}$. For image inpainting,
$\bm{A}$ is a diagonal matrix of 1’s and 0’s, where the 0’s correspond to
sample indices of missing data [guillemot:14:iio]. In MRI, the system matrix
is often approximated as a diagonal matrix times a discrete Fourier transform
matrix, though more accurate models are often needed [fessler:10:mbi]. In some
settings one can learn $\bm{A}$ [golub:80:aao], or at least parts of $\bm{A}$
[ying:07:jir], as part of the estimation process. Although the bilevel method
generalizes to learning $\bm{A}$, the majority of papers in the field assume
$\bm{A}$ is known; Section 6 discusses a few exceptions.

Using the system model (1.3), if $\bm{n}$ were known and $\bm{A}$ were
invertible, we could simply compute
${\bm{\hat{x}}}={\bm{x}}^{\text{true}}={\bm{A}}^{\text{-}1}({\bm{y}}-{\bm{n}})$.
However, $\bm{n}$ is random and, while we may be able to model its
characteristics, we never know it exactly. Further, the system matrix,
$\bm{A}$, is often not invertible because the reconstruction problem is
frequently under-determined, with fewer knowns than unknowns ($M<N$).
Therefore, we must include prior assumptions about ${\bm{x}}^{\text{true}}$ to
make the problem feasible. These assumptions about ${\bm{x}}^{\text{true}}$
are captured in the second, regularization term in (2.1), which depends on
$\bm{\upgamma}$. The following section further discusses regularizers.

In sum, image reconstruction involves finding ${\bm{\hat{x}}}$ that matches
the collected data and satisfies a set of prior assumptions. The data-fit term
encourages ${\bm{\hat{x}}}$ to be a good match for the data; without this
term, there would be no need to collect data. The regularization term
encourages ${\bm{\hat{x}}}$ to match the prior assumptions. Finally, the
tuning parameter, $\beta$, controls the relative importance of the two terms.
The cost function can be minimized using different optimization techniques
depending on the form of each term.

This section is a very short overview of image reconstruction methods. See
[mccann:2019:biomedicalimagereconstruction] for a more thorough review of
biomedical image reconstruction.

###  2.2 Sparsity-Based Regularizers

The regularization, or prior assumption, term in (2.1) often involves
assumptions about sparsity [eldar:12:cs]. The basic idea behind sparsity-based
regularization is that the true signal is sparse in some representation, while
the noise or corruption is not. Thus, one can use the representation to
separate the noise and signal, and then keep only the sparse signal component.

The design problem therefore becomes determining what representation best
sparsifies the signal. There are two main types of sparsity-based regularizers
corresponding to two representational assumptions: synthesis and analysis
[elad:07:avs, ravishankar:20:irf]; Fig. 2.2 depicts both. While both are
popular, this review concentrates on analysis regularizers, which are more
widely represented in the bilevel image reconstruction literature. This
section briefly compares the analysis and synthesis formulations. Here we
simplify the formulas by considering ${\bm{A}}=\bm{I}$; the discussion
generalizes to reconstruction by including $\bm{A}$. For more thorough
discussions of analysis and synthesis regularizers, see [elad:07:avs,
nam:2013:cosparseanalysismodel, ravishankar:20:irf].

Figure 2.2: Depiction of synthesis and analysis sparsity. Under the synthesis
model of sparsity (left), $\bm{x}$ is a linear combination of a few dictionary
atoms. The dictionary, $\bm{D}$, is typically wide, with more atoms (columns)
than elements in $\bm{x}$. Under the analysis model of sparsity (right),
$\bm{x}$ is orthogonal to many filters. The filter matrix, $\bm{\Omega}$, is
typically tall, with more filters (rows) than elements in $\bm{x}$.

####  2.2.1 Synthesis Regularizers

Synthesis regularizers model a signal being composed of building blocks, or
“atoms.” Small subsets of the atoms span a low dimensional subspace and the
sparsity assumption is that the signal requires using only a few of the atoms.
More formally, the synthesis model is ${\bm{y}}={\bm{x}}+{\bm{n}}$, where the
signal ${\bm{x}}={\bm{D}}{\bm{z}}$ and $\bm{z}$ is a sparse vector. The
columns of ${\bm{D}}\in\mathbb{F}^{N\times K}$ contain contain the $K$
dictionary atoms and form a low dimensional subspace for the signal. If
$\bm{D}$ is a wide matrix ($N<K$), the dictionary is over-complete and it is
easier to represent a wide range of signals with a given number of dictionary
atoms. The dictionary is complete when $\bm{D}$ is square (and full rank) and
under-complete if $\bm{D}$ is tall (an uncommon choice).

Assuming one knows or has already learned $\bm{D}$, one can use the sparsity
synthesis assumption to denoise a noisy signal $\bm{y}$ by optimizing

$\displaystyle{\bm{\hat{x}}}={\bm{D}}\cdot(\underbrace{\operatorname*{argmin}_{\bm{z}}\dfrac{1}{2}\left\lVert{\bm{y}}-{\bm{D}}{\bm{z}}\right\rVert^{2}+\phi({\bm{z}})}_{\hat{{\bm{z}}}}).$
(2.2)

The estimation procedure involves finding the sparse codes, $\hat{{\bm{z}}}$,
from which the image is synthesized via
${\bm{\hat{x}}}={\bm{D}}\hat{{\bm{z}}}$. Common sparsity-inducing functions,
$\phi$, are the 0-norm and 1-norm. The 2-norm is occasionally used as $\phi$,
but it does not yield true sparse codes and it over-penalizes large values
[elad:10].

As written in (2.2), the synthesis formulation constrains the signal,
$\bm{x}$, to be in the range of $\bm{D}$. This “strict synthesis” model can be
undesirable in some applications, e.g., when one is not confident in the
quality of the dictionary. An alternative formulation is

$\displaystyle{\bm{\hat{x}}}$
$\displaystyle=\operatorname*{argmin}_{\bm{x}}\dfrac{1}{2}\left\lVert{\bm{y}}-{\bm{x}}\right\rVert^{2}+\beta
R({\bm{x}}),$ $\displaystyle R({\bm{x}})$
$\displaystyle=\min_{\bm{z}}\dfrac{1}{2}\left\lVert{\bm{x}}-{\bm{D}}{\bm{z}}\right\rVert^{2}+\phi({\bm{z}}),$
(2.3)

which no longer constrains $\bm{x}$ to be exactly in the range of $\bm{D}$.
One can also learn $\bm{D}$ while solving (2.3) [peyre:11:aro].

Both synthesis denoising forms have equivalent sparsity constrained versions;
one can replace $\phi$ a characteristic function that is 0 within some desired
set and infinite outside it, e.g.,

$\displaystyle\phi({\bm{z}})=\begin{cases}0&\text{ if
}\left\lVert{\bm{z}}\right\rVert_{0}\leq\alpha\\\ \infty&\text{ else,
}\end{cases}$ (2.4)

for some sparsity constraint given by the hyperparameter
$\alpha\in\mathbb{N}$.

See [candes:2006:robustuncertaintyprinciples, elad:10] for discussions of when
the synthesis model can guarantee accurate recovery of signals. The
minimization problem in (2.3) is called sparse coding and is closely related
to the LASSO problem [tibshirani:1996:regressionshrinkageselection]. One can
think of the entire dictionary $\bm{D}$ as a hyperparameter that can be
learned with a bilevel method [zhou:17:bmb].

####  2.2.2 Analysis Regularizers

Analysis regularizers model a signal as being sparsified when mapped into
another vector space by a linear transformation, often represented by a set of
filters. More formally, an analysis model assumes the signal satisfies
${\bm{\Omega}}{\bm{x}}={\bm{z}}$ for a sparse coefficient vector $\bm{z}$.
Often the rows of the matrix ${\bm{\Omega}}\in\mathbb{F}^{K\times N}$ are
thought of as filters and the rows of $\bm{\Omega}$ where
$[{\bm{\Omega}}{\bm{x}}]_{k}=0$ span a subspace to which $\bm{x}$ is
orthogonal. The analysis operator is called over-complete if $\bm{\Omega}$ is
tall ($N<K$), complete if $\bm{\Omega}$ is square (and full rank), and under-
complete if $\bm{\Omega}$ is wide.

A particularly common analysis regularizer is based on a discretized version
of total variation (TV) [rudin:92:ntv], and uses finite difference filters
(or, more generally, filters that approximate higher-order derivatives). The
finite difference filters sparsify any piece-wise constant (flat) regions in
the signal, leaving the edges that are often approximately sparse in natural
images. Other common analysis regularizers include the discrete Fourier
transform (DFT), curvelets, and wavelet transforms
[candes:2011:compressedsensingcoherent].

The literature is less consistent in analysis regularizer vocabulary, and
$\bm{\Omega}$ has been called an analysis dictionary, an analysis operator, a
filter matrix, and a cosparse operator. The term “cosparse” comes from the
sparsity holding in the codomain of the transformation
$T\\{{\bm{x}}\\}={\bm{\Omega}}{\bm{x}}$. The cosparsity of $\bm{x}$ with
respect to $\bm{\Omega}$ is the number of zeros in ${\bm{\Omega}}{\bm{x}}$ or
$K-\left\lVert{\bm{\Omega}}{\bm{x}}\right\rVert_{0}$
[nam:2013:cosparseanalysismodel]. Correspondingly, “cosupport” describes the
indices of the rows where ${\bm{\Omega}}{\bm{x}}=0$. We find the phrase
“analysis operator” intuitive for general $\bm{\Omega}$’s and “filter matrix”
more descriptive when referring to the specific (common) case when the rows of
$\bm{\Omega}$ are dictated by a set of convolutional filters.

Assuming one knows, or has already learned, $\bm{\Omega}$, one can use the
analysis sparsity assumption to denoise a noisy signal, $\bm{y}$, by
optimizing

${\bm{\hat{x}}}=\operatorname*{argmin}_{\bm{x}}\dfrac{1}{2}\left\lVert{\bm{y}}-{\bm{x}}\right\rVert^{2}+\beta\phi({\bm{\Omega}}{\bm{x}}).$
(2.5)

An alternative version is

$\displaystyle{\bm{\hat{x}}}=\operatorname*{argmin}_{\bm{x}}\dfrac{1}{2}\left\lVert{\bm{y}}-{\bm{x}}\right\rVert^{2}+\beta
R({\bm{x}})$ (2.6) $\displaystyle
R({\bm{x}})=\min_{\bm{z}}\dfrac{1}{2}\left\lVert{\bm{\Omega}}{\bm{x}}-{\bm{z}}\right\rVert^{2}+\phi({\bm{z}}).$

As in the synthesis case, both analysis formulations have equivalent sparsity
constrained forms using a $\phi$ as in (2.4).

See [candes:2011:compressedsensingcoherent] for an error bound on the
estimated signal ${\bm{\hat{x}}}$ when using a 1-norm as the sparsity
function.

####  2.2.3 Comparing Analysis and Synthesis Approaches

The analysis and synthesis models are equivalent when the dictionary and
analysis operator are invertible, with ${\bm{D}}={\bm{\Omega}}^{\text{-}1}$
[elad:07:avs]. Furthermore, in the denoising scenario where the system matrix
$\bm{A}$ is identity, the two are almost equivalent in the under-complete case
[elad:07:avs], with the lack of full equivalence stemming from the analysis
form not constraining $\bm{x}$ to be in the range space $\bm{D}$. However, the
two models diverge in the more common over-complete case.

Whether analysis-based or synthesis-based regularizers are preferable is an
open question, and the answer likely depends on the application and the
relative importance of reconstruction accuracy and speed [elad:07:avs].
Synthesis regularization is perhaps easier to interpret because of its
constructive nature, and the synthesis approach used to be “widely considered
to provide superior results” [elad:07:avs, 950]. However, [elad:07:avs] goes
on to show that an analysis regularizer produced more accurate reconstructed
images in experiments on real images. Later analysis-based results also show
competitive, if not superior, quality results when compared to similar
synthesis models [hawe:13:aol,
ravishankar:2013:learningsparsifyingtransforms]. See [fessler:20:omf] for a
survey of optimization methods for MRI reconstruction and a comparison of the
computational challenges for cost functions with synthesis and analysis-based
regularizers.

The analysis and synthesis regularizers in (2.2) and (2.6) quickly yield
infeasibly large operators as the signal size increases. In practice, both
approaches are usually implemented with patch-based formulations. For the
synthesis approach, the patches typically overlap and there is an averaging
effect. Analysis regularizers that have rows corresponding to filters, called
the convolutional analysis model, extend very naturally to a global image
regularizer. For example, in the lower-level cost function of our running
filter learning example (Ex), we can define an analysis regularizer matrix as
follows:

${\bm{\Omega}}=\begin{bmatrix}{\bm{C}}_{1}\\\ \vdots\\\
{\bm{C}}_{K}\end{bmatrix}\in\mathbb{F}^{KN\times N}.$

See [chen:2014:insightsanalysisoperator] and [pfister:2019:learningfilterbank]
for discussion of the connections between global models and patch-based models
for analysis regularizers. This survey focuses on bilevel learning of analysis
regularizers.

###  2.3 Brief History of Analysis Regularizer Learning

In 2003, haber:2003:learningregularizationfunctionals proposed using bilevel
methods to learn part of the regularizer in inverse problems. The authors
motivate the use of bilevel methods through the task-based nature, nothing
that “the choice of good regularization operators strongly depends on the
forward problem.” They consider learning tuning parameters, space-varying
weights, and regularization operators (comparable to defining $\phi$), all for
regularizers based on penalizing the energy in the derivatives of the
reconstructed image. Their framework is general enough to handle learning
filters. Ref. [haber:2003:learningregularizationfunctionals] was published a
few years earlier than the other bilevel methods we consider in this review
and was not cited in most other early works;
[afkham:2021:learningregularizationparameters] calls it a “groundbreaking, but
often overlooked publication.”

In 2005, roth:2005:fieldsexpertsframework proposed the Field of Experts (FoE)
model to learn filters. Many papers on bilevel methods for filter learning
cite FoE as a starting or comparison point. The FoE model is a translation-
invariant analysis operator model, built on convolutional filters. It is
motivated by the local operators and presented as a Markov random field model,
with the order of the field determined by the filter size.

Under the FoE model, the negative log111 By taking the log of the probability
model in [roth:2005:fieldsexpertsframework], the connection between the FoE
and the regularization term in the lower-level of the running filter learning
example (Ex) is more evident. of the probability of a full image, $\bm{x}$, is
proportional to

$\displaystyle\sum_{k}\beta_{k}\;\phi.({\bm{c}}_{k}\circledast{\bm{x}})\text{
where }\phi(z)=\text{log}\left(1+\dfrac{1}{2}z^{2}\right).$ (2.7)

This (non-convex) choice of sparsity function $\phi$ stems from the Student-t
distribution. Ref. [roth:2005:fieldsexpertsframework] learns the filters and
filter-dependent tuning parameters such that the model distribution is as
close as possible (defined using Kullback-Leibler divergence) to the training
data distribution.

In 2007, tappen:2007:learninggaussianconditional proposed a different model
based on convolutional filters: the Gaussian Conditional Random Field (GCRF)
model. Rather than using a sparsity promoting regularizer, the GCRF uses a
quadratic function for $\phi$. The authors introduce space-varying weights,
$\bm{W}$, so that the quadratic model does not overly penalize sharp features
in the image. The general idea behind ${\bm{W}}$ is to use the given (noisy)
image to guess where edges occur, and correspondingly penalize those areas
less to avoid blurring edges. The likelihood for GCRF model is thus (to within
a proportionality constant and monotonic function transformations):

$\displaystyle\sum_{k}\left\lVert{\bm{c}}_{k}\circledast{\bm{x}}-e_{k}\\{{\bm{x}}\\}\right\rVert^{2}_{{\bm{W}}_{k}},$

where the term $e_{k}\\{{\bm{x}}\\}$ captures the estimated value of the
filtered image. For example, [tappen:2007:learninggaussianconditional] used
one averaging filter and multiple differencing filters for the
${\bm{c}}_{k}$’s. The corresponding estimated values are $\bm{x}$ for the
averaging filter and zero for the differencing filters.

The filters, ${\bm{c}}_{k}$, are pre-determined in the GCRF model; the learned
element is how to form the weights as a function of image features.
Specifically, each ${\bm{W}}_{k}$ is formed as a linear combination of the
(absolute) responses to a set of edge-detecting filters, with the linear
combination coefficients learned from training data. Rather than maximizing
the likelihood of training data as in [roth:2005:fieldsexpertsframework],
[tappen:2007:learninggaussianconditional] learns these coefficients to
minimize the (corner-rounded) $l_{1}$ norm of the error of the predicted
image, which is a form of bilevel learning even though not described with that
terminology.

Apparently one of the first papers to explicitly propose using bilevel methods
to learn filters appeared in 2009, where samuel:2009:learningoptimizedmap
considered a bilevel formulation where the upper-level loss was mean square
error on training data and the lower-level cost was a denoising task based on
filter sparsity equivalent to (Ex). The method builds on the FoE model, using
the same $\phi$ as in [roth:2005:fieldsexpertsframework], but now learns the
filters using a bilevel formulation rather than by maximizing a likelihood.

In 2011, peyre:2011:learninganalysissparsity proposed a similar bilevel method
to learn analysis regularizers. The authors generalized the denoising task to
use an analysis operator matrix and a wider class of sparsifying functions.
Their results concentrate on the convolutional filter case with a corner-
rounded 1-norm for $\phi$.

Both [samuel:2009:learningoptimizedmap] and
[peyre:2011:learninganalysissparsity] focus on introducing the bilevel method
for analysis regularizer learning, with denoising or inpainting as
illustrations. Section 4 further discusses the methodology of both papers.
Many of the bilevel based papers in this review build on one or both of their
efforts. The rest of the review will summarize other bilevel based papers;
here, we continue to highlight some of papers in the non-bilevel thread of the
literature for context and comparison.

ophir:2011:sequentialminimaleigenvalues proposed another approach to learning
an analysis operator. The method learns the operator one row at a time by
searching for vectors orthogonal to the training signals. Algorithm parameters
were chosen empirically without an upper-level loss function as a guide.

Between 2011 [yaghoobi:2011:analysisoperatorlearning] and 2013
[yaghoobi:2013:constrainedovercompleteanalysis],
yaghoobi:2011:analysisoperatorlearning were among the first to formally
present analysis operator learning as an optimization problem. Their
conference paper [yaghoobi:2011:analysisoperatorlearning] considered noiseless
training data and proposed learning an analysis operator as

$\operatorname*{argmin}_{\bm{\Omega}}\|{\bm{\Omega}}{\bm{X}}^{\text{true}}\|_{1}\text{
s.t. }{\bm{\Omega}}\in\mathcal{S}$ (2.8)

for some constrained set $\mathcal{S}$. Each column of
${\bm{X}}^{\text{true}}$ contains a training sample. The authors discussed
varying options for $\mathcal{S}$, including a row norm, full rank, and tight
frame constrained set.

Without any constraint on $\bm{\Omega}$, the trivial solution to (2.8) would
be to learn the zero matrix, which is not informative for any problem such as
image denoising. Section 7.3 discusses in more detail the need for constraints
and the various options proposed for filter learning.

Ref. [yaghoobi:2013:constrainedovercompleteanalysis] extends (2.8) to the
noisy case where one does not have access to ${\bm{X}}^{\text{true}}$. The
proposed cost function is

$\operatorname*{argmin}_{{\bm{\Omega}},\,{\bm{X}}}\left\lVert{\bm{\Omega}}{\bm{X}}\right\rVert_{1}+\frac{\beta}{2}\left\lVert{\bm{Y}}-{\bm{X}}\right\rVert^{2}\text{
s.t. }{\bm{\Omega}}\in\mathcal{S},$ (2.9)

where each column of $\bm{Y}$ contains a noisy data vector. Ref.
[yaghoobi:2013:constrainedovercompleteanalysis] minimized (2.9) by alternating
updating $\bm{X}$, using alternating direction method of multipliers (ADMM),
and $\bm{\Omega}$, using a projected subgradient method for various constraint
sets $\mathcal{S}$, especially Parseval tight frames.

Also in 2013, ravishankar:2013:learningsparsifyingtransforms made a
distinction between the analysis model, where one models
${\bm{y}}={\bm{x}}+{\bm{n}}$ with ${\bm{z}}={\bm{\Omega}}{\bm{x}}$ being
sparse, and the transform model, where
${\bm{\Omega}}{\bm{y}}={\bm{z}}+{\bm{n}}$ where $\bm{z}$ is sparse. The
analysis version models the measurement as being a cosparse signal plus noise;
the transform version models the measurement as being approximately cosparse.
Another perspective on the distinction is that, if there is no noise, the
analysis model constrains $\bm{y}$ to be in the range space of $\bm{\Omega}$,
while there is no such constraint on the transform model. The corresponding
transform learning problem is

$\displaystyle\operatorname*{argmin}_{{\bm{\Omega}}}\min_{\bm{Z}}\left\lVert{\bm{\Omega}}{\bm{Y}}-{\bm{Z}}\right\rVert^{2}_{2}+$
$\displaystyle R({\bm{\Omega}})\quad\text{ s.t.
}\left\lVert{\bm{Z}}_{i}\right\rVert_{0}\leq\alpha\;\forall i,$ (2.10)

where $i$ indexes the columns of $\bm{Z}$. Ref.
[ravishankar:2013:learningsparsifyingtransforms] considers only square
matrices $\bm{\Omega}$. The regularizer, $R$, promotes diversity in the rows
of $\bm{\Omega}$ to avoid trivial solutions, similar to the set constraint in
(2.9).

A more recent development is directly modeling the convolutional structure
during the learning process. In 2020,
[chun:2020:convolutionalanalysisoperator] proposed Convolutional Analysis
Operator Learning (CAOL) to learn convolutional filters without patches. The
CAOL cost function is

$\displaystyle\operatorname*{argmin}_{[{\bm{c}}_{1},\ldots,{\bm{c}}_{K}]}\sum_{k=1}^{K}\min_{\bm{z}}\dfrac{1}{2}\left\lVert{\bm{c}}_{k}\circledast{\bm{x}}-{\bm{z}}\right\rVert^{2}_{2}+\beta\left\lVert{\bm{z}}\right\rVert_{0}\text{
s.t. }[{\bm{c}}_{1}\ldots{\bm{c}}_{K}]\in\mathcal{S}.$ (2.11)

Unlike the previous cost functions, which typically require patches, CAOL can
easily handle full-sized training images $\bm{x}$ due to the nature of the
convolutional operator.

While model-based methods were being developed in the signal processing
literature, convolutional neural network (CNN) models were being advanced and
trained in the machine learning and computer vision literature [haykin:96:nne]
[hwang:97:tpp] [lucas:18:udn]. The filters used in CNN models like U-Nets
[ronneberger:15:unc] can be thought of as having analysis roles in the earlier
layers, and synthesis roles in the final layers [ye:18:dcf]. See also
[wen:20:tlf] for further connections between analysis and transform models
within CNN models. CNN training is usually supervised, and the supervised
approach of bilevel learning of filters strengthens the relationships between
the two approaches. A key distinction is that CNN models are generally feed-
forward computations, whereas bilevel methods of the form (LL) have a cost
function formulation. See Section 7 for further discussion of the parallels
between CNNs and bilevel methods.

###  2.4 Summary

This background section focused on the lower-level problem: image
reconstruction with a sparsity-based regularizer. After defining the problem
and the need for regularization, Section 2.3 reviewed the history of analysis
regularizer learning and included many examples of methods to learn
hyperparameters.

Bilevel methods are just one, task-based way to learn such hyperparameters.
Section 7.3 further expands on this point, but we can already see benefits of
the task-based nature of bilevel methods. Without the bilevel approach,
filters are learned so that they best sparsify training data. These
sparsifying filters can then be used in a prior for image reconstruction
tasks. However, they are learned to sparsify, not necessarily to best
reconstruct. In contrast, the bilevel approach aims to learn filters that best
reconstruct images (or whatever other task is desired), even if those filters
are not the ones that best sparsify. Although this may seem like a subtle
distinction, [chambolle:2021:learningconsistentdiscretizations] shows that
different (synthesis) filters work better for image denoising versus image
inpainting.

Having provided some background on the lower-level cost function and motivated
bilevel methods, this review now turns to defining the upper-level loss
function and surveying methods of hyperparameter optimization.

##  Chapter 3 Background: Hyperparameter Optimization

Most inverse problems involve at least one hyperparameter. For example, the
general reconstruction cost function (2.1) requires choosing the tuning
parameter $\beta$ that trades-off the influence of the data-fit and
regularization terms. The field of hyperparameter optimization is quite large
and can include categorical hyperparameters, such as which optimizer to use;
conditional hyperparameters, where certain hyperparameters are relevant only
if others take on certain values; and integer or real-valued hyperparameters
[feurer:2019:chapterhyperparameteroptimization]. Here, we focus on learning
real-valued, continuous hyperparameters.

Figure 3.1:  Example reconstructed simulated MRI images that demonstrate the
importance of tuning parameters. (a) The original image,
${\bm{x}}^{\text{true}}\in R^{N}$, is a SheppLogan phantom [shepp:74:tfr] and
$N$ is the number of pixels. (b) A simplistic reconstruction
$\frac{1}{N}{\bm{A}}^{\prime}{\bm{y}}$ of the noisy, undersampled data,
$\bm{y}$. This image is used as initialization, ${\bm{x}}^{(0)}$, for the
following reconstructions. (c-e) Reconstructed images, found by optimizing
$\operatorname*{argmin}_{\bm{x}}\dfrac{1}{2}\|{\bm{A}}{\bm{x}}-{\bm{y}}\|^{2}_{2}+10^{\beta}N\phi({\bm{C}}{\bm{x}})$,
where $\bm{C}$ is an operator that takes vertical and horizontal finite
differences. The reconstructed images correspond to (c) $\beta=-6$, resulting
in an image that contains ringing artifacts, (d) $\beta=-3$, resulting in a
visually appealing ${\bm{\hat{x}}}$, and (e) $\beta=1$, resulting in a blurred
image. The demonstration code and more details about the reconstruction set-up
are available on github [fessler:2020:mirtdemo01recon].

A hyperparameter’s value can greatly influence the properties of the minimizer
and a tuned hyperparameter typically improves over a default setting
[feurer:2019:chapterhyperparameteroptimization]. Fig. 3.1 illustrates how
changing a tuning parameter can dramatically impact the visual quality of the
reconstructed image. If $\beta$ is too low, not enough weight is on the
regularization term, and the minimizer is likely to be corrupted by noise in
the measurements. If $\beta$ is too high, the regularization term dominates,
and the minimizer will not align with the measurements.

Generalizing to an arbitrary learning problem that could have multiple
hyperparameters, the goal of hyperparameter optimization is to find the “best”
set of hyperparameters, $\hat{{\bm{\upgamma}}}$ to meet a goal, described by a
loss function $\ell$. Specifically, we wish to solve

$\hat{{\bm{\upgamma}}}=\operatorname*{argmin}_{{\bm{\upgamma}}\in\Gamma}\mathbb{E}\left[\ell({\bm{\upgamma}})\right],$
(3.1)

where $\Gamma$ is the set of all possible hyperparameters and the expectation
is taken with respect to the distribution of the input data. If evaluating
$\ell$ uses the output of another optimization problem, e.g.,
${\bm{\hat{x}}}$, then (3.1) is a bilevel problem as defined in (UL).

There are two key tasks in hyperparameter optimization. The first is to
quantify how good a hyperparameter is; this step is equivalent to defining
$\ell$ in (3.1). The second step is finding a good hyperparameter, which is
equivalent to designing an optimization algorithm to minimize (3.1). The next
two sections address each of these tasks in turn.

###  3.1 Image Quality Metrics

This section concentrates on the part of the upper-level loss function that
compares the reconstructed image, ${\bm{\hat{x}}}({\bm{\upgamma}})$, to the
true image, ${\bm{x}}^{\text{true}}$. As mentioned in Section 1, bilevel
methods rarely require additional regularization for $\bm{\upgamma}$, but it
is simple to add a regularization term to any of the loss functions if useful
for a specific application. To discuss only the portion of the loss function
that measures image quality, we use the notation
$\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))=l({\bm{\hat{x}}},\,{\bm{x}}^{\text{true}})$.

Picking a loss function is part of the engineering design process. No single
loss function is likely to work in all scenarios; users must decide on the
loss function that best fits their system, data, and goals. Consequently,
there are a wide variety of loss functions proposed in the literature and some
approaches combined multiple loss functions
[you:2018:structuresensitivemultiscaledeep,
hammernik:2020:machinelearningimage].

One important decision criteria when selecting a loss function is the end
purpose of the image. Much of the image quality assessment (IQA) literature
focuses on metrics for images of natural scenes and is often motivated by
applications where human enjoyment is the end-goal [wang:04:iqa,
wang:2011:reducednoreferenceimage]. In contrast, in the medical image
reconstruction field, image quality is not the end-goal, but rather a means to
achieving a correct diagnosis. Thus, the perceptual quality is less important
than the information content.

There are two major classes of image quality metrics in the IQA literature,
called full reference and no reference IQA111There are also reduced-reference
image quality metrics, but we will not consider those here.. The principles
are somewhat analogous to supervised and unsupervised approaches in the
machine learning literature. This section discusses some of the most common
full reference and no reference loss functions; see
[zhang:2012:comprehensiveevaluationfull] for a comparison of 11 full-reference
IQA metrics and [zhang:2020:blindimagequality] for additional no-reference IQA
metrics.

Perhaps surprisingly, the bilevel filter learning literature contains few
examples of loss functions other than MSE or slight variants (see Section
6.2). While this is likely at least partially due to the computational
requirements of bilevel methods (see Section 4 and 5), exploring additional
loss functions is an interesting future direction for bilevel research.

####  3.1.1 Full Reference IQA

Full reference IQA metrics assume that you have a noiseless image,
${\bm{x}}^{\text{true}}$, for comparison. Some of the simplest (and most
common) full reference loss functions are:

  * •

Mean square error (or $\ell_{2}$ error):
$l_{\mathrm{MSE}}({\bm{\hat{x}}},{\bm{x}}^{\text{true}})=\frac{1}{N}\left\lVert{\bm{\hat{x}}}-{\bm{x}}^{\text{true}}\right\rVert^{2}_{2}$

  * •

Mean absolute error (or $\ell_{1}$ error):
$l_{\mathrm{MAE}}({\bm{\hat{x}}},{\bm{x}}^{\text{true}})=\frac{1}{N}\left\lVert{\bm{\hat{x}}}-{\bm{x}}^{\text{true}}\right\rVert_{1}$

  * •

Signal to Noise Ratio (SNR, commonly expressed in dB):
$l_{\mathrm{SNR}}({\bm{\hat{x}}},{\bm{x}}^{\text{true}})=10\operatorname{log}\left(\frac{\left\lVert{\bm{x}}^{\text{true}}\right\rVert^{2}_{2}}{\left\lVert{\bm{\hat{x}}}-{\bm{x}}^{\text{true}}\right\rVert^{2}_{2}}\right)$

  * •

Peak SNR (PSNR, in dB):
$l_{\mathrm{PSNR}}({\bm{\hat{x}}},{\bm{x}}^{\text{true}})=10\operatorname{log}\left(\frac{N\left\lVert{\bm{x}}^{\text{true}}\right\rVert_{\infty}}{\left\lVert{\bm{\hat{x}}}-{\bm{x}}^{\text{true}}\right\rVert^{2}_{2}}\right)$.

The Euclidean norm is also frequently used as the data-fit term for
reconstruction.

MSE (and the related metrics SNR and PSNR) are common in the signal processing
field; they are intuitive and easy to use because they are differentiable and
operate point-wise. However, these measures do not align well with human
perceptions of image quality [mason:2020:comparisonobjectiveimage,
zhang:2012:comprehensiveevaluationfull]. For example, scaling an image by 2
leads to the same visual quality but causes 100% MSE. Fig. 3.2 shows a clean
image and five images with different degradations. All five degraded images
have almost equivalent MSEs, but humans judge their qualities as very
different.

Figure 3.2:  Example distortions that yield images with identical normalized
MSE values:
$\left\lVert{\bm{x}}^{\text{true}}-{\bm{x}}\right\rVert/\left\lVert{\bm{x}}^{\text{true}}\right\rVert=0.17$.
(a) The original image, ${\bm{x}}^{\text{true}}$, is a SheppLogan phantom
[shepp:74:tfr]. The remaining images are displayed with the same colormap and
have the following distortions: (b) blurred with an averaging filter, (c)
additive, white Gaussian noise, (d) salt and pepper noise, and (e) a constant
value added to every pixel.

Tuning parameters using MSE tends to lead to images that are overly-smoothed,
sacrificing high frequency information [gholizadehansari:20:dlf, seif:18:ebl].
High frequency details are particularly important for perceptual quality as
they correspond to edges in images. Therefore, some authors use the MSE on
edge-enhanced versions of images to discourage solutions that blur edges. For
example, [ravishankar:2011:mrimagereconstruction] used a “high frequency error
norm” metric consisting of the MSE of the difference of ${\bm{\hat{x}}}$ and
${\bm{x}}^{\text{true}}$ after applying a Laplacian of Gaussian (LoG) filter.

Another common full-reference IQA is Structural SIMilarity (SSIM)
[wang:2004:imagequalityassessment] that attempts to address the issues with
MSE discussed above. SSIM is defined in terms of the local luminance,
contrast, and structure in images. A multiscale extension of SSIM, called MS-
SSIM, considers information at different resolutions
[wang:2003:multiscalestructuralsimilarity]. The method computes the contrast
and structure measures of SSIM for downsampled versions of the input images
and then defines MS-SSIM as the product of the luminance at the original scale
and the contrast and structure measures at each scale. However, SSIM and MS-
SSIM may not correlate well with human observer performance on radiological
tasks [renieblas:17:ssi].

Recent works, e.g., [bosse:2018:deepneuralnetworks,
zhang:2020:blindimagequality], consider using (deep) CNN models for IQA. CNN
methods are increasingly popular and their use as a model for the human visual
system [lindsay:2020:convolutionalneuralnetworks] makes them an attractive
tool for assessing images. For example, [bosse:2018:deepneuralnetworks]
proposed a CNN with convolutional and pooling layers for feature extraction
and fully connected layers for regression. They used VGG
[simonyan:2015:verydeepconvolutional], a frequently-cited CNN design with
$3\times 3$ convolutional kernels, as the basis of the feature extraction
portion of their network. Ref. [bosse:2018:deepneuralnetworks] showed that
deeper networks with more learnable parameters were able to better predict
image quality. However, datasets of images with quality labels remain
relatively scarce, making it difficult to train deep networks.

####  3.1.2 No Reference IQA

No reference, or unsupervised, IQA metrics attempt to quantify an image’s
quality without access to a noiseless version of the image. These metrics rely
on modeling statistical characteristics of images or noise.

Many no reference IQA metrics assume the noise distribution is known. The
discrepancy principle is a classic example that uses an assumed noise
distribution to characterize the expected relation between the reconstructed
image and the noisy data. For additive zero-mean white Gaussian noise with
known variance $\sigma^{2}$, the discrepancy principle uses the fact that the
expected MSE in the data space is the noise variance [phillips:62:atf]:

$\mathbb{E}\left[\frac{1}{M}\left\lVert{\bm{A}}{\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{y}}\right\rVert^{2}_{2}\right]=\sigma^{2}.$

The discrepancy principle can be used as a stopping criteria in machine
learning methods or as a loss function, e.g.,

$\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))=\left(\frac{1}{M}\left\lVert{\bm{A}}{\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{y}}\right\rVert^{2}_{2}-\sigma^{2}\right)^{2}.$

However, images of varying quality can yield the same noise estimate, as seen
in Fig. 3.2. Related methods have been developed for Poisson noise as well
[hebert:92:sbm].

Paralleling MSE’s popularity among supervised loss metrics, Stein’s Unbiased
Risk Estimator (SURE) [stein:1981:estimationmeanmultivariate] is an unbiased
estimate of MSE that does not require noiseless images. Let
${\bm{y}}={\bm{x}}^{\text{true}}+{\bm{n}}$ denote a signal plus noise
measurement where $\bm{n}$ is, as above, Gaussian noise with known variance
$\sigma^{2}$. The SURE estimate of the MSE of a denoised signal,
${\bm{\hat{x}}}$, is

$\displaystyle\frac{1}{N}\left\lVert{\bm{y}}-{\bm{\hat{x}}}({\bm{y}})\right\rVert^{2}_{2}-\sigma^{2}+\frac{2\sigma^{2}}{N}\nabla_{y}{\bm{\hat{x}}}({\bm{y}}),$
(3.2)

where we write ${\bm{\hat{x}}}$ as a function of $\bm{y}$ to emphasize the
dependence. For large signal dimensions, $N$, such as is common in image
reconstruction problems, the law of large numbers suggests SURE is a fairly
accurate approximation of the true MSE.

However, it is often impractical to evaluate the divergence term in (3.2), due
to computational limitations or not knowing the form of
${\bm{\hat{x}}}({\bm{y}})$. A Monte-Carlo approach to estimating the
divergence [ramani:2008:montecarlosureblackbox] uses the following key
equation:

$\displaystyle\nabla_{\bm{y}}{\bm{\hat{x}}}({\bm{y}})=\lim_{\epsilon\rightarrow
0}\mathbb{E}\left[{\bm{b}}^{\prime}\cdot\frac{{\bm{\hat{x}}}({\bm{y}}+\epsilon{\bm{b}})-{\bm{\hat{x}}}({\bm{y}})}{\epsilon}\right],$
(3.3)

where $\bm{b}$ is a independent and identically distributed (i.i.d.) random
vector with zero mean, unit variance, and bounded higher order moments.
Theoretical and empirical arguments show that a single noise vector can well-
approximate the divergence [ramani:2008:montecarlosureblackbox], so only two
calls to the lower-level solver ${\bm{\hat{x}}}({\bm{y}})$ are required. This
method treats the lower-level problem like a blackbox, thus allowing one to
estimate the divergence of complicated functions, including those that may not
be differentiable.

See [soltanayev:2018:trainingdeeplearning, kim:20:uto, zhussip:19:tdl] for
examples of applying the Monte-Carlo estimation of SURE to train deep neural
networks, and [zhang:2020:bilevelnestedsparse,
deledalle:2014:steinunbiasedgradient] for two examples of learning a tuning
parameter using a bilevel approach with SURE as the upper-level loss function.
For extensions to inverse problems and noise from exponential families, see
[eldar:2009:generalizedsureexponential].

While SURE and the discrepancy principle are popular no-reference metrics in
the signal processing literature, there are many additional no-reference
metrics in the image quality assessment literature. These metrics typically
depend on modeling one (or more) of three things
[wang:2011:reducednoreferenceimage]:

  * •

image source characteristics,

  * •

image distortion characteristics, e.g., blocking artifact from JPEG
compression, and/or

  * •

human visual system perceptual characteristics.

As an example of a strategy that can capture both image source and human
visual system characteristics, natural scene222 Natural scenes are those
captured by optical cameras (not created by computer graphics or other
artificial processes) and are not limited to outdoor scenes.  statistics
characterize the distribution of various features in natural scenes, typically
using some filters [mittal:2013:makingcompletelyblind,
wang:2011:reducednoreferenceimage]. If a feature reliably follows a specific
statistical pattern in natural images but has a noticeably different
distribution in distorted images, one can use that feature to assign quality
scores to images. Some IQA metrics attempt to first identify the type of
distortion and measure features specific to that distortion, while others use
the same features for all images.

In addition to their use in full-reference IQA, CNN models have be trained to
perform no-reference IQA [kang:2014:convolutionalneuralnetworks,
bosse:2018:deepneuralnetworks]. For example,
kang:2014:convolutionalneuralnetworks proposes a CNN model that extracts small
(32 $\times$ 32) patches from images, estimates quality for each one, and
averages the scores over all patches to get the quality score for the entire
image. Briefly, their method involves local contrast normalization for each
patch, applying (learned) convolutional filters to extract features, maximum
and minimum pooling, and fully connected layers with rectified linear units
(ReLUs). As with most no reference IQAs,
[kang:2014:convolutionalneuralnetworks] trained their CNN on a dataset of
(human encoded) image quality scores (see
[laboratoryforimagevideoengineering::imagevideoquality] for a commonly used
collection of publicly available test images with quality scores). Unlike most
other IQA approaches, [kang:2014:convolutionalneuralnetworks] used
backpropagation to learn all the CNN weights rather than learning a
transformation from (handcrafted) features to quality scores.

Interestingly, some of the no-reference IQA metrics
[mittal:2013:makingcompletelyblind, kang:2014:convolutionalneuralnetworks,
wang:2011:reducednoreferenceimage] approach the performance of the full-
reference IQAs in terms of their ability to match human judgements of image
quality. This observation suggests that there is room to improve full-
reference IQA metrics and that assessing image quality is a very challenging
problem!

###  3.2 Parameter Search Strategies

After selecting a metric to measure how good a hyperparameter is, the next
task is devising a strategy to find the best hyperparameter according to that
metric. Search strategies fall into three main categories: (i) model-free,
$\ell$-only; (ii) model-based, $\ell$-only; and (iii) gradient-based, using
both $\ell$ and $\nabla\ell$. Model-free strategies do not assume any
information about about the hyperparameter landscape, whereas model-based
strategies use historical $\ell$ evaluations to predict the loss-function at
untested hyperparameter values.

The following sections describe common model-free and model-based
hyperparameter search strategies that only use $\ell$. See
[dempe:2020:bileveloptimizationadvances, Ch. 13 and Ch. 20.6] for discussion
of additional gradient-free methods for bilevel problems, e.g., population-
based evolutionary algorithms, and
[larson:2019:derivativefreeoptimizationmethods] for a general discussion of
derivative-free optimization methods.

The third class of hyperparameter optimization schemes are approaches based on
gradient descent of a bilevel problem. The high-level strategy in bilevel
approaches is to calculate the gradient of the upper-level loss function,
$\ell$, with respect to $\bm{\upgamma}$ and then use any gradient descent
method to minimize $\bm{\upgamma}$. Although this approach can be
computationally challenging, it generalizes well to a large number of
hyperparameters. Section 4 and Section 5 discuss this point further and go
into depth on different methods for computing this gradient.

####  3.2.1 Model-free

The most common search strategy is probably an empirical search, where a
researcher tries different hyperparameter combinations manually. A punny, but
often accurate, term for this manual search is GSD: grad[uate] student descent
[gencoglu:2019:harksidedeep]. bergstra:12:rsf hypothesizes that manual search
is common because it provides some insight as the user must evaluate each
option, it requires no overhead for implementation, and it can perform
reliably in very low dimensional hyperparameter spaces.

Grid search is a more systematic alternative to manual search. When there are
only one or two continuous hyperparameters, or the possible set of
hyperparameters, ${\bm{\Gamma}}$, is small, a grid search (or exhaustive
search) strategy may suffice to find the optimal value,
$\hat{{\bm{\upgamma}}}$, to within the grid spacing. However, the complexity
of grid search grows exponentially with the number of hyperparameters.
Regularizers frequently have many hyperparameters, so one must use a more
sophisticated search strategy.

One popular approach is random search, which bergstra:12:rsf shows is superior
to a grid search, especially when some hyperparameters are more important than
others. There are also variations on random search, such as using Poisson disk
sampling theory to explore the hyperparameter space [muniraju:18:crs]. The
simplicity of random search makes it popular, and, even if one uses a more
complicated search strategy, random search can provide a useful baseline or an
initialization strategy. However, random search, like grid search, suffers
from the curse of dimensionality, and is less effective as the hyperparameter
space grows.

Another group of model-free blackbox strategies are population-based methods
such as evolutionary algorithms. A popular population-based method is the
covariance matrix adaption evolutionary strategy (CMA-ES)
[beyer:2001:theoryevolutionstrategies]. In short, every iteration, CMA-ES
involves sampling a multivariate normal distribution to create a number of
“offspring” samples. Mimicking natural selection, these offspring are judged
according to some fitness function, a parallel to the upper-level loss
function. The fittest offspring determine the update to the normal
distribution and thus “pass on” their good characteristics to the next
generation.

####  3.2.2 Model-based

Model-based search strategies assume a model (or prior) for the hyperparameter
space and use only loss function evaluations (no gradients). This section
discusses two common model-based strategies: Bayesian approaches and trust
region methods.

##### Bayesian Approaches

Bayesian approaches fit previous hyperparameter trials’ results to a model to
select the hyperparameters that appear most promising to evaluate next
[klein:17:fbh]. For example, a common model for the hyperparameters is the
Gaussian Process prior. Given a few hyperparameter and cost function points, a
Bayesian method involves the following steps.

  1. 1.

Find the mean and covariance functions for the Gaussian Process. The mean
function will generally interpolate the sampled points. The covariance
function is generally expressed as a kernel function, often using squared
exponential functions [frazier:2018:tutorialbayesianoptimization].

  2. 2.

Create an acquisition function. The acquisition function captures how
desirable it is to sample (“acquire”) a hyperparameter setting. Thus, it
should be large (desirable) for hyperparameter values that are predicted to
yield small loss function values or that have high enough uncertainty that
they may yield low losses. The design of the acquisition function thus trades-
off between exploring new areas of the hyperparameter landscape with high
uncertainty and a more locally focused exploitation of the current best
hyperparameter settings. See [frazier:2018:tutorialbayesianoptimization] for a
discussion of specific acquisition function designs.

  3. 3.

Maximize the acquisition function (typically designed to be easy to optimize)
to determine which hyperparameter point to sample next.

  4. 4.

Evaluate the loss function at the new hyperparameter candidate.

These steps repeat for a given amount of time or until convergence.

##### Trust Region Method

Another derivative-free optimization method that uses only loss function
evaluations is a trust-region method. This section describes the specific
trust-region method as presented in
[ehrhardt:2021:inexactderivativefreeoptimization] (see references therein for
previous, similar methods).

The derivative-free, trust-region method (TRM) [conn:2000:trustregionmethods]
is similar to Bayesian optimization in that it involves fitting an easier to
optimize function to the loss function of interest, $\ell$, and then
minimizing the easier, surrogate function (the “model”). Thus, TRM requires
only function evaluations, not gradients, to construct and then minimize the
model. However, unlike most Bayesian optimization-based approaches, TRM uses a
local (often quadratic) model for $\ell$ around the current iterate, rather
than a surrogate that fits all previous points. In taking a step based on this
local information, TRM is more similar to gradient-based approaches.

The trust-region captures how well the model matches the observed $\ell$
values and determines the maximum step at every iteration. The “goodness” of
the model is typically quantified as the ratio of the actual decrease in
$\ell$ (based on observed function evaluations) to the predicted decrease
(based on the model). If this ratio is relatively large (close to one), then
the model is a good approximation of $\ell$ and the trust-region grows for the
next iteration. If this ratio is close to zero, then the observed decrease is
much less than predicted and the trust-region shrinks.

Recall that evaluating $\ell$ is typically expensive in bilevel problems as
each upper-level function evaluation involves optimizing the lower-level cost.
Thus, even constructing the model for a TRM can be expensive. To mitigate this
computational complexity, [ehrhardt:2021:inexactderivativefreeoptimization]
incorporated a dynamic accuracy component, with the accuracy for the lower-
level cost initially set relatively loose (leading to rough estimates of
$\ell$) but increasing with the upper-level iterations (leading to refined
estimates of $\ell$ as the algorithm nears a stationary point). One can use
any optimization method for the lower-level cost;
[ehrhardt:2021:inexactderivativefreeoptimization] used a gradient method for
the lower-level optimization method with well-known convergence results to
facilitate establishing convergence and computational complexity results.

The upper-level loss function considered in
[ehrhardt:2021:inexactderivativefreeoptimization] is additively separable and
quadratic:

$\displaystyle\ell({\bm{\upgamma}})$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}_{j}({\bm{\upgamma}}))=\frac{1}{J}\sum_{j=1}^{J}\underbrace{\left(r({\bm{\upgamma}}\,;{\bm{\hat{x}}}_{j}({\bm{\upgamma}}))\right)^{2}}_{\text{Equivalently,
}r_{j}({\bm{\upgamma}})^{2}},$

where $r$ is typically
${\bm{\hat{x}}}_{j}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}_{j}$. (Although we
define the sum to be over the number of training samples, this expression
easily generalizes to include a regularization term on $\bm{\upgamma}$ by
defining an additional $r_{J+1}$ term.)

Given a current value for the hyperparameters, ${\bm{\upgamma}}^{*}$, the TRM
models the local upper-level loss function by creating a linear model such
that

$r_{j}({\bm{\upgamma}}^{*}+\bm{\delta})\approx
m_{j}(\bm{\delta})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=r_{j}({\bm{\upgamma}}^{*})+{\bm{g}}_{j}^{\prime}\bm{\delta},$

where ${\bm{g}}_{j}\in\mathbb{F}^{R}$ approximates $\nabla
r_{j}({\bm{\upgamma}}^{*})$. Then, the overall model for the upper-level
problem is quadratic:

$\displaystyle\ell({\bm{\upgamma}}^{*}+\bm{\delta})$
$\displaystyle\approx\sum_{j}m^{2}_{j}(\bm{\delta})=\ell({\bm{\upgamma}}^{*})+\frac{1}{J}\sum_{j}2r_{j}({\bm{\upgamma}}^{*}){\bm{g}}_{j}^{\prime}\bm{\delta}+({\bm{g}}_{j}^{\prime}\bm{\delta})^{2}.$
(3.4)

One can estimate the gradients, ${\bm{g}}_{j}$, by interpolating a set of $R$
samples (recall ${\bm{\upgamma}}\in\mathbb{F}^{R}$) of the upper-level loss
function. This process involves choosing a set of interpolating points,
$\\{{\bm{\upgamma}}^{*}+\bm{\delta}^{(1)},\ldots,{\bm{\upgamma}}^{*}+\bm{\delta}^{(R)}\\}$,
(approximately) evaluating $r$ at each one, then solving

$\displaystyle\underbrace{\begin{bmatrix}r_{j}({\bm{\upgamma}}^{*}+\bm{\delta}^{(1)})-r_{j}({\bm{\upgamma}}^{*})\\\
\vdots\\\
r_{j}({\bm{\upgamma}}^{*}+\bm{\delta}^{(R)})-r_{j}({\bm{\upgamma}}^{*})\end{bmatrix}}_{R}=\underbrace{\begin{bmatrix}\left(\bm{\delta}^{(1)}\right)^{\prime}\\\
\vdots\\\ \left(\bm{\delta}^{(R)}\right)^{\prime}\end{bmatrix}}_{R\times
R}{\bm{g}}_{j}$

for $j\in[1\ldots J]$.

After forming the quadratic model for the upper-level loss, the TRM minimizes
the model (3.4) within some trust region, which is a simple convex-constrained
quadratic problem After computing $\hat{\bm{\delta}}$, the TRM accepts the
step and updates the hyperparameters
(${\bm{\upgamma}}^{(i+1)}={\bm{\upgamma}}^{(i)}+\hat{\bm{\delta}}$) if the
actual reduction (based on the estimated loss function values) to predicted
reduction (based on the quadratic model of the loss function) ratio is large
enough. Otherwise, if the ratio is low, the update step is rejected and the
trust region shrinks.

While the TRM appears to involve $R$ upper-level function evaluations every
iteration to construct the gradient estimates, after an initialization, one
can generally reuse samples, gradually replacing old samples with the samples
at new hyperparameter iterates. Ref.
[ehrhardt:2021:inexactderivativefreeoptimization] discusses requirements on
the interpolation set to guarantee a good geometry and conditions for re-
setting the interpolation sample if the model is not sufficiently accurate.

A main result from [ehrhardt:2021:inexactderivativefreeoptimization] is a
bound on the number of iterations to reach an $\epsilon$-optimal point
(defined as
$\min_{u}\|\nabla_{\bm{\upgamma}}\ell({{\bm{\upgamma}}^{(u)}})\|<\epsilon$,
where $u$ indexes the upper-level iterates). The bound derivation assumes (i)
$\Upphi$ is differentiable in $\bm{x}$, (ii) $\Upphi$ is $\mu$-strongly
convex, i.e.,
$\Upphi({\bm{x}})-\frac{\mu}{2}\left\lVert{\bm{x}}\right\rVert^{2}$ is convex,
(iii) the derivative of $\Upphi$ is Lipschitz continuous, and (iv) the first
and second derivative of the lower-level cost with respect to $\bm{x}$ exist
and are continuous. These requirements are satisfied by the example filter
learning problem (Ex), when $\bm{A}$ has full column rank, and more generally
when there are certain constraints on the hyperparameters. The iteration bound
is a function of the following:

  * •

the tolerance $\epsilon$,

  * •

the trust region parameters (parameters that control the increase and decrease
in trust region size based on the actual to predicted reduction, the starting
trust region size, and the minimum possible trust region size),

  * •

the initialization for $\bm{\upgamma}$, and

  * •

the maximum possible error between the gradient of the upper-level loss
function and the gradient of the model for the upper-level loss within a trust
region (when the gradient of $\ell$ is Lipschitz continuous, this bound is the
corresponding Lipschitz constant).

The number of iterations required to reach such an $\epsilon$-optimal point is
$\mathcal{O}\left(\frac{1}{\epsilon^{2}}\right)$
[ehrhardt:2021:inexactderivativefreeoptimization] and the number of required
upper-level loss function evaluations depends more than linearly on $N$
[roberts:2021:inexactdfobilevel]. That growth with $N$ impedes its use in
problems with many hyperparameters, though new techniques such as
[cartis:2021:scalablesubspacemethods] may be able to decrease or remove the
dependency.

###  3.3 Summary

Turning from the discussion of the lower-level problem in Section 2, this
section concentrated on the other two aspects of bilevel problems: the upper-
level loss function and the optimization strategy.

The loss function defines what a “good” hyperparameter is, typically using a
metric of image quality to compare ${\bm{\hat{x}}}({\bm{\upgamma}})$ to a
clean, training image, ${\bm{x}}^{\text{true}}$. MSE is the most common upper-
level loss function. It is well known from statistical estimation that the
estimator that minimizes MSE is the conditional mean,
${\bm{\hat{x}}}({\bm{y}})=\mathbb{E}\left[{\bm{x}}|{\bm{y}}\right]$. Thus, if
MSE is the true metric of interest, then lower-level problems should be
designed to try to approximate $\mathbb{E}\left[{\bm{x}}|{\bm{y}}\right]$
closely. Yet lower-level formulations in most bilevel papers are not described
as conditional mean estimators or approximations thereof. Section 3.1
discussed many other full reference and no reference options, including ones
motivated by human judgements of perceptual quality, from the image quality
assessment literature; Section 6.2 gives examples of bilevel methods that use
some of these other loss functions.

The second half of this section concentrated on model-free and model-based
hyperparameter search strategies. The grid search, CMA-ES, and trust region
methods described above all scale at least linearly with the number of
hyperparameters. Similarly, Bayesian optimization is best-suited for small
hyperparameter dimensions; [frazier:2018:tutorialbayesianoptimization]
suggests it is typically used for problems with 20 or fewer hyperparameters.

The remainder of this review considers gradient-based strategies for
hyperparameter optimization. The main benefit of gradient-based methods is
that they can scale to the large number of hyperparameters that are commonly
used in machine learning applications. Correspondingly, the main drawback of a
gradient-based method over the methods discussed in this section is the
implementation complexity and the differentiability requirement. Sections 4
and 5 discuss multiple options for gradient-based methods.

##  Chapter 4 Gradient Based Bilevel Methodology: The Groundwork

When the lower-level optimization problem (LL) has a closed-form solution,
${\bm{\hat{x}}}$, one can substitute that solution into the upper-level loss
function (UL). In this case, the bilevel problem is equivalent to a single
level problem and one can use classic single-level optimization methods to
minimize the upper-level loss. (See [kunisch:2013:bileveloptimizationapproach]
for analysis and discussion of some simple bilevel problems with closed-form
solutions for ${\bm{\hat{x}}}$.) This review focuses on the more typical
bilevel problems that lack a closed-form solution for ${\bm{\hat{x}}}$.

Although there are a wide variety of optimization methods for this challenging
category of bilevel problems, many methods are built on gradient descent of
the upper-level loss. The primary challenge with gradient-based methods is
that the gradient of the upper-level function depends on a variable that is
itself the solution to an optimization problem involving the hyperparameters
of interest. This section describes two main approaches for overcoming this
challenge. The first approach uses the fact that the gradient of the lower-
level cost function is zero at the minimizer to compute an exact gradient at
the exact minimizer. The second approach uses knowledge of the update scheme
for the lower-level cost function to calculate the exact gradient for an
approximation to the minimizer after a specific number of lower-level
optimization steps.

With this (approximation of the) gradient of the lower-level optimization
variable with respect to the hyperparameters, one can compute the gradient of
the upper-level loss function with respect to the hyperparameters,
$\bm{\upgamma}$. Section 5 uses the building blocks from this section to
explain various bilevel methods based on this gradient.

###  4.1 Set-up

Recall from Section 1.2 that a generic bilevel problem is

$\displaystyle\operatorname*{argmin}_{\bm{\upgamma}}\;$
$\displaystyle\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))\text{
where }$
$\displaystyle{\bm{\hat{x}}}({\bm{\upgamma}})=\operatorname*{argmin}_{{\bm{x}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}}).$
(4.1)

The way (4.1) expresses the lower-level problem implies a unique minimizer,
${\bm{\hat{x}}}$. This is satisfied when $\Upphi$ is a strictly convex
function of ${\bm{x}}$, which is typically the case when we use a strictly
convex regularizer. However, the following derivations also hold in the case
of non-unique or local minimizers. For more discussion,
[holler:2018:bilevelapproachparameter] defines optimistic and pessimistic
versions of the bilevel problem for the case of multiple lower-level
solutions.

For simplicity, hereafter we focus on the case $\mathbb{F}=\mathbb{R}$. Using
the chain rule, the gradient of the upper-level loss function with respect to
the hyperparameters is

$\displaystyle\nabla\ell({\bm{\upgamma}})$
$\displaystyle=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))+\left(\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})\right)^{\prime}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}})),$
(4.2)

where on the right hand side $\nabla_{\bm{\upgamma}}$ and $\nabla_{\bm{x}}$
denote partial derivatives w.r.t. the first and second arguments of
$\ell({\bm{\upgamma}}\,;\,{\bm{x}})$, respectively. We typically select the
loss function such that it is easy to compute these partials. For example, if
$\ell$ is the MSE training loss, i.e.,
$\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))=\frac{1}{2}\left\lVert{\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}\right\rVert^{2}_{2}$,
then

$\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))=0\text{
and
}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))={\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}.$

The following sections survey methods to find
$\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})$, which is the more
challenging piece in (4.2).

###  4.2 Minimizer Approach

The first approach finds
$\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})$ by assuming the
gradient of $\Upphi$ at the minimizer is zero. There are two ways to arrive at
the final expression: the implicit function theorem (IFT) perspective (as in
[samuel:2009:learningoptimizedmap,
gould:2016:differentiatingparameterizedargmin]) and the Lagrangian/KKT
transformation perspective (as in [chen:2014:insightsanalysisoperator,
holler:2018:bilevelapproachparameter]). This section presents both
perspectives in sequence. The end of the section summarizes the required
assumptions and discusses computational complexity and memory requirements.

The first step in both perspectives is to assume we have computed
${\bm{\hat{x}}}({\bm{\upgamma}})$ and that the lower-level problem 4.1 is
unconstrained (e.g., no non-negativity or box constraints). Therefore, the
gradient of $\Upphi$ with respect to $\bm{x}$ and evaluated at
${\bm{\hat{x}}}$ must be zero:

$\displaystyle\nabla_{\bm{x}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})\bigg{\rvert}_{{\bm{x}}={\bm{\hat{x}}}({\bm{\upgamma}})}=\nabla_{\bm{x}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})=\bm{0}.$
(4.3)

After this point, the two perspectives diverge.

####  4.2.1 Implicit Function Theorem Perspective

In the IFT perspective, we apply the IFT (_cf_. [fessler:96:mav]) to define a
function $h$ such that
${\bm{\hat{x}}}({\bm{\upgamma}})=h({\bm{y}},{\bm{\upgamma}})$. If we could
write $h$ explicitly, then the bilevel problem could be converted to an
equivalent single-level. However, per the IFT, we do not need to define $h$,
we only state that such an $h$ exists. Combining this definition with (4.3)
yields

$\displaystyle\bm{0}$
$\displaystyle=\nabla_{{\bm{x}}}\Upphi(h({\bm{y}},{\bm{\upgamma}});{\bm{\upgamma}}).$
(4.4)

Using the chain rule, we differentiate both sides of (4.4) with respect to
$\bm{\upgamma}$. The $\bm{I}$ in the equation below follows from the chain
rule because $\nabla_{\bm{\upgamma}}{\bm{\upgamma}}=\bm{I}$. We then rearrange
terms to solve for the desired quantity, noting that
$\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})=\nabla_{\bm{\upgamma}}h({\bm{y}},{\bm{\upgamma}})$.
Thus, evaluating all terms at ${\bm{\hat{x}}}$ leads to a gradient expression
for the lower-level problem:

$\displaystyle 0=$
$\displaystyle\nabla_{{\bm{x}}{\bm{x}}}\Upphi(h({\bm{y}},{\bm{\upgamma}});{\bm{\upgamma}})\nabla_{\bm{\upgamma}}h({\bm{y}},{\bm{\upgamma}})+\bm{I}\cdot\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi(h({\bm{y}},{\bm{\upgamma}});{\bm{\upgamma}})$
$\displaystyle\nabla_{\bm{\upgamma}}h({\bm{y}},{\bm{\upgamma}})=$
$\displaystyle-[\nabla_{{\bm{x}}{\bm{x}}}\Upphi(h({\bm{y}},{\bm{\upgamma}});{\bm{\upgamma}})]^{-1}\cdot\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi(h({\bm{y}},{\bm{\upgamma}});{\bm{\upgamma}})$
$\displaystyle\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})=$
$\displaystyle-[\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}})]^{-1}\cdot\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}}).$
(4.5)

When $\Upphi$ is strictly convex, the Hessian of $\Upphi$ is positive definite
and $\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}})$ is
invertible.

Substituting (4.5) into (4.2) yields the following expression for the gradient
of the upper-level loss function with respect to $\bm{\upgamma}$:

$\displaystyle\nabla\ell({\bm{\upgamma}})$
$\displaystyle=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}({\bm{\upgamma}}))-\left(\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}})\right)^{\prime}{\left(\nabla_{{\bm{x}}{\bm{x}}}\Upphi\left({\bm{\hat{x}}}\,;{\bm{\upgamma}}\right)\right)^{\text{-}1}}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}).$

Ref. [gould:2016:differentiatingparameterizedargmin] provides example simple
bilevel problems that verify that the IFT gradient agrees with the analytic
gradient obtained by substituting the closed-form solution of the lower-level
problem into the upper-level problem and differentiating directly.

####  4.2.2 KKT Conditions

In the Lagrangian perspective, (4.3) is treated as a constraint on the upper-
level problem, creating a single-level problem with $N$ equality
constraints111 The KKT transformation (4.6) relates bilevel optimization to
mathematical programs with equilibrium constraints (MPEC). Some authors use
the terms interchangeably; see [dempe:2020:bileveloptimizationadvances, Ch.
12]. :

$\displaystyle\operatorname*{argmin}_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}})\text{
subject to }\nabla_{\bm{x}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})={\bm{0}}_{N}.$
(4.6)

The corresponding Lagrangian is

$\displaystyle L({\bm{x}},{\bm{\upgamma}},{\bm{\nu}})$
$\displaystyle=\ell({\bm{\upgamma}}\,;{\bm{x}})+{\bm{\nu}}^{T}\nabla_{\bm{x}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})$

where ${\bm{\nu}}\in\mathbb{F}^{N}$ is a vector of Lagrange multipliers
associated with the $N$ equality constraints in (4.6). The Lagrange
reformulation is generally well-posed because many bilevel problems, such as
(Ex), satisfy the linear independence constraint qualification (LICQ)
[dempe:2003:annotatedbibliographybilevel, scholtes:2001:howstringentlinear].
LICQ requires that the matrix of derivatives of the constraint have full row
rank [scholtes:2001:howstringentlinear], i.e.

$\text{rank}\left(\begin{bmatrix}\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})&\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})\end{bmatrix}\right)=N.$

Strict convexity of $\Upphi({\bm{x}}\,;{\bm{\upgamma}})$ is therefore a
sufficient condition for LICQ to hold. (Note the similarity to the IFT
perspective, where strict convexity is sufficient for the Hessian to be
invertible.)

The first KKT condition states that, at the optimal point, the gradient of the
Lagrangian with respect to $\bm{x}$ must be $\bm{0}$. We can use this fact to
solve for the optimal Lagrangian multiplier, $\hat{{\bm{\nu}}}$:

$\displaystyle\nabla_{\bm{x}}L({\bm{\hat{x}}},{\bm{\upgamma}},\hat{{\bm{\nu}}})$
$\displaystyle=\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})+\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})\hat{{\bm{\nu}}}=\bm{0}$
$\displaystyle\hat{{\bm{\nu}}}$
$\displaystyle=\text{-}{\left(\nabla_{{\bm{x}}{\bm{x}}}\Upphi\left({\bm{\hat{x}}}\,;{\bm{\upgamma}}\right)\right)^{\text{-}1}}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}).$

Substituting the expression for $\hat{{\bm{\nu}}}$ into the gradient of the
Lagrangian with respect to $\bm{\upgamma}$ yields

$\displaystyle\nabla_{\bm{\upgamma}}L({\bm{\hat{x}}},{\bm{\upgamma}},\hat{{\bm{\nu}}})$
$\displaystyle=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})+\left(\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})\right)^{\prime}\hat{{\bm{\nu}}}$
$\displaystyle=$
$\displaystyle\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})-\left(\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})\right)^{\prime}{\left(\nabla_{{\bm{x}}{\bm{x}}}\Upphi\left({\bm{\hat{x}}}\,;{\bm{\upgamma}}\right)\right)^{\text{-}1}}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}),$

which is equivalent to (4.8).

Ref. [holler:2018:bilevelapproachparameter] generalized the Lagrangian
approach to the case where the forward model is defined only implicitly, e.g.,
as the solution to a differential equation. The authors write the lower-level
problem as

${\bm{\hat{x}}}=\operatorname*{argmin}_{{\bm{x}}}\min_{\tilde{{\bm{y}}}}\left\lVert{\bm{y}}-\tilde{{\bm{y}}}\right\rVert^{2}_{2}+R({\bm{x}})\text{
s.t. }e(\tilde{{\bm{y}}},{\bm{x}})=0,$ (4.7)

where the constraint function, $e$, incorporates the implicit system model.
For example, when the forward model is linear (${\bm{A}}{\bm{x}}$), taking
$e(\tilde{{\bm{y}}},{\bm{x}})=\left\lVert{\bm{A}}{\bm{x}}-\tilde{{\bm{y}}}\right\rVert^{2}_{2}$
shows the equivalence of the approach here to the one in
[holler:2018:bilevelapproachparameter].

####  4.2.3 Summary of Minimizer Approach

In summary, the upper-level gradient expression for the minimizer approach
(i.e., when one “exactly” minimizes the lower-level cost function) is

$\displaystyle\nabla\ell({\bm{\upgamma}})$
$\displaystyle=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})-\left(\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}})\right)^{\prime}{\left(\nabla_{{\bm{x}}{\bm{x}}}\Upphi\left({\bm{\hat{x}}}\,;{\bm{\upgamma}}\right)\right)^{\text{-}1}}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}}).$
(4.8)

Thus, for a given loss function and cost function, calculating the gradient of
the upper-level loss function (with respect to $\bm{\upgamma}$) requires the
following components all evaluated at ${\bm{x}}={\bm{\hat{x}}}$:
$\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}})\in\mathbb{F}^{R}$,
$\nabla_{{\bm{x}}{\bm{\upgamma}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})\in\mathbb{F}^{N\times
R}$,
$\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})\in\mathbb{F}^{N\times
N}$, and $\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}})\in\mathbb{F}^{N}$.

Continuing the specific example of learning filter coefficients and tuning
parameters (Ex), the components are:

$\displaystyle\nabla_{{\bm{x}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$
$\displaystyle={\bm{A}}^{\prime}({\bm{A}}{\bm{x}}-{\bm{y}})+e^{\beta_{0}}\sum_{k=1}^{K}e^{\beta_{k}}\tilde{{\bm{c}}}_{k}\circledast{{\dot{\phi}}}.({\bm{c}}_{k}\circledast{\bm{x}};\epsilon)$
$\displaystyle\nabla_{{\bm{x}}\beta_{k}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$
$\displaystyle=e^{\beta_{0}+\beta_{k}}\tilde{{\bm{c}}}_{k}\circledast{{\dot{\phi}}}.({\bm{c}}_{k}\circledast{\bm{\hat{x}}})$
$\displaystyle\nabla_{{\bm{x}}c_{k,s}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$
$\displaystyle=e^{\beta_{0}+\beta_{k}}\left({{\dot{\phi}}}.(({\bm{c}}_{k}\circledast{\bm{\hat{x}}})^{\langle{\bm{s}}\rangle})+\tilde{{\bm{c}}}_{k}\circledast\left({\ddot{\phi}}.({\bm{c}}_{k}\circledast{\bm{\hat{x}}})\odot{\bm{\hat{x}}}^{\langle{\bm{s}}\rangle}\right)\right)$
$\displaystyle\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$
$\displaystyle={\bm{A}}^{\prime}{\bm{A}}+e^{\beta_{0}}\sum_{k}e^{\beta_{k}}{\bm{C}}_{k}^{\prime}\text{diag}({\ddot{\phi}}.({\bm{c}}_{k}\circledast{\bm{\hat{x}}})){\bm{C}}_{k}$
$\displaystyle\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}})$
$\displaystyle=0$
$\displaystyle\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})$
$\displaystyle={\bm{\hat{x}}}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}.$ (4.9)

Here, the notation ${\bm{x}}^{\langle{\bm{i}}\rangle}$ means circularly
shifting the vector $\bm{x}$ by $\bm{i}$ elements, and $c_{k,{\bm{s}}}$
denotes the ${\bm{s}}$th element of the $k$th filter ${\bm{c}}_{k}$, where
${\bm{s}}$ is a tuple that indexes each dimension of ${\bm{c}}_{k}$. Appendix
A.1 gives examples of using the ${\bm{x}}^{\langle{\bm{i}}\rangle}$ notation
and derives $\nabla_{c_{k,s}}\left(\tilde{{\bm{c}}}_{k}\circledast
f.({\bm{c}}_{k}\circledast{\bm{x}})\right)$, which is the key step to
expressing
$\nabla_{{\bm{x}}c_{k,{\bm{s}}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$. The
other gradients easily follow from
$\nabla_{{\bm{x}}}\Upphi({\bm{\hat{x}}}\,;{\bm{\upgamma}})$ using standard
gradient tools for matrix expressions [petersen:2012:matrixcookbook].

The minimizer approach to finding $\nabla\ell({\bm{\upgamma}})$ uses the
following assumptions:

  1. 1.

Both the upper and lower optimization problems have no inequality constraints.

  2. 2.

${\bm{\hat{x}}}$ is the minimizer to the lower-level cost function, not an
approximation of the minimizer. This constraint ensures that (4.3) holds.

  3. 3.

The cost function $\Upphi$ is twice-differentiable in $\bm{x}$ and
differentiable with respect to $\bm{x}$ and $\bm{\upgamma}$. This condition
disqualifies the 0-norm and 1-norm as choices for the function $\phi$.

  4. 4.

The Hessian of the lower-level cost function,
$\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{x}}\,;{\bm{\upgamma}})$, is invertible;
this is guaranteed when $\Upphi$ is strictly convex.

The first condition technically excludes applications like CT imaging, where
the image is typically constrained to be non-negative. However, non-negativity
constraints are rarely required when good regularizers are used, so the
resulting non-constrained image can still be useful in practice
[fessler:96:mav].

The second constraint is often the most challenging since the lower-level
problem typically uses an iterative algorithm that runs for a certain number
of iterations or until a given convergence criteria is met. As previously
noted, if there were a closed-form solution for ${\bm{\hat{x}}}$, then we
would not have needed to use the IFT or Lagrangian to find the partial
derivative of ${\bm{\hat{x}}}$ with respect to $\bm{\upgamma}$. Since one
usually does not reach the exact minimizer, the calculated gradient will have
some error in it, depending on how close the final iterate is to the true
minimizer ${\bm{\hat{x}}}$. Thus, the practical application of this method is
more accurately called Approximate Implicit Differentiation (AID)
[ji:2021:bileveloptimizationconvergence,
grazzi:2020:iterationcomplexityhypergradient]. Section 4.5 further discusses
gradient accuracy.

####  4.2.4 Computational Costs

The largest cost in computing the gradient of the upper-level loss using (4.8)
is often finding (an approximation of) ${\bm{\hat{x}}}$. However, this cost is
difficult to quantify, as the IFT approach is agnostic to the lower-level
optimization methodology. To compare the bilevel gradient methods, we will
later assume the cost is comparable to the gradient descent calculations used
in the unrolled approach (described in the following section). However, this
is an over-estimation of the cost, as the IFT approach is not constrained to
smooth lower-level updates, and one can use optimization methods with, e.g.,
warm starts and restarts to reduce this cost.

When the lower-level problem satisfies the assumptions above, and assuming one
has already found ${\bm{\hat{x}}}$, a straight-forward approach to computing
the gradient (4.8) would be dominated by the $\mathcal{O}\left(N^{3}\right)$
operations required to compute the Hessian’s inverse. For many problems, $N$
is large, and that matrix inversion is infeasible due to computation or memory
requirements. Instead, as described in
[do:2007:efficientmultiplehyperparameter], one can use a conjugate gradient
(CG) method to compute the matrix-vector product

${\left(\nabla_{{\bm{x}}{\bm{x}}}\Upphi\left({\bm{\hat{x}}}\,;{\bm{\upgamma}}\right)\right)^{\text{-}1}}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{\hat{x}}})$
(4.10)

because the Hessian is symmetric and positive definite (see assumption #4 in
the previous section). For a generic $\bm{A}$, each CG iteration requires
multiplying the Hessian by a vector, which is $\mathcal{O}\left(N^{2}\right)$.

CG takes $N$ iterations to converge fully (ignoring finite numerical
precision), so the final complexity is still $\mathcal{O}\left(N^{3}\right)$
in general. However, the Hessian often has a special structure that simplifies
computing the matrix-vector product. Consider the running example of learning
filters per (Ex). The Hessian, as given in (4.9), multiplied with any vector
${\bm{v}}\in\mathbb{F}^{N}$ is

$\displaystyle\nabla_{{\bm{x}}{\bm{x}}}\Upphi({\bm{\hat{x}}};{\bm{\upgamma}},{\bm{y}})\cdot{\bm{v}}=$
$\displaystyle\underbrace{{\bm{A}}^{\prime}({\bm{A}}{\bm{v}})}_{2N^{2}}+$
$\displaystyle
e^{\beta_{0}}\sum_{k}e^{\beta_{k}}\underbrace{{\bm{C}}_{k}^{\prime}\cdot}_{NS}\overbrace{\text{diag}({\ddot{\phi}}.(\underbrace{{\bm{c}}_{k}\circledast{\bm{\hat{x}}}}_{NS}))\cdot}^{N}\underbrace{({\bm{C}}_{k}{\bm{v}})}_{NS}.$
(4.11)

The annotations show the multiplications required for each component, where we
used the simplifying assumption that the number of measurements matches the
number of unknowns ($M=N$).

As written, (4.11) does not made any assumptions on $\bm{A}$, so the first
term is still computationally expensive. If $\bm{A}$ is the identity matrix
(as in denoising), the $N^{2}$ term could instead be zero cost. If
${\bm{A}}^{\prime}{\bm{A}}$ is circulant, e.g., if $\bm{A}$ is a MRI sampling
matrix that can be written in terms of a discrete Fourier transform, then the
cost is $N\operatorname{log}\left(N\right)$. More generally, the computational
cost for one (of $N$) iterations of CG is
$\mathcal{O}\left(c_{{\bm{A}}}N\right)$ where $c_{{\bm{A}}}\in[0,N]$ is some
constant dependent on the structure of $\bm{A}$.

For the second addend in (4.11), we assume that $S\ll N$, so direct
convolution is most efficient and the matrix-vector product requires
$\mathcal{O}\left(NS\right)$ multiplies. When the filters are relatively
large, one can use Fourier transforms for the filtering, and the cost is
$\mathcal{O}\left(N\text{log}(N)\right)$. The final cost of the Hessian-vector
product for (Ex) is $\mathcal{O}\left(c_{{\bm{A}}}N+RN\right)$. This cost
includes a multiplication by $K$ to account for the sum over all filters,
which simplifies since $SK$ is222 The full parameter dimension includes the
filters and tuning parameters, so $R=S(K+1)+1$.  $\mathcal{O}\left(R\right)$.

###  4.3 Translation to a Single Level

Before discussing the other widely used approach to calculating the gradient
of the upper-level loss, we summarize a specialized approach for 1-norm
regularizers. Like the minimizer approach described above, this approach
assumes we have computed an (almost) exact minimizer of the lower-level cost
function. It writes the minimizer as an (almost everywhere) differentiable
function in terms of that ${\bm{\hat{x}}}$, then substitutes this expression
into the upper-level loss to create a single-level optimization problem that
is suitable for one hyperparameter update step.

Ref. [sprechmann:2013:supervisedsparseanalysis] proposed the translation to a
single-level approach to solve a bilevel problem with both synthesis and
analysis operators. mccann:2020:supervisedlearningsparsitypromoting recently
presented a version specific to analysis operators. Here, we summarize the
approach in [mccann:2020:supervisedlearningsparsitypromoting] that avoids an
auxiliary variable from splitting the one-norm sparsity penalty into two terms
in [sprechmann:2013:supervisedsparseanalysis].

The specific bilevel problem considered in
[mccann:2020:supervisedlearningsparsitypromoting] is:

$\displaystyle\operatorname*{argmin}_{{\bm{\upgamma}}}\sum_{j}\dfrac{1}{2}\|{\bm{\hat{x}}}_{j}({\bm{\upgamma}})-{\bm{x}}^{\text{true}}_{j}\|^{2}_{2}$
$\displaystyle\hat{{\bm{x}}}_{j}({\bm{\upgamma}})=$
$\displaystyle\operatorname*{argmin}_{{\bm{x}}}\dfrac{1}{2}\left\lVert{\bm{A}}{\bm{x}}-{\bm{y}}_{j}\right\rVert^{2}_{2}+\beta\left\lVert{\bm{\Omega}}_{\bm{\upgamma}}{\bm{x}}\right\rVert_{1},$
(4.12)

where ${\bm{\Omega}}_{\bm{\upgamma}}$ is a matrix constructed based on
$\bm{\upgamma}$, e.g., as a convolutional matrix as in our running example.
(In the results, [mccann:2020:supervisedlearningsparsitypromoting]
parameterizes ${\bm{\Omega}}_{\bm{\upgamma}}$ to have $3\times 3$
convolutional filters as the learnable parameters.) We write $\bm{\Omega}$ in
the following discussion to simplify notation. Eqn. (4.12) slightly
generalizes [mccann:2020:supervisedlearningsparsitypromoting] by considering
an inverse problem, rather than just denoising where ${\bm{A}}=\bm{I}$.

As in the minimizer approach, the first step is to compute
${\bm{\hat{x}}}({\bm{\upgamma}})$ for the current guess of $\bm{\upgamma}$,
e.g., using a fast proximal gradient method [kim:18:aro] or ADMM (as used in
[mccann:2020:supervisedlearningsparsitypromoting]). One can then re-write the
1-norm in (4.12) based on the sign pattern of
${\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})$ since

$\left\lVert{\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})\right\rVert_{1}=\sum_{i\in\mathcal{I}_{+}({\bm{\upgamma}})}[{\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})]_{i}-\sum_{i\in\mathcal{I}_{-}({\bm{\upgamma}})}[{\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})]_{i},$

where $\mathcal{I}_{+}({\bm{\upgamma}})$ and
$\mathcal{I}_{-}({\bm{\upgamma}})$ denote the set of indices where
${\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})$ is positive and negative,
respectively. Ref. [mccann:2020:supervisedlearningsparsitypromoting] defines a
diagonal sign matrix,
${\bm{S}}({\bm{\upgamma}})=\text{diag}(\mathrm{sign}\left({\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})\right))$,
having positive and negative diagonal elements at the appropriate indices. The
lower-level problem (4.12) is thus equivalent to

$\displaystyle\hat{{\bm{x}}}_{j}({\bm{\upgamma}})=\operatorname*{argmin}_{{\bm{x}}}\dfrac{1}{2}\left\lVert{\bm{A}}{\bm{x}}-{\bm{y}}_{j}\right\rVert^{2}_{2}+\beta{\bm{1}}^{\prime}{\bm{S}}({\bm{\upgamma}}){\bm{\Omega}}{\bm{x}},\text{
s.t. }[{\bm{\Omega}}{\bm{x}}]_{\mathcal{I}_{0}({\bm{\upgamma}})}={\bm{0}},$
(4.13)

where $\mathcal{I}_{0}({\bm{\upgamma}})$ denotes the set of indices where
$[{\bm{\Omega}}{\bm{\hat{x}}}({\bm{\upgamma}})]_{i}=0$.

The solution to the rewritten optimization problem (4.13) is the same as the
original version (4.12) with the 1-norm. (The rewriting process requires
${\bm{\hat{x}}}({\bm{\upgamma}})$, so one cannot use the equivalence to solve
the lower-level problem.) However, the rewritten problem has a simple
structure: it is a quadratic cost function with a linear equality constraint.
Therefore, using KKT conditions, if ${\bm{A}}^{\prime}{\bm{A}}$ is positive
definite, then the minimizer of (4.13) is the solution to a linear system of
equations

$\begin{bmatrix}{\bm{A}}^{\prime}{\bm{A}}&{\bm{\Omega}}_{0}^{\prime}({\bm{\upgamma}})\\\
{\bm{\Omega}}_{0}({\bm{\upgamma}})&{\bm{0}}\end{bmatrix}\begin{bmatrix}{\bm{\hat{x}}}({\bm{\upgamma}})\\\
{\bm{\nu}}\end{bmatrix}=\begin{bmatrix}{\bm{A}}^{\prime}{\bm{y}}-\beta{\bm{\Omega}}^{\prime}{\bm{S}}^{\prime}({\bm{\upgamma}}){\bm{1}}\\\
\bm{0}\end{bmatrix},$

where ${\bm{\Omega}}_{0}$ denotes the rows of $\bm{\Omega}$ corresponding to
the indices $\mathcal{I}_{0}({\bm{\upgamma}})$ and $\bm{\nu}$ is a vector of
Lagrange multipliers. See [boyd:04, Sect. 10.1.1] for discussion of existence
and uniqueness of solutions to this linear system.

Because of the discontinuity of the sign function in (4.13),
${\bm{\hat{x}}}({\bm{\upgamma}})$ is not everywhere differentiable in
$\bm{\upgamma}$. However, [mccann:2020:supervisedlearningsparsitypromoting]
states that ${\bm{\hat{x}}}({\bm{\upgamma}})$ is differentiable everywhere
except a set of measure zero when ${\bm{A}}=\bm{I}$ and when
${\bm{\Omega}}_{0}$ has full row rank. Then
[mccann:2020:supervisedlearningsparsitypromoting] uses matrix calculus to
derive a gradient expression for ${\bm{\hat{x}}}({\bm{\upgamma}})$ that is
suitable for standard automatic differentiation software tools.

In summary, the translation to a single level approach involves computing
${\bm{\hat{x}}}$, creating a closed-form expression for ${\bm{\hat{x}}}$, and
then differentiating the closed-form expression to compute the desired
gradient, $\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})$. In terms of
computation, the approach requires optimizing the lower-level cost and then
two calls to a CG method as part of calculating the gradient, making the cost
similar to the minimizer approach.

###  4.4 Unrolled Approaches

Another popular approach to finding
$\nabla_{\bm{\upgamma}}{\bm{\hat{x}}}({\bm{\upgamma}})$ is to assume that the
lower-level cost function is approximately minimized by applying $T$
iterations of some (sub)differentiable optimization algorithm, where we write
the update step at iteration $t\in[1\ldots T]$ as

${\bm{x}}^{(t)}=\Psi({\bm{x}}^{(t-1)}\,;{\bm{\upgamma}}),$

for some mapping $\Psi:\mathbb{F}^{N}\mapsto\mathbb{F}^{N}$ that should have
the fixed-point property
$\Psi({\bm{\hat{x}}}({\bm{\upgamma}})\,;{\bm{\upgamma}})={\bm{\hat{x}}}({\bm{\upgamma}})$.
For example, GD has
$\Psi({\bm{x}}\,;{\bm{\upgamma}})={\bm{x}}-\alpha\nabla\Upphi({\bm{x}}\,;{\bm{\upgamma}})$
for some step size $\alpha$. (Algorithms with momentum depend on earlier
iterates like ${\bm{x}}^{(t-2)}$ as well, but we omit that dependence to
simplify notation.) In contrast to the two approaches described above, this
“unrolled” approach no longer assumes the solution to the lower-level problem
is an exact minimizer. Instead, the unrolled approach reformulates the bilevel
problem (LL) as

$\displaystyle\operatorname*{argmin}_{\bm{\upgamma}}$
$\displaystyle\underbrace{\ell\left({\bm{\upgamma}}\,;\,{\bm{x}}^{(T)}({\bm{\upgamma}})\right)}_{\ell({\bm{\upgamma}})}\text{
s.t. }$ (4.14)
$\displaystyle{\bm{x}}^{(t)}({\bm{\upgamma}})=\Psi({\bm{x}}^{(t-1)}\,;{\bm{\upgamma}}),\quad\forall
t\in[1\ldots T],$

where ${\bm{x}}^{(0)}$ is an initialization, e.g., the pseudo-inverse of
$\bm{A}$ times $\bm{y}$. One can then take the (sub)gradient of a finite
number $T$ of iterations of $\Psi$, hoping that ${\bm{x}}^{T}$ approximately
minimizes the lower-level function, $\Upphi$.

The chain rule for derivatives is the foundation of the unrolled method. The
gradient of interest, $\nabla\ell({\bm{\upgamma}})$, depends on the gradient
of the optimization algorithm step with respect to $\bm{x}$ and
$\bm{\upgamma}$. For readability, define the following gradient matrices for
the $t$th unrolled iteration

$\displaystyle{\bm{H}}_{t}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\nabla_{{\bm{x}}}\Psi\left({\bm{x}}^{(t-1)}\,;{\bm{\upgamma}}\right)\in\mathbb{F}^{N\times
N}\text{ and }{\bm{J}}_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\nabla_{{\bm{\upgamma}}}\Psi\left({\bm{x}}^{(t-1)}\,;{\bm{\upgamma}}\right)\in\mathbb{F}^{N\times
R},$

for $t\in[1,T]$. We use these letters because, when using gradient descent as
the optimization algorithm, $\nabla_{\bm{x}}\Psi({\bm{x}}\,;{\bm{\upgamma}})$
is closely related to the Hessian of $\Upphi$ and
$\nabla_{\bm{\upgamma}}\Psi({\bm{x}}\,;{\bm{\upgamma}})$ is the Jacobian of
the gradient. Thus, when $\Psi$ corresponds to GD, an unrolled approach
involves computing the same quantities as required by the IFT approach (4.8).

By the chain rule, the gradient of (4.14) is

$\displaystyle\nabla\ell({\bm{\upgamma}})=$
$\displaystyle\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+\left(\sum_{t=1}^{T}\left({\bm{H}}_{T}\cdots{\bm{H}}_{t+1}\right){\bm{J}}_{t}\right)^{\prime}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})\in\mathbb{F}^{R}.$
(4.15)

One can derive this gradient expression using a reverse or forward
perspective, with parallels to back-propagation through time and real-time
recurrent learning respectively [franceschi:2017:forwardreversegradientbased].
We present the update here in terms of $\bm{x}$ for simplicity, but the idea
easily extends to updates in terms of a state vector that allows one to
include momentum terms, weights, and other accessory variables in
$\bm{\upgamma}$ [franceschi:2017:forwardreversegradientbased].

Unlike the minimizer approach, where the goal is to run the lower-level
optimization until (close to) convergence, most unrolling methods set $T$ in
advance. This set number of lower-level iterations mimics the set depth of
CNNs and allows a precise estimate of how much computational effort each
lower-level optimization takes. However, if $T$ is relatively small, the
lower-level optimization may not near convergence, and the unrolled method is
only loosely tied to the original bilevel optimization problem. To maintain a
strong connection to the bilevel problem while avoiding setting $T$ larger
than necessary for convergence, [antil:2020:bileveloptimizationdeep] used a
convergence criterion to determine the number of $\Psi$ iterations rather than
pre-specifying a number of iterations. Section 7.1 further discusses the
relation between bilevel problems and unrolling methods in the broader
literature [monga:21:aui].

####  4.4.1 Computational and Memory Complexity

To compare the forward and reverse approaches to gradient computation for
unrolled methods, we introduce notation for an ordered product of matrices. We
indicate the arrangement of the multiplications by the set endpoints,
$s\in[s_{1}\leftrightarrow s_{2}]$ with the left endpoint, $s_{1}$,
corresponding to the index for the left-most matrix in the product and the
right endpoint, $s_{2}$, corresponding to the right-most matrix. Thus, for any
sequence of square matrices $\\{{\bm{A}}\\}_{i}$:

$\displaystyle\prod_{s\in\left[t\leftrightarrow
T\right]}{\bm{A}}_{s}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}={\bm{A}}_{t}{\bm{A}}_{t+1}\cdots{\bm{A}}_{T}=\left({\bm{A}}_{T}^{\prime}{\bm{A}}_{T-1}^{\prime}\cdots{\bm{A}}_{t}^{\prime}\right)^{\prime}=\left(\prod_{s\in\left[T\leftrightarrow
t\right]}{\bm{A}}_{s}^{\prime}\right)^{\prime}.$

The above double arrow notation does not indicate order of operations. In the
following notation the arrow direction does not affect the product result
(ignoring finite precision effects), but rather signifies the direction
(order) of calculation:

$\displaystyle\prod_{s\in\left[T\leftarrow t\right]}{\bm{A}}_{s}$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}={\bm{A}}_{T}\left({\bm{A}}_{T-1}\cdots\left({\bm{A}}_{t+1}\left({\bm{A}}_{t}\right)\right)\right)$
$\displaystyle\prod_{s\in\left[T\rightarrow t\right]}{\bm{A}}_{s}$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\left(\left(\left({\bm{A}}_{T}{\bm{A}}_{T-1}\right)\cdots\right){\bm{A}}_{t+1}\right){\bm{A}}_{t}.$

We use a similar arrow notation to denote the order that terms are computed
for sums; as above, the order is only important for computational
considerations and does not affect the final result.

Figure 4.1:  Reverse mode computation of the unrolled gradient from (4.1). The
first gradient computation requires ${\bm{x}}^{(T)}$, so all computations
occur after the lower-level optimization algorithm is complete. The final
gradient is
$\nabla\ell({\bm{\upgamma}})=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+{\bm{r}}$.

Using this notation, the reverse gradient calculation of (4.15) is

$\displaystyle\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+\sum_{t\in[T\rightarrow
1]}{\bm{J}}_{t}{}^{\prime}\left(\prod_{s\in[(t+1)\leftarrow
T]}{\bm{H}}_{s}^{\prime}\right)\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)}).$
(4.16)

This expression requires $\prod_{s\in[(T+1)\leftarrow
T]}{\bm{H}}_{s}^{\prime}=\bm{I}$, because ${\bm{H}}_{T+1}$ is not defined. For
example, for $T=3$, we have

$\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(3)})+\underbrace{{\bm{J}}_{3}^{\prime}(\bm{I}){\bm{g}}}_{t=3}+\underbrace{{\bm{J}}_{2}^{\prime}\left({\bm{H}}_{3}^{\prime}\right){\bm{g}}}_{t=2}+\underbrace{{\bm{J}}_{1}^{\prime}\left({\bm{H}}_{2}^{\prime}{\bm{H}}_{3}^{\prime}\right){\bm{g}}}_{t=1},$

where $\bm{g}$ is shorthand for
$\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})$ here. This version is
called reverse as all computations (arrows) begin at the end, $T$.

The primary benefit of the reverse mode comes from the ability to group
$\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})$ with the right-most
${\bm{H}}_{T}$, such that all products are matrix-vector products, as seen in
Fig. 4.1 Further, one can save the matrix-vector products for use during the
next iteration and avoid duplicating the computation. Continuing the example
for $T=3$, we have

$\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(3)})+\underbrace{{\bm{J}}_{3}^{\prime}(\bm{I}){\bm{g}}}_{t=3}+\underbrace{{\bm{J}}_{2}^{\prime}(\overbrace{{\bm{H}}_{3}^{\prime}{\bm{g}}}^{\bm{\Delta}})}_{t=2}+\underbrace{{\bm{J}}_{1}^{\prime}({\bm{H}}_{2}^{\prime}\overbrace{\left({\bm{H}}_{3}^{\prime}{\bm{g}}\right)}^{\bm{\Delta}})}_{t=1},$

where one only needs to compute $\bm{\Delta}$ once. This ability to rearrange
the parenthesis to compute matrix-vector products greatly decreases the
computational requirement compared to matrix-matrix products. Excluding the
costs of the optimization algorithm steps and forming the ${\bm{H}}_{s}$ and
${\bm{J}}_{t}$ matrices (these costs will be the same in the forward mode
computation), reverse mode requires $\mathcal{O}\left(T\right)$ Hessian-vector
multiplies and $\mathcal{O}\left(TNR\right)$ additional multiplies. The trade-
off is that reverse mode requires storing all $T$ iterates, ${\bm{x}}^{(t)}$,
so that one can compute the corresponding Hessians and Jacobians from them as
needed, and thus has a memory complexity $\mathcal{O}\left(TN\right)$.

Figure 4.2: Forward mode computation of the unrolled gradient from (4.17). The
intermediate computation matrix, $\bm{Z}$, is initialized to zero
(${\bm{Z}}_{0}=\bm{0}$) then updated every iteration. The final gradient is
$\nabla\ell({\bm{\upgamma}})=\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+{\bm{Z}}_{T}^{\prime}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})$.

The forward mode calculation of (4.15), depicted in Fig. 4.2, has all
computations (arrows) starting at the earlier iterate:

$\displaystyle\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+\left(\sum_{t\in[1\rightarrow
T]}\left(\prod_{s\in[T\leftarrow(t+1)]}{\bm{H}}_{s}\right){\bm{J}}_{t}\right)^{\prime}\nabla_{\bm{x}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)}).$
(4.17)

As before, ${\bm{H}}_{T+1}$ is not defined, so we take
$\prod_{s\in[T\leftarrow(T+1)]}{\bm{H}}_{s}=\bm{I}$. For example, for $T=3$ we
have

$\displaystyle\nabla_{\bm{\upgamma}}$
$\displaystyle\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+\left(\underbrace{\left(({\bm{H}}_{3}{\bm{H}}_{2}){\bm{J}}_{1}\right)^{\prime}}_{t=1}+\underbrace{\left(({\bm{H}}_{3}){\bm{J}}_{2}\right)^{\prime}}_{t=2}+\underbrace{\left((\bm{I}){\bm{J}}_{3}\right)^{\prime}}_{t=3}\right){\bm{g}}.$

How the forward mode avoids storing $\bm{x}$ iterates is evident after
rearranging the parenthesis to avoid duplicate calculations, as illustrated in
Fig. 4.2. Continuing the example for $T=3$, we have

$\nabla_{\bm{\upgamma}}\ell({\bm{\upgamma}}\,;{\bm{x}}^{(T)})+\left[\underbrace{{\bm{H}}_{3}\left(\overbrace{{\bm{H}}_{2}\underbrace{\left({\bm{H}}_{1}\cdot\bm{0}+{\bm{J}}_{1}\right)}_{{\bm{Z}}_{1}}+{\bm{J}}_{2}}^{{\bm{Z}}_{2}}\right)+{\bm{J}}_{3}}_{{\bm{Z}}_{3}}\right]^{\prime}{\bm{g}},$

where
${\bm{Z}}_{s}={\bm{H}}_{s}{\bm{Z}}_{s-1}+{\bm{J}}_{s}\in\mathbb{F}^{N\times
R}$ stores the intermediate calculations. The above formula also illustrates
why ${\bm{H}}_{1}$ is not needed in (4.14);
$\nabla_{\bm{\upgamma}}{\bm{x}}^{(0)}=\bm{0}$ is the last element from
applying the chain rule.

There is no way to rearrange the terms in the forward mode formula to achieve
matrix-vector products (while preserving the computation order). Therefore,
the computation requirement is much higher at $\mathcal{O}\left(TR\right)$
Hessian-vector multiplications. The corresponding benefit of the forward mode
method is that it does not require storing iterates, thus decreasing (in the
common case when $T>R$) the memory requirement to $\mathcal{O}\left(NR\right)$
for storing the intermediate matrix ${\bm{Z}}_{s}$ during calculation.

As with the IFT/Lagrangian approach, the computational complexity of the
unrolled approach is lower than the generic bound when we consider the
specific example of learning convolutional filters according to (Ex).
Nevertheless, the general comparison that reverse mode takes more memory but
less computation holds true. See Tab. 4.1 for a full comparison of the
computational and memory complexities.

Neural network training methods commonly use the reverse method,
backpropagation, to compute (4.15) because it requires less computation than
forward mode. However, in some problems with large dimensions, storing all the
iterates is infeasible. A strategy to trade-off the memory and computation
requirements is checkpointing, which stores $\bm{x}$ every few iterations.
Checkpointing is an active research area; see
[dauvergne:2006:dataflowequationscheckpointing] for an overview.

###  4.5 Summary