∎ aFEMclipboard 11institutetext: A. Chambolle 22institutetext: CEREMADE, CNRS
& Université Paris Dauphine, PSL Research University, Paris
22email<EMAIL_ADDRESS>ORCID: 0000-0002-9465-4659
33institutetext: R. Tovey 44institutetext: MOKAPLAN, INRIA Paris, Paris
44email<EMAIL_ADDRESS>ORCID: 0000-0001-5411-2268

# “FISTA” in Banach spaces with adaptive discretisations

Antonin Chambolle  Robert Tovey

###### Abstract

FISTA is a popular convex optimisation algorithm which is known to converge at
an optimal rate whenever a minimiser is contained in a suitable Hilbert space.
We propose a modified algorithm where each iteration is performed in a subset
which is allowed to change at every iteration. Sufficient conditions are
provided for guaranteed convergence, although at a reduced rate depending on
the conditioning of the specific problem. These conditions have a natural
interpretation when a minimiser exists in an underlying Banach space. Typical
examples are L1-penalised reconstructions where we provide detailed
theoretical and numerical analysis.

###### Keywords:

Convex optimization Multiscale Multigrid Sparsity Lasso

††journal: Computational Optimization and Applications

##  1 Introduction

The Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) was proposed by
Beck and Teboulle Beck2009  as an extension of Nesterov’s fast gradient method
Nesterov2004  and is now a very popular algorithm for minimising the sum of
two convex functions. We write this as the problem of computing

$\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\qquad\text{such
that}\qquad\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\coloneqq\operatorname{f}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])+\operatorname{g}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]),$
(1)

for a Hilbert space $\mathds{H}$ where
$\operatorname{f}\colon\mathds{H}\to\mathds{R}$ is a convex differentiable
function with $L$-Lipschitz gradient and
$\operatorname{g}\colon\mathds{H}\to\overline{\mathds{R}}$ is a “simple”
convex function, whose “proximity operator” is easy to compute. Throughout
this work we assume that $\operatorname{E}$ is bounded below so that the
infimum is finite. The iterates of the FISTA algorithm will be denoted
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}\in\mathds{H}$. If, moreover the
infimum is achieved, it has been shown that
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
converges at the optimum rate of $n^{-2}$ Beck2009 , and later (after a small
modification) the convergence of the iterates was also shown in a general
Hilbert space setting Chambolle2015 . Many further works have gone on to
demonstrate faster practical convergence rates for slightly modified variants
of FISTA Tao2016 ; Liang2017 ; Alamo2019 .

In this work we address the case where the minimiser possibly fails to exist
or lies in a larger space where $\mathds{H}$ is dense. There is much overlap
between the techniques used in this work and those used in the literature of
inexact optimisation, however, our interpretation is relatively novel. In
particular, we emphasise the infinite-dimensional setting where errors come
from “discretisation”, rather than random or decaying errors in $\mathds{H}$,
which enables two new perspectives:

  * •

Analytically, we prove new rates of convergence for FISTA when the minimum
energy is not achieved (at least not in $\mathds{H}$). The exact rate can be
computed by quantifying coercivity and regularity properties of
$\operatorname{E}$. If there isn’t a minimiser in $\mathds{H}$, then this rate
is strictly slower than $n^{-2}$.

  * •

Numerically, we allow the optimisation domain to change on every iteration.
This enables us to understand how FISTA behaves with adaptive discretisations.
Adaptive finite-element methods are known to improve the efficiency of, for
example, approximating the solutions of PDEs. Our analytical results show how
to combine such tools with FISTA without reducing the guaranteed rate of
convergence, and our numerical results confirm much improved time and computer
memory efficiency in the Lasso example (Section 6).

All the examples in this work, discussed from Section 5 onward, consider
$\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}\operatorname{\;s.t.\;}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])<\infty\\}$
to be contained in some ambient Banach space $\mathds{U}$. The idea is that
FISTA provides a minimising sequence in $\mathds{H}\cap\mathds{U}$, but
further properties like rate of convergence (of $\operatorname{E}$ or the
iterates) must come from the topology of $\mathds{U}$. It will not be
necessary for $\mathds{H}\hookrightarrow\mathds{U}$ to be a continuous
embedding, nor in fact the full inclusion $\mathds{H}\subset\mathds{U}$.

Some other works for FISTA-like algorithms include Jiang2012 ; Villa2013 . Of
particular note, our stability estimate for FISTA in Theorem 4.1 is very
similar to (Schmidt2011, , Prop 2) and (Aujol2015, , Prop 3.3). This is then
used to analyse the convergence properties in our more general Banach space
setting, but where all sources of inexactness come from subspace
approximations. The ideas in Parpas2017  are similar although in application
to the proximal gradient method with an additional smoothing on the functional
$\operatorname{g}$. The permitted refinement steps are also more broad in our
work. Very recent work in Yu2021  proposes a “Multilevel FISTA” algorithm
which allows similar coarse-to-fine refinement strategies, although only a
finite number. We also allow for non-uniform refinement with a posteriori
strategies.

###  1.1 Outline

This work is organised as follows. Section 2 defines notation and the generic
form of our proposed refining FISTA algorithm, Algorithm 1. The main
theoretical contribution of this work is the convergence analysis of Algorithm
1 which is split into two parts: first we outline the proof structure in
Section 3, then we state the specific results in the case of FISTA in Section
4. The main results are Theorems 4.2/4.3 which extend the convergence of FISTA
to cases with un-attained minima with uniform/adaptively chosen subspaces
$\mathds{U}^{n}$ respectively.

Section 5 presents some general results for the application of Algorithm 1 in
Banach spaces and Section 6 gives a much more detailed discussion of adaptive
refinement for Lasso minimisation. In particular, we describe how to choose
efficient refining discretisations to approximate
$\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$,
estimate the convergence of $\operatorname{E}$, and identify the support of
the minimiser. The numerical results in Section 7 demonstrate these techniques
in four different models demonstrating the apparent sharpness of our
convergence rates and the computational efficiency of adaptive
discretisations.

##  2 Definitions and notation

We consider optimisation of (1) over a Hilbert space
$(\mathds{H},\left\langle\cdot,\cdot\right\rangle,{\left\lVert\cdot\right\rVert})$.
In the more analytical section (Sections 3 and 4) it will be more convenient
to use the translated energy

$\operatorname{E}_{0}\colon\mathds{H}\to\mathds{R},\qquad\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\coloneqq\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])-\inf_{\widetilde{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}\in\mathds{H}}\operatorname{E}(\widetilde{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]})$
(2)

so that
$\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=0$,
although access to this function is not assumed for numerical examples.

The proposed generalised FISTA algorithm is stated in Algorithm 1 for an
arbitrary choice of closed convex subsets $\mathds{U}^{n}\subset\mathds{H}$
for $n\in\mathds{N}$. The only difference from standard FISTA is that on
iteration $n$, all computations are performed in the subset $\mathds{U}^{n}$.
If $\mathds{U}^{n}=\mathds{H}$, then we recover the original algorithm. More
generally, the idea is that $\mathds{U}^{n}$ are “growing”, for example
$\mathds{U}^{n}\subset\mathds{U}^{n+1}$, but this assumption is not necessary
in most of the results.

Without loss of generality we will assume $L=1$, i.e. $\nabla\operatorname{f}$
is 1-Lipschitz. To get the general statement of any of the results which
follow, replace $\operatorname{E}$ with $\frac{\operatorname{E}}{L}$. In
particular,

${\left\lVert\nabla\operatorname{f}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])-\nabla\operatorname{f}(\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1])\right\rVert}\leq{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\right\rVert}$
(3)

for all
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0],\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\mathds{H}$
and $\operatorname{g}$ is called “simple” if it is proper, convex, weakly
lower-semicontinuous, and

$\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}}\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\right\rVert}^{2}+\operatorname{g}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
(4)

is exactly computable for all
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\mathds{H}$ and all
$\widetilde{\mathds{U}}\in\\{\mathds{U}^{n}\\}_{n=0}^{\infty}$. Closed subsets
of $\mathds{H}$ are locally weakly compact, therefore this argmin is always
non-empty.

One defining property of the FISTA algorithm is an appropriate choice of
inertia, dictated by $t_{n}$. In particular, we will say that
$(t_{n})_{n=0}^{\infty}$ is a _FISTA stepsize_ if

$t_{0}=1,\qquad t_{n}\geq 1,\qquad\text{and}\qquad\rho_{n}\coloneqq
t_{n}^{2}-t_{n+1}^{2}+t_{n+1}\geq 0\qquad\text{ for all }n=0,1,\ldots.$ (5)

The precise constants associated to a given rate are given in the statements
of the theorems but, for convenience, are otherwise omitted from the text. For
sequences $(a_{n})_{n=0}^{\infty}$,$(b_{n})_{n=0}^{\infty}$ we will use the
notation:

$\displaystyle a_{n}\lesssim b_{n}\qquad$ $\displaystyle\iff\qquad\exists
C,N>0\operatorname{\;s.t.\;}a_{n}\leq Cb_{n}\text{ for all }n>N,$
$\displaystyle a_{n}\simeq b_{n}\qquad$ $\displaystyle\iff\qquad a_{n}\lesssim
b_{n}\lesssim a_{n}.$

For $n\in\mathds{N}$ we use the abbreviation $[n]=\\{1,2,\ldots,n\\}$. When
the subdifferential of $\operatorname{E}$ is set-valued, we will use the
short-hand

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\coloneqq\inf_{\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
(6)

for any specified norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$.

Algorithm 1 Refining subset FISTA

1:Choose $(\mathds{U}^{n})_{n\in\mathds{N}}$,
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}\in\mathds{U}^{0}$ and some FISTA
stepsize choice $(t_{n})_{n\in\mathds{N}}$

2:$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{0}\leftarrow\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0},n\leftarrow
0$

3:repeat

4:
$\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}_{n}\leftarrow(1-\tfrac{1}{t_{n}})\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}+\tfrac{1}{t_{n}}\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$

5:
$\displaystyle\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n+1}\leftarrow\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n+1}}\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}_{n}+\nabla\operatorname{f}(\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}_{n})\right\rVert}^{2}+\operatorname{g}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
$\triangleright$ Only modification, $\mathds{U}^{n+1}\subset\mathds{U}$

6:
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n+1}\leftarrow(1-t_{n})\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}+t_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n+1}$

7: $n\leftarrow n+1$

8:until some stopping criterion is met

##  3 General proof recipe

In this section we give an intuitive outline of the full proof for convergence
of Algorithm 1 before giving formal theorems and proofs in the next section.
First we recall the classical FISTA convergence guarantee given by
(Chambolle2015, , Thm 3.1); if there exists
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$,
then

$t_{N}^{2}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})+\sum^{N-1}_{n=1}\rho_{n}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\tfrac{1}{2}{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{N}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}^{2}\leq\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}^{2}$
(7)

for any FISTA stepsize choice $t_{N}\simeq N$ such that $\rho_{n}\geq 0$.

##### Step 1: Quantifying the stability

The first step is to generalise (7) to account for the adapting subsets
$\mathds{U}^{n}$. In the notation of Algorithm 1, Theorem 4.1 shows that

$t_{N}^{2}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})+\sum_{n=1}^{N-1}\rho_{n}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\tfrac{1}{2}{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{N}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{N}\right\rVert}^{2}\leq\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{0}\right\rVert}^{2}+\tfrac{{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{N}\right\rVert}^{2}-{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{0}\right\rVert}^{2}}{2}+\sum^{N}_{n=1}t_{n}\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n})+\left\langle\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n-1},\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n-1}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\right\rangle$
(8)

for any $\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\in\mathds{U}^{n}$. The
similarities to (7) are clear. If $\mathds{U}^{n}=\mathds{H}$, then we can
choose
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$
and the two estimates agree. These extra terms in (8) quantify the robustness
to changing of discretisation.

##### Step 2: Quantifying the scaling properties

To show that the extra terms in (8) are small, we need to quantify the
approximation properties of $\mathds{U}^{n}$. The idea is that there is a
sequence $\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\in\mathds{U}^{n}$,
$n\in\mathds{N}$ such that
${\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\right\rVert}$ grows
slowly and $\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n})$
decreases quickly. To quantify this balance, we introduce a secondary sequence
$n_{0}<n_{1}<\ldots$ and constants
$a_{\operatorname{U}},a_{\operatorname{E}}\geq 1$ such that for each
$k\in\mathds{N}$

$n\leq
n_{k}\implies{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\right\rVert}\lesssim
a_{\operatorname{U}}^{k},\qquad n\geq
n_{k}\implies\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n})\lesssim
a_{\operatorname{E}}^{-k}.$ (9)

A canonical example would be
$\mathds{U}^{n}=\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}\operatorname{\;s.t.\;}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}\leq
a_{\operatorname{U}}^{k}\\}$ for $n\in[n_{k},n_{k+1})$, then
$a_{\operatorname{E}}$ reflects the smoothness of $\operatorname{E}_{0}$. The
choice of exponential scaling is introduced to improve stability of Algorithm
1. It is natural if we consider the $\mathds{U}^{n}$ to be the subspace of
functions discretised on a uniform mesh. If that mesh is sequentially refined,
then the resolution of the mesh will be of order $h^{k}$ after $k$ refinements
and for some $h<1$. The integer $n_{k}$ is then the time at which the mesh has
refined $k$ times. The trade-off between $a_{\operatorname{E}}$ and
$a_{\operatorname{U}}$ dictates the final convergence rate of the algorithm.
If $a_{\operatorname{U}}>1$, then we cannot guarantee the original $n^{-2}$
rate of convergence.

##### Step 3: Generalising the convergence bound

In this step we combine the FISTA stability estimate with the subset
approximation guarantees to provide a sharper estimate of stability with
respect to the parameters $a_{\operatorname{E}}$ and $a_{\operatorname{U}}$.
For example, if for each $k\in\mathds{N}$

$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n_{k}}\text{
for each }n=n_{k},n_{k}+1,\ldots,n_{k+1}-1,$

then many terms on the right-hand side of (8) telescope to 0. The result of
this is presented in Lemma 3. The key idea is that the stability error in (8)
has $K\ll N$ terms, rather than $N$.

##### Step 4: Sufficiently fast growth

In Step 3 we develop a convergence bound, now we wish to show that it is only
worse than the classical (7) by a constant factor. In particular, it is
equivalent to either run Algorithm 1 for $N$ iterations, or the classical
FISTA algorithm for $N$ iterations on the fixed subset $\mathds{U}^{N}$. The
estimate from (7) provides the estimate
$N^{2}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{N}\right\rVert}^{2}=O(a_{\operatorname{U}}^{2K})$
for $N\leq n_{K}$. Lemma 4 shows that Algorithm 1 can achieve the same order
of approximation, so long as $\mathds{U}^{n}$ grow sufficiently quickly (in
particular $n_{k}^{2}\lesssim
a_{\operatorname{E}}^{k}a_{\operatorname{U}}^{2k}$).

##### Step 5: Sufficiently slow growth

The result of Step 3 is sufficient to prove convergence, but not yet a rate.
If the subsets grow too quickly, then the influence of
${\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}\right\rVert}\to\infty$
will slow the rate of convergence. If $n_{k}$ is too large, then we overfit to
the discrete problem, but if $n_{k}$ is too small, then FISTA converges
slowly. Lemma 5 balances these two factors in an optimal way ($n_{k}^{2}\simeq
a_{\operatorname{E}}^{k}a_{\operatorname{U}}^{2k}$) for Algorithm 1 resulting
in a convergence rate of

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{a_{\operatorname{U}}^{2K}}{N^{2}}\lesssim\frac{N^{2\kappa}}{N^{2}}$

for all $N\in\mathds{N}$ and $\kappa=\frac{2\log a_{\operatorname{U}}}{\log
a_{\operatorname{E}}+2\log a_{\operatorname{U}}}\in[0,1)$. In particular, if
the minimum is attained in $\mathds{H}$, then we recover the classical rate
with $\kappa=0$.

##### Step 6: Adaptivity

Up to this point we have implicitly focused on the case where $\mathds{U}^{n}$
(and $n_{k}$) are chosen a priori. The main challenge for adaptive choice of
$\mathds{U}^{n}$ is to guarantee (9) from Step 3 using a posteriori estimates.
Combined with the partial telescoping requirement in Step 3, a natural choice
is
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1}$
for $n\in[n_{k},n_{k+1})$, i.e. the value of $n_{k}$ is chosen to be $n+1$
once the iterate $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ is observed.
Theorem 4.3 shows that a sufficient condition is

$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1},\qquad{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1}\right\rVert}\lesssim
a_{\operatorname{U}}^{k},\quad\text{and}\quad\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1})\lesssim
a_{\operatorname{E}}^{-k}.$

Convergence is most stable if the approximation spaces $\mathds{U}^{n}$
satisfy a monotone inclusion, breaking the monotonicity requires more care.
The only non-trivial property to verify is the energy gap
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$.
Lemma 6 proposes some sufficient conditions to guarantee the same overall rate
of convergence as in Step 3,

$\min_{n\leq
N}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\lesssim\frac{N^{2\kappa}}{N^{2}}$

for all $N\in\mathds{N}$, with the same $\kappa\in[0,1)$ from Step 3. The
penalty for accelerating the change of discretisation is a potential loss of
stability or monotonicity in
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$, although
this behaviour has not been seen in numerical experiments.

##  4 Proof of convergence

In this section we follow the recipe motivated in Section 3 to prove
convergence of two variants of Algorithm 1. Each of the main theorems and
lemmas will be stated with a sketch proof in this section. The details of the
proofs are either trivial or very technical and are therefore placed in
Section A to preserve the flow of the argument.

###  4.1 Computing the convergence bound

For Step 3 of Section 3 we look to replicate the classical bound of the form
in (7) for Algorithm 1. The proofs in this step follow the classical arguments
Beck2009 ; Chambolle2015  very closely. Throughout this section we consider a
sequence $(\mathds{U}^{n})_{n\in\mathds{N}}$ which generate the iterates
$(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})_{n\in\mathds{N}}$ in Algorithm 1
such that

$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{n}\in\mathds{U}^{n+1}\subset\mathds{H}\quad\text{where
}\mathds{U}^{n}\text{ is a closed, convex subset for all }n\in\mathds{N}.$
(10)

####  4.1.1 Single iterations

We first wish to understand a single iteration of Algorithm 1. This is done
through the following two lemmas.

######  Lemma 1 (equivalent to (Chambolle2015, , Lemma 3.1))

thm: descent lemma Suppose $\nabla\operatorname{f}$ is 1-Lipschitz, for any
$\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}\in\mathds{H}$ define

$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\coloneqq\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}+\nabla\operatorname{f}(\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]})\right\rVert}^{2}+\operatorname{g}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]).$

Then, for all $\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]\in\mathds{U}^{n}$, we have

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])+\tfrac{1}{2}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]\right\rVert}^{2}\leq\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2])+\tfrac{1}{2}{\left\lVert\overline{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]\right\rVert}^{2}.$

The proof is exactly the same as in Chambolle2015  on the subset
$\mathds{U}^{n}$. Applying Lemma 1 to the iterates from Algorithm 1 gives a
more explicit inequality.

######  Lemma 2 ((Chambolle2015, , (17)), (Beck2009, , Lemma 4.1))

thm: one step FISTA Let
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\in\mathds{U}^{n}$ be chosen
arbitrarily and
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$/$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$
be generated by Algorithm 1 for all $n\in\mathds{N}$. For all $n>0$, it holds
that

$\Copy{thm:eq:onestepFISTA}{t_{n}^{2}(\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}))-(t_{n}^{2}-t_{n})(\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n-1})-\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}))\leq\tfrac{1}{2}\left[{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n-1}\right\rVert}^{2}-{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}\right\rVert}^{2}\right]+\left\langle\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}-\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n-1},\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\right\rangle.}$
(11)

The proof is given in Theorem A.1 and is a result of the convexity of
$\operatorname{E}_{0}$ and $\mathds{U}_{n}$ for a well chosen
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]$ in Lemma 1.

####  4.1.2 Generic convergence bound

Lemma 2 gives us an understanding of a single iteration of Algorithm 1,
summing over $n$ then gives our generic convergence bound for any variant of
Algorithm 1.

######  Theorem 4.1 (analogous to (Chambolle2015, , Thm 3.2), (Beck2009, , Thm
4.1))

thm: mini FISTA convergence Fix a sequence of subsets
$(\mathds{U}^{n})_{n\in\mathds{N}}$ satisfying (10), arbitrary
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}\in\mathds{U}^{0}$, and FISTA
stepsize choice $(t_{n})_{n\in\mathds{N}}$. Let
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ and
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$ be generated by Algorithm 1, then,
for any choice of $\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\in\mathds{U}^{n}$
and $N\in\mathds{N}$ we have

$\Copy{thm:eq:miniFISTAconvergence}{t_{N}^{2}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})+\sum_{n=1}^{N-1}\rho_{n}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\frac{{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{N}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{N}\right\rVert}^{2}}{2}\leq\frac{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{0}\right\rVert}^{2}-{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{0}\right\rVert}^{2}+{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{N}\right\rVert}^{2}}{2}+\sum^{N}_{n=1}t_{n}\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n})+\left\langle\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n-1},\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n-1}-\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}\right\rangle.}$
(12)

The proof is given in Theorem A.2. This result is the key approximation for
showing convergence of FISTA with changing subsets. In the classical setting,
we have $\mathds{U}^{n}=\mathds{H}$,
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{0}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
and the extra terms on the right-hand side collapse to 0.

If there exists a minimiser
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$,
then the natural choice in (12) is
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$
for some projection $\mathsf{\Pi}_{n}\colon\mathds{H}\to\mathds{U}^{n}$,
however, there are simple counter-examples which give
$\operatorname{E}_{0}(\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})=\infty$
and so this inequality becomes useless. For example, if
$\operatorname{f}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])={\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{L^{2}([0,1])}^{2}$,
$\operatorname{g}$ is the indicator on the set
$\mathds{D}=\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in
L^{1}([0,1])\operatorname{\;s.t.\;}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0](x)\geq
x\\}$, and $\mathsf{\Pi}_{n}$ is the $L^{2}$ projection onto a set of
piecewise constant functions, then
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}=x\mapsto x$. On the other hand,
suppose one of the pixels of the discretisation is $[x_{0}-h,x_{0}+h]$, then

$\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\left(x_{0}+\tfrac{h}{2}\right)=\operatorname*{argmin}_{\IfEqCase{5}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k5]\in\mathds{R}}\int_{x_{0}-h}^{x_{0}+h}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(x)-\IfEqCase{5}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k5])^{2}\mathop{}\\!\mathrm{d}x=\operatorname*{argmin}_{\IfEqCase{5}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k5]\in\mathds{R}}\int_{x_{0}-h}^{x_{0}+h}(x-\IfEqCase{5}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k5])^{2}\mathop{}\\!\mathrm{d}x=x_{0}<x_{0}+\tfrac{h}{2}.$

In particular
$\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\notin\mathds{D}$
therefore
$\operatorname{E}_{0}(\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})=\infty$.
The choice
$\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}=\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
is much more robust and allows us to apply Algorithm 1 more broadly. The
penalty for this flexibility is a more complicated analysis; each time the
subset changes, because
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}\in\mathds{U}_{n}$, the system
receives a “shock” proportional to
${\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}\right\rVert}{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n}-\mathsf{\Pi}_{n}\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]_{n+1}\right\rVert}$.

###  4.2 Convergence bound with milestones

In standard FISTA, the right-hand side of (12) is a constant. The following
lemma minimises the growth of the “constant” as a function of $N$ by partially
telescoping the sum on the right-hand side. Before progressing to the content
of Step 3, we will first formalise the definition of the constants
$a_{\operatorname{U}}$ and $a_{\operatorname{E}}$ introduced in Step 3.

###### Definition 1

Fix $a_{\operatorname{U}},a_{\operatorname{E}}\geq 1$ and a sequence
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{H}$. We say
that $(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$
is an _$(a_{\operatorname{U}},a_{\operatorname{E}})$ -minimising sequence of
$\operatorname{E}$ _ if

${\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rVert}\lesssim
a_{\operatorname{U}}^{k}\qquad\text{and}\qquad\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})\lesssim
a_{\operatorname{E}}^{-k}$

for all $k\in\mathds{N}$.

In this section we will simply assume that such sequences exist and in Section
5 we will give some more general examples.

###### Lemma 3

thm: mini exponential FISTA convergence Let
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$,
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$ be generated by Algorithm 1 with
$(\mathds{U}^{n})_{n\in\mathds{N}}$ satisfying (10),
$(n_{k}\in\mathds{N})_{k\in\mathds{N}}$ be a monotone increasing sequence, and
choose

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1}$

for each $k\in\mathds{N}$. If such a sequence exists, then for all
$K\in\mathds{N}$, $n_{K}\leq N<n_{K+1}$ we have

$\Copy{thm:eq:miniexponentialFISTAconvergence}{t_{N}^{2}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})+\sum_{n=1}^{N-1}\rho_{n}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\frac{{\left\lVert\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{N}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{K}\right\rVert}^{2}}{2}\leq
C+\frac{{\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{K}\right\rVert}^{2}}{2}+\frac{(N+1)^{2}-n_{K}^{2}}{2}\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{K})\\\
+\sum_{k=1}^{K}\frac{n_{k}^{2}-n_{k-1}^{2}}{2}\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k-1})+\left\langle\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n_{k}-1},\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k-1}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rangle}$

thm:end: mini exponential FISTA convergence where
$C=\frac{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{0}\right\rVert}^{2}-{\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{0}\right\rVert}^{2}}{2}$.

The proof is given in Lemma 11. The introduction of $n_{k}$ has greatly
compressed the expression of Theorem 4.1. On the right-hand side, we now only
consider $\operatorname{E}_{0}$ evaluated on the sequence
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$ and there are $K$
elements to the sum rather than $N$.

###  4.3 Refinement without overfitting

The aim of Step 3 is to show that $n$ iterations of Algorithm 1 is no slower
(up to a constant factor) than $n$ iterations of classical FISTA on the space
$\mathds{U}^{n}$. In other words, we would like to ensure that

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\lesssim\frac{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rVert}^{2}}{n^{2}}$
(13)

uniformly for $n\in[n_{k},n_{k+1})$. If this condition is not satisfied, then
it indicates that computational effort has been wasted by a poor choice of
subsets. This can be interpreted as an overfitting to the discretisation of
$\operatorname{E}_{0}|_{\mathds{U}^{n}}$ rather than the desired function
$\operatorname{E}_{0}|_{\mathds{H}}$. Combining the assumptions given by
Definition 1 and the result of Lemma 3, the following lemma proves the
convergence of Algorithm 1 provided that the refinement times $n_{k}$ are
sufficiently small (i.e. $\mathds{U}^{n}$ grows sufficiently quickly).

###### Lemma 4

thm: sufficiently fast Suppose $\mathds{U}^{n},\
\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n},\
\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$ and $n_{k}$ satisfy the conditions
of Lemma 3 and
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$
forms an
$(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ with

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1}.$

If either:

  * •

$a_{\operatorname{U}}>1$ and $n_{k}^{2}\lesssim
a_{\operatorname{E}}^{k}a_{\operatorname{U}}^{2k}$,

  * •

or $a_{\operatorname{U}}=1$,
$\sum_{k=1}^{\infty}n_{k}^{2}a_{\operatorname{E}}^{-k}<\infty$, and
$\sum_{k=1}^{\infty}{\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k+1}\right\rVert}<\infty,$

then

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{a_{\operatorname{U}}^{2K}}{N^{2}}\qquad\text{for
all}\qquad n_{K}\leq N<n_{K+1}.$

The proof is given in Lemma 12. We make two observations of the optimality of
Lemma 4:

  * •

The convergence guarantee for $N\in[n_{K},n_{K+1})$ iterations of classical
FISTA in the space $\mathds{U}^{N}$ is

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{K}\right\rVert}^{2}}{N^{2}}+\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{N}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\lesssim\frac{a_{\operatorname{U}}^{2K}}{N^{2}}+a_{\operatorname{E}}^{-K}.$

This is equivalent to Lemma 4 after the assumptions on $n_{k}$.

  * •

If $\mathds{H}$ is finite dimensional, then the condition
$a_{\operatorname{U}}=1$ is almost trivially satisfied. Norms in finite
dimensions are equivalent and any discretisation can be achieved with a finite
number of refinements (i.e. the sums over $k$ are finite).

###  4.4 Convergence rate

In Lemma 4 we show that
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$ converges at
a rate depending on $k$ and $n$, so long as $k$ grows sufficiently quickly. On
the other hand, as $k$ grows, the rate becomes worse and so we need to also
put a lower limit on the growth of $n_{k}$. The following lemma completes Step
3 by computing the global convergence rate of
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$ when $k$
grows at the minimum rate which is consistent with Lemma 4.

As a special case, note that if $a_{\operatorname{U}}=1$ then Lemma 4 already
gives the optimal $O(N^{-2})$ convergence rate. This is in fact a special case
of that shown in (Aujol2015, , Prop 3.3). If the minimum is achieved in
$\mathds{H}$, then it is not possible to refine “too quickly” and the
following lemma is not needed.

###### Lemma 5

thm: sufficiently slow Suppose $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ and
$n_{k}$ are sequences satisfying

$\forall N\in[n_{K},n_{K+1}),\
\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{a_{\operatorname{U}}^{2K}}{N^{2}}\qquad\text{where}\qquad
n_{K}^{2}\gtrsim a_{\operatorname{E}}^{K}a_{\operatorname{U}}^{2K},$

then

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{1}{N^{2(1-\kappa)}}\qquad\text{
where }\qquad\kappa=\frac{\log a_{\operatorname{U}}^{2}}{\log
a_{\operatorname{E}}+\log a_{\operatorname{U}}^{2}}.$

The proof is given in Lemma 13.

####  4.4.1 FISTA convergence with a priori discretisation

We can summarise Lemmas 3 to 5 into a single theorem stating the convergence
guarantees when $\mathds{U}^{n}$ and $n_{k}$ are chosen a priori.

###### Theorem 4.2

thm: exponential FISTA convergence Let
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$ be
an $(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ and choose any $\mathds{U}^{n}$ satisfying (10) such that

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1}$

for all $k\in\mathds{N}$. Compute $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$
and $\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$ by Algorithm 1.

Suppose that either:

  * •

$a_{\operatorname{U}}>1$ and $n_{k}^{2}\simeq
a_{\operatorname{E}}^{k}a_{\operatorname{U}}^{2k}$, or

  * •

$a_{\operatorname{U}}=1$,
$\sum_{k=1}^{\infty}n_{k}^{2}a_{\operatorname{E}}^{-k}<\infty$ and
$\sum_{k=1}^{\infty}{\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k+1}\right\rVert}<\infty$,

then

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\lesssim\frac{1}{N^{2(1-\kappa)}}\qquad\text{
where }\qquad\kappa=\frac{\log a_{\operatorname{U}}^{2}}{\log
a_{\operatorname{E}}+\log a_{\operatorname{U}}^{2}}\qquad\text{uniformly for
}N\in\mathds{N}.$

Analytically, this theorem gives new rates of convergence for FISTA when the
minimiser is not achieved in $\mathds{H}$. Indeed for the original algorithm
($\mathds{U}^{n}=\mathds{H}$), if $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{0}=0$
for simplicity and
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$ is
_any_ $(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ exists, the result of Lemma 3 is

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})\leq\inf_{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]\in\mathds{H}}\frac{{\left\lVert\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]\right\rVert}^{2}+N^{2}\operatorname{E}_{0}(\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2])}{2t_{N}^{2}}\leq\min_{k\in\mathds{N}}\frac{{\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rVert}^{2}+N^{2}\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})}{2t_{N}^{2}}\lesssim\min_{k\in\mathds{N}}\frac{a_{\operatorname{U}}^{2k}+N^{2}a_{\operatorname{E}}^{-k}}{N^{2}}\lesssim
N^{-2(1-\kappa)}.$ (14)

In this sense, we could say that $\operatorname{E}_{0}$ converges at the rate
$N^{-2(1-\kappa)}$ if and only if such a sequence exists. Nothing is lost (or
gained) analytically by choosing $\mathds{U}_{n}\subsetneq\mathds{H}$.

Numerically, it is easy to implement the strategy of Theorem 4.2 and requires
very little knowledge of how to estimate
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$. So long as
$a_{\operatorname{U}}$ and $a_{\operatorname{E}}$ can be computed
analytically, one can choose
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$ implicitly to be the
discrete minimisers of some “uniform” discretisations (e.g.
$\mathds{U}^{n}=\\{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}\leq
k\\}$ or finite element spaces with uniform mesh) to achieve the stated
convergence rate.

####  4.4.2 FISTA convergence with adaptivity

There are two properties of the sequence $(\mathds{U}^{n})_{n\in\mathds{N}}$
which we may wish to decide adaptively: the refinement times $n_{k}$ and the
discretising spaces $\\{\mathds{U}^{n}\operatorname{\;s.t.\;}n_{k}\leq
n<n_{k+1}\\}$. We will refer to these as temporal and spatial adaptivity
respectively.

Lemma 4 gives a sufficient condition on $n_{k}$ for converging at the rate
$O(N^{2(\kappa-1)})$, but it is not necessary. Indeed for $n\leq n_{k}$ we
have

$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\geq\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=O(a_{\operatorname{E}}^{-k})=O(n^{2(\kappa-1)}),$

which suggests that to converge faster than $n^{2(\kappa-1)}$ requires
choosing smaller $n_{k}$. As an example, in Section 7.2 we will see Algorithm
1 can converge at a near-linear rate, although this is not possible without
adaptive refinement times. On the other hand, choice of spatial adaptivity has
no impact on rate but can impact computational efficiency. It will be
permitted to use greedy discretisation techniques so long as it is sufficient
to estimate $\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$
accurately.

Theorem 4.2 already allows for spatial adaptivity, so we focus on temporal
adaptivity. Lemma 4 suggests that a good refinement time strategy is to choose
$n_{k}$ to be the minimal integer such that
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1})\lesssim
a_{\operatorname{E}}^{-k}$. However, the value of $\operatorname{E}_{0}$ may
be hard to estimate and so we retain a “backstop” condition which guarantees
that convergence is no slower than the rate given by Theorem 4.2. In the non-
classical case of $a_{\operatorname{U}}>1$, we provide the following theorem.

###### Theorem 4.3

thm: stronger exponential FISTA convergence Let
$(\mathds{U}^{n}\subset\mathds{H})_{n\in\mathds{N}}$ be a sequence of subsets
satisfying (10), compute $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ and
$\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]_{n}$ by Algorithm 1. Suppose that there
exists a monotone increasing sequence $n_{k}\in\mathds{N}$ such that

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\coloneqq\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1}$

for all $k\in\mathds{N}$.

If $(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$
is an $(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ with $a_{\operatorname{U}}>1$ and $n_{k}^{2}\lesssim
a_{\operatorname{E}}^{k}a_{\operatorname{U}}^{2k}$, then

$\min_{n\leq
N}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})=\min_{n\leq
N}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\lesssim\frac{1}{N^{2(1-\kappa)}}\qquad\text{
where }\qquad\kappa=\frac{\log a_{\operatorname{U}}^{2}}{\log
a_{\operatorname{E}}+\log a_{\operatorname{U}}^{2}}$

uniformly for $N\in\mathds{N}$.

The proof is given in Theorem A.3. If we directly compare Theorems 4.2 and
4.3, both are a direct result of Lemma 4 assuming a specific choice of $n_{k}$
or $\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$ respectively. We
note that the convergence rate is the same in both theorems but the price for
better adaptivity (i.e. only an upper bound on $n_{k}$) is a slightly weaker
stability guarantee (now convergence of $\min_{n\leq
N}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$). In Theorem
4.2, as in the original FISTA algorithm, the sequence
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$ is not
monotone but the magnitude of oscillation is guaranteed to decay in time. This
behaviour is lost in Theorem 4.3. Although we do not prove it here, it can be
shown that the stronger condition

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{U}^{n_{k}}\cap\mathds{U}^{n_{k}+1}\cap\ldots\cap\mathds{U}^{n_{k+1}-1}\cap\ldots\cap\mathds{U}^{N}$
(15)

is sufficient to restore the stronger last-iterate guarantee on
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{N})$. Again,
monotonicity of $\mathds{U}^{n}$ corresponds with improved stability of
Algorithm 1.

To enable a more practical implementation of Theorem 4.3, the following lemma
describes several refinement strategies which provide sufficient condition for
$\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})\lesssim
a_{\operatorname{E}}^{-k}$.

###### Lemma 6

thm: practical refinement criteria Let
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$ be a
sequence in $\mathds{H}$ with
${\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rVert}\lesssim
a_{\operatorname{U}}^{k}$. Suppose
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\widetilde{\mathds{U}}^{k}\coloneqq\mathds{U}^{n_{k}}$
and denote
$\operatorname{E}_{0}(\widetilde{\mathds{U}}^{k})\coloneqq\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}^{k}}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$.
Any of the following conditions are sufficient to show that
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$ is an
$(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$:

  1. 1.

Small continuous gap refinement:
$\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
for all $k\in\mathds{N}$, some
$\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]>0$.

  2. 2.

Small discrete gap refinement:
$\operatorname{E}_{0}(\widetilde{\mathds{U}}^{k})\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
and
$\operatorname{E}_{0}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})-\operatorname{E}_{0}(\widetilde{\mathds{U}}^{k-1})\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
for all $k>0$, some
$\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]>0$.

Otherwise, suppose there exists a Banach space
$(\mathds{U},{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|})$
which contains each $\widetilde{\mathds{U}}^{k}$,
$\sup_{k\in\mathds{N}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}<\infty$,
and the sublevel sets of $\operatorname{E}$ are
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$-bounded.
With the subdifferential
$\partial\operatorname{E}\colon\mathds{U}\rightrightarrows\mathds{U}^{*}$, it
is also sufficient if either:

  1. 3.

Small continuous gradient refinement:
$\sup_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\inf_{\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\partial\operatorname{E}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})}\frac{|\left\langle\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1],\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rangle|}{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
for all $k\in\mathds{N}$, some
$\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]>0$.

  2. 4.

Small discrete gradient refinement:
$\operatorname{E}_{0}(\widetilde{\mathds{U}}^{k})\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
and
$\sup_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0],\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}\in\widetilde{\mathds{U}}^{k}}\inf_{\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\mathds{V}^{k}}\frac{|\left\langle
v,\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}\right\rangle|}{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}\leq\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]a_{\operatorname{E}}^{-k}$
for all $k\in\mathds{N}$, some
$\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]>0$,
where
$\mathds{V}^{k}\coloneqq\partial(\operatorname{E}|_{\widetilde{\mathds{U}}^{k}})(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})$.

The proof is given in Lemma 14. The refinement criteria described by Lemma 6
can be split into two groups. Cases (1) and (3) justify that any choice of
$\mathds{U}^{n_{k}}$ satisfies the required conditions, so long as
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{U}^{n_{k}}$.
In cases (2) and (4), $\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$
is sufficient to choose the refinement time $n_{k}$, but an apriori bound is
required on $\operatorname{E}_{0}(\widetilde{\mathds{U}}^{k})$. In these cases
one could, for example, choose $\widetilde{\mathds{U}}^{k}$ to be a uniform
discretisation with a priori estimates.

Another splitting of the criteria is into gap and gradient computations.
Typically, gradient norms (in (4) and (5)) should be easier to estimate than
function gaps because they only require local knowledge rather than global,
i.e. $\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$
rather than an estimate of
$\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$.
Implicitly, the global information comes from an extra condition on
$\operatorname{E}$ to assert that sublevel sets are bounded.

##  5 General examples

We consider the main use of Algorithm 1 to be where there exists a Banach
space
$(\mathds{U},{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|})$
such that
$\mathds{U}\supset\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}\operatorname{\;s.t.\;}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])<\infty\\}$
and

$\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$

for some $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\mathds{U}$. The cases
where $\mathds{H}$ has finite dimension or is separable are more
straightforward; if the total number of refinements is finite (i.e.
$\mathds{U}^{n}=\mathds{U}^{N}$ for all $n\geq N$, some $N\in\mathds{N}$),
then $a_{\operatorname{U}}=1$. This holds for most finite dimensional problems
as well as the countable example discussed in detail in Section 6. In this
section we give explicit computations of $a_{\operatorname{U}}$ and
$a_{\operatorname{E}}$ in the setting where $\mathds{H}=L^{2}(\Omega)$ for
some domain $\Omega\subset\mathds{R}^{d}$ and the subsets $\mathds{U}^{n}$
will be finite dimensional finite-element–like spaces, as defined below.

###### Definition 2

Suppose
${\left\lVert\cdot\right\rVert}_{q}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
(i.e. $\mathds{U}\subset L^{q}(\Omega)$) for some $q\in[1,\infty]$ and
connected, bounded, measurable domain $\Omega\subset\mathds{R}^{d}$. We say
that a collection $\mathds{M}$ is a _mesh_ if

$\bigcup_{\omega\in\mathds{M}}\omega\supset\Omega\qquad\text{and}\qquad|\omega\cap\omega^{\prime}|=0\qquad\text{for
all $\omega,\omega^{\prime}\in\mathds{M},\ \omega\neq\omega^{\prime}$.}$

Furthermore, we say a sequence of meshes $(\mathds{M}^{k})_{k\in\mathds{N}}$
is _consistent_ if there exists $\omega_{0}\subset\Omega$ such that

$\forall\omega\in\mathds{M}^{k}\quad\exists(\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]_{\omega},\vec{\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]}_{\omega})\in\mathds{R}^{d\times
d}\times\mathds{R}^{d}\quad\text{such
that}\quad\vec{x}\in\omega_{0}\iff\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]_{\omega}\vec{x}+\vec{\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]}_{\omega}\in\omega.$

Fix $h\in(0,1)$, linear subspaces
$\widetilde{\mathds{U}}^{k}\subset\mathds{H}$, and consistent meshes
$\mathds{M}^{k}$. We say that the sequence
$(\widetilde{\mathds{U}}^{k})_{k\in\mathds{N}}$ is an _$h$ -refining sequence
of finite element spaces_ if there exists
$c_{\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]}>0$
such that:

$\forall(\widetilde{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]},\omega)\in\widetilde{\mathds{U}}^{k}\times\mathds{M}^{k},\quad\operatorname{det}(\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]_{\omega})\geq
c_{\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]}h^{kd}\quad\text{and}\quad\exists\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}^{0}\quad\text{such
that}\quad\forall\vec{x}\in\omega_{0},\
\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0](\vec{x})=\widetilde{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}(\IfEqCase{0}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k0]_{\omega}\vec{x}+\vec{\IfEqCase{1}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k1]}_{\omega}).$

We say that $(\widetilde{\mathds{U}}^{k})_{k\in\mathds{N}}$ is _of order $p$_
if for any
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
there exists a sequence
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$ such
that

$\forall
k\in\mathds{N},\qquad\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\widetilde{\mathds{U}}^{k}\quad\text{and}\quad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\lesssim_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}}h^{kp}.$
(16)

We allow the implicit constant to have any dependence on
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$ so long as it is finite. For
example, in the case of Sobolev spaces we would expect an inequality of the
form
${\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}_{W^{0,2}}\lesssim
h^{kp}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}_{W^{p,2}}$
Strang1972 .

###### Remark 1

To clarify this definition with an example, suppose we wish to approximate
$L^{q}(\Omega)$ with piecewise linear finite elements with a triangulated
mesh. Then, $\omega_{0}\subset\Omega$ is a single triangle of diameter $O(h)$
and all meshes $\mathds{M}^{k}$ must be triangulations of $\Omega$ with cell
volumes scaling no faster than $O(h^{kd})$. The function
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]$ from the $h$-refining property is an
arbitrary linear element, so that each
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}^{k}$ is linear
on each $\omega\in\mathds{M}^{k}$, which leads to an order $p=2$ if
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in W^{1,2}(\Omega)$.

We note that any piecewise polynomial finite element (or spline) space can be
used to form a $h$-refining sequence of subspaces. Wavelets with a compactly
supported basis behave like a multi-resolution finite element space as there
is always overlap in the supports of basis vectors. Similarly, a Fourier basis
does satisfy the scaling properties, but each basis vector has global support.
Both of these exceptions are important and could be accounted for with further
analysis but we focus on the more standard finite element case. In order to
align these discretisation properties with the assumptions of Theorems 4.2 and
4.3, we make the following observation.

###### Lemma 7

Fix
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
and $p^{\prime},q^{\prime}>0$. If a sequence
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\mathds{H}$ satisfies

${\left\lVert\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right\rVert}\lesssim
h^{-kq^{\prime}}\quad\text{and}\quad\operatorname{E}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\lesssim
h^{kp^{\prime}},$

then $(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$
is an $(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ for $a_{\operatorname{U}}=h^{-q^{\prime}}$ and
$a_{\operatorname{E}}=h^{-p^{\prime}}$.

This is precisely rewriting the statement of Definition 1 into terms of
resolution $h$. The following theorem links $p$ and $q$ from Definition 2 with
$p^{\prime}$ and $q^{\prime}$ from Lemma 7.

###### Theorem 5.1

Suppose $\mathds{H}=L^{2}(\Omega)$ for some connected, bounded domain
$\Omega\subset\mathds{R}^{d}$ and
${\left\lVert\cdot\right\rVert}_{q}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
for some $q\in[1,\infty]$. For $p\geq 0$ and $h\in(0,1)$, if
$(\widetilde{\mathds{U}}^{k})_{k\in\mathds{N}}$ is an $h$-refining sequence of
finite element spaces of order $p$, then
$(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})_{k\in\mathds{N}}$ is
an $(a_{\operatorname{U}},a_{\operatorname{E}})$-minimising sequence of
$\operatorname{E}$ for

$\displaystyle a_{\operatorname{U}}$ $\displaystyle\leq\begin{cases}1&\text{if
}q\geq 2,\\\ \sqrt{h^{-d}}&q<2\text{ and
}\sup_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}^{0}}\frac{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{L^{\infty}(\omega_{0})}}{{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{L^{2}(\omega_{0})}}<\infty\end{cases},$
$\displaystyle a_{\operatorname{E}}$
$\displaystyle\geq\begin{cases}h^{-2p}&\text{if }\nabla\operatorname{E}\text{
is
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$-Lipschitz
at }\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*},\\\ h^{-p}\qquad&\text{if
}\operatorname{E}\text{ is
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$-Lipschitz
at }\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*},\\\
1&\text{otherwise.}\end{cases}$

The proof of this theorem is in Appendix B. Note that
$\sup_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\widetilde{\mathds{U}}^{0}}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{L^{\infty}(\omega_{0})}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{L^{2}(\omega_{0})}^{-1}$
is finite whenever $\widetilde{\mathds{U}}^{0}\subset L^{\infty}(\Omega)$ is
finite dimensional, so this is not a very strong assumption. The main take-
home for this theorem is that the computation of $a_{\operatorname{U}}$ and
$a_{\operatorname{E}}$ is typically very simple and clear given a particular
choice of
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
and $\operatorname{E}$. We also briefly remark that the Lipschitz constants in
this lemma do not need to be valid globally, only on the sequence
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}$. The same result holds
under a local-Lipschitz assumption, for example on the ball of radius
$\sup_{k\in\mathds{N}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
which is finite whenever $p\geq 0$.

##  6 L1 penalised reconstruction

The canonical example for FISTA is the LASSO problem with a quadratic data
fidelity and L1 regularisation. In this section we develop the necessary
analytical tools for the variant with general smooth fidelity term which will
be used for numerical results in Section 7. We consider three forms which will
be referred to as the continuous, countable, and discrete problem depending on
whether the space $\mathds{U}$ is $\mathcal{M}([0,1]^{d})$,
$\ell^{1}(\mathds{R})$, or $\mathds{R}^{M}$ respectively. We choose
$\mathds{H}$ to be $L^{2}([0,1]^{d}),\ \ell^{2}(\mathds{R}),$ or
$\mathds{R}^{M}$ correspondingly. Let
$\mathsf{A}\colon\mathds{U}\cap\mathds{H}\to\mathds{R}^{m}$ be a linear
operator represented by the kernels $\psi_{j}\in\mathds{H}$ such that

$\forall\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}\cap\mathds{H},\
j=1,\ldots,m,\qquad(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])_{j}=\left\langle\psi_{j},\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rangle.$
(17)

In the continuous case we will assume the additional smoothness $\psi_{j}\in
C^{1}([0,1]^{d})$. In Section 6.5 we will formally define and estimate several
operator semi-norms for $\mathsf{A}$ of this form, for example Lemma 8
confirms that $\mathsf{A}$ is continuous on $\mathds{H}$ (without loss of
generality ${\left\lVert\mathsf{A}\right\rVert}\leq 1$). In each case, the
energy we consider is written as

$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\eta)+\mu{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
(18)

for some $\mu>0$ where
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}={\left\lVert\cdot\right\rVert}_{1}$.
We assume $\operatorname{f}\in C^{1}(\mathds{R}^{m})$ is convex, bounded from
below, and $\nabla\operatorname{f}$ is 1-Lipschitz. Let
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$,
which is non-empty so long as $\psi_{j}\in C([0,1]^{d})$, see the proof of
(Bredies2013, , Prop. 3.1) when $\operatorname{f}$ is quadratic.

The aim of this section is to develop all of the necessary tools for
implementing Algorithm 1 on the energy (18) using the convergence guarantees
of either Theorem 4.2 or Theorem 4.3. This includes computing the rates
$a_{\operatorname{U}}$ and $a_{\operatorname{E}}$, estimating the continuous
gap $\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$, and
developing an efficient refinement choice for $\mathds{U}^{n}$. Below we will
just describe the form of $\widetilde{\mathds{U}}^{k}$ under the assumption
that $\mathds{U}^{n}\subset\widetilde{\mathds{U}}^{k}$ is chosen adaptively
for $n=n_{k-1}+1,\ldots,n_{k}$. The index $k$ refers to the scale or
resolution and $n$ refers to the iteration number of the reconstruction
algorithm.

###  6.1 Continuous case

We start by estimating rates in the case $\mathds{U}=\mathcal{M}(\Omega)$
where $\Omega=[0,1]^{d}$. In this case we choose $\widetilde{\mathds{U}}^{k}$
to be the span of all piecewise constant functions on a mesh of squares with
maximum side length $2^{-k}$ (i.e. $h=\tfrac{1}{2}$) and

$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\coloneqq\sum_{\omega\in\mathds{M}^{k}}\frac{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\omega)}{|\omega|}\mathds{1}_{\omega}\quad\text{where}\quad\mathds{1}_{\omega}(\vec{x})=\begin{cases}1&\vec{x}\in\omega\\\
0&\text{else}\end{cases}.$

By construction
$\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\in\widetilde{\mathds{U}}^{k}$,
however note that for any $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in
L^{1}(\Omega)$ and Dirac mass $\delta$ supported in $(0,1)^{d}$,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\delta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\sup_{\varphi\in
C(\Omega),{\left\lVert\varphi\right\rVert}_{L^{\infty}}\leq
1}\left\langle\varphi,\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\delta\right\rangle={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\delta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\geq
1=h^{0}.$ (19)

Because of this, application of Theorem 5.1 with $p=0$ gives
$a_{\operatorname{U}}=2^{\frac{d}{2}}$ but only $a_{\operatorname{E}}\geq 1$.
To improve our estimate of $a_{\operatorname{E}}$ requires additional
assumptions on $\mathsf{A}$. Note that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\sum_{\omega\in\mathds{M}^{k}}|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\omega)|\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$,
therefore we have

$\displaystyle\operatorname{E}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$
$\displaystyle=\operatorname{f}(\mathsf{A}\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\eta)-\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}-\eta)+\mu\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}-{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\right)$
(20)
$\displaystyle\leq\nabla\operatorname{f}(\mathsf{A}\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\eta)\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\mathsf{A}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$
(21)
$\displaystyle\leq\left[{\left\lVert\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}-\eta)\right\rVert}_{\ell^{2}}+{\left\lVert\mathsf{A}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\right\rVert}_{\ell^{2}}\right]{\left\lVert\mathsf{A}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\right\rVert}_{\ell^{2}}.$
(22)

as $\operatorname{f}$ is convex with 1-Lipschitz gradient. Clearly
${\left\lVert\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}-\eta)\right\rVert}_{\ell^{2}}$
is a constant. For the other term, for all $\vec{r}\in\mathds{R}^{m}$ denote
$\varphi\coloneqq\mathsf{A}^{*}\vec{r}$, then note that

$\vec{r}\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\mathsf{A}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})=\left\langle\varphi,\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rangle=\sum_{\omega\in\mathds{M}^{k}}\int_{\omega}\varphi(\vec{x})\mathop{}\\!\mathrm{d}[\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}]=\sum_{\omega\in\mathds{M}^{k}}|\omega|^{-1}\iint_{\omega^{2}}[\varphi(\vec{x})-\varphi(\vec{y})]\mathop{}\\!\mathrm{d}\vec{x}\mathop{}\\!\mathrm{d}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\vec{y}).$
(23)

With the pointwise bound
$|\varphi(\vec{x})-\varphi(\vec{y})|\leq\operatorname{diam}(\omega){\left\lVert\nabla\varphi\right\rVert}_{L^{\infty}}=\sqrt{d}2^{-k}{\left\lVert\nabla[\mathsf{A}^{*}\vec{r}]\right\rVert}_{L^{\infty}}$,
we deduce the estimate

${\left\lVert\mathsf{A}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\right\rVert}_{\ell^{2}}=\sup_{\vec{r}\in\mathds{R}^{m}}{\left\lVert\vec{r}\right\rVert}_{\ell^{2}}^{-1}\left\langle\mathsf{A}^{*}\vec{r},\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k}-\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rangle\leq\sqrt{d}2^{-k}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\sup_{\vec{r}\in\mathds{R}^{m}}{\left\lVert\vec{r}\right\rVert}_{\ell^{2}}^{-1}{\left\lVert\nabla[\mathsf{A}^{*}\vec{r}]\right\rVert}_{L^{\infty}}.$
(24)

In Lemma 9 we will show that this last term, which we denote the semi-norm
$|\mathsf{A}^{*}|_{\ell^{2}\to C^{1}}$, is bounded by
$\sqrt{m}\max_{j\in[m]}{\left\lVert\nabla\Psi_{j}\right\rVert}_{\infty}$. We
conclude that
$\operatorname{E}(\widetilde{\IfEqCase{2}{{0}{u}{1}{v}{2}{w}}[u2]}_{k})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\lesssim
2^{-k}$. In particular, this computation confirms two things. Firstly that the
scaling constant is $a_{\operatorname{E}}=2$, and secondly that the required
smoothness to achieve a good rate with Algorithm 1 is that
$\mathsf{A}^{*}\colon\mathds{R}^{m}\to C^{1}(\Omega)$ is a bounded operator.
This accounts for using the weaker topology of $\mathcal{M}(\Omega)$ rather
than $L^{1}(\Omega)$.

Inserting the computed rates into Theorem 4.2 or Theorem 4.3 gives the
guaranteed convergence rate

$\kappa=\frac{\log a_{\operatorname{U}}^{2}}{\log a_{\operatorname{E}}+\log
a_{\operatorname{U}}^{2}}=\frac{d}{1+d}\quad\implies\quad\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])\lesssim
n^{-2(1-\kappa)}=n^{-\frac{2}{1+d}}.$ (25)

This rate can be used to infer the required resolution at each iteration, in
particular on iteration $n$ with
$n^{2}\simeq(a_{\operatorname{E}}a_{\operatorname{U}}^{2})^{k}$ we expect the
resolution to be

$2^{-k}=\left(a_{\operatorname{E}}a_{\operatorname{U}}^{2}\right)^{\frac{k}{1+d}}\simeq
n^{-\frac{2}{1+d}}.$ (26)

###  6.2 Countable and discrete case

We now extend the rate computations to the case when
$\mathds{U}=\ell^{1}(\mathds{R})$, or a finite dimensional subspace. The key
fact here is that, even when $\mathds{U}$ is infinite dimensional, it is known
(e.g. (Unser2016, , Thm 6) and (Boyer2019, , Cor 3.8)) that there exists
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\operatorname*{argmin}_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
with at most $m$ non-zeros. If this is the case, then
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\in\ell^{2}(\mathds{R})$, indeed
${\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}_{\ell^{2}}\leq\sqrt{m}{\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}_{\ell^{1}}$.
This makes the estimates of $a_{\operatorname{E}}/a_{\operatorname{U}}$ much
simpler than in the continuous case as we can stay in the finite-dimensional
Hilbert-space setting.

For countable dimensions we consider discretisation subspaces of the form

$\widetilde{\mathds{U}}^{k}=\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\ell^{1}(\mathds{R})\operatorname{\;s.t.\;}i\notin
J_{k}\implies\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{i}=0\\}$

for some sets $J_{k}\subset\mathds{N}$, i.e. infinite vectors with finitely
many non-zeros. The key change in analysis from the continuous case is
${\left\lVert\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}\right\rVert}<\infty$, so
$a_{\operatorname{U}}=1$ and the expected rate of $n^{-2}$, independent of
$a_{\operatorname{E}}$ or any additional properties of $\mathsf{A}$. The
number of refinements will also be finite, therefore $n_{k}=\infty$ for some
$k$, the remaining conditions of Theorems 4.2 and 4.3 hold trivially.

###  6.3 Refinement metrics

Lemma 6 shows that adaptive refinement can be performed based on estimates of
the function gap or the subdifferential. In this subsection we provide
estimates for the forth case of Lemma 6 which can be easily computed. In this
case we consider
$\partial\operatorname{E}\colon\mathds{H}\rightrightarrows\mathds{H}$ so that
subdifferentials are well behaved, for example for explicit computation
assuming validity of the chain/sum rules for differentiation.

####  6.3.1 Bounds for discretised functionals

We start by computing estimates for discretised energies. This covers the
cases when either the continuous/countable energy is projected onto
$\mathds{U}^{n}$, or $\mathds{U}$ is finite dimensional. For notation we will
use the continuous case, to recover the other cases just replace continuous
indexing with discrete (i.e.
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0](\vec{x})\leadsto\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{i}$).

Let $\mathsf{\Pi}_{n}\colon\mathds{H}\to\mathds{U}^{n}$ denote the orthogonal
projection. We consider the discretised function
$\operatorname{E}|_{\mathds{U}^{n}}\colon\mathds{U}^{n}\to\mathds{R}$ and its
subdifferential
$\partial_{n}\operatorname{E}(\cdot)=\mathsf{\Pi}_{n}\partial\operatorname{E}(\cdot)$
on $\mathds{U}^{n}$. In our case, the behaviour of
$\operatorname{E}|_{\mathds{U}^{n}}$ is equivalent to replacing
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]$ with
$\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]$, and $\mathsf{A}^{*}$
with $\mathsf{\Pi}_{n}\mathsf{A}^{*}$.

##### Discrete gradient

We can use $\mathsf{\Pi}_{n}$ to compute the discrete subdifferential at
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}\in\mathds{U}^{n}$:

$\displaystyle\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})(\vec{x})$
$\displaystyle=[\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)](\vec{x})+\begin{cases}\\{+\mu\\}&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})>0\\\
[-\mu,\mu]&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})=0\\\
\\{-\mu\\}&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})<0\end{cases}$
(27)
$\displaystyle\eqqcolon[\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)](\vec{x})+\mu\mathsf{\Pi}_{n}\operatorname{sign}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x}))$
(28)

where we define $s+\mu[-1,1]=[s-\mu,s+\mu]$ for all $s\in\mathds{R}$, $\mu\geq
0$.

As
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}={\left\lVert\cdot\right\rVert}_{1}$,
the natural metric for $\partial_{n}\operatorname{E}$ is
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}={\left\lVert\cdot\right\rVert}_{\infty}$
which we can estimate

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max_{\vec{x}\in\Omega}\min_{\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]}\left\\{|\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]|\operatorname{\;s.t.\;}\IfEqCase{1}{{0}{u}{1}{v}{2}{w}}[u1]\in\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)(\vec{x})+\mu\mathsf{\Pi}_{n}\operatorname{sign}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x}))\right\\}$
(29)
$\displaystyle=\max_{\vec{x}\in\Omega}\begin{cases}|[\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)(\vec{x})+\mu|&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})>0\\\
|[\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)(\vec{x})-\mu|&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})<0\\\
\max\left(|\mathsf{\Pi}_{n}\mathsf{A}^{*}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)(\vec{x})|-\mu,0\right)&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})=0\end{cases}$
(30)

which can be used directly in Lemma 6.

##### Discrete gap

We now move on to the discrete gap,
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$.
This can be computed with a dual representation (e.g. Duval2017a ),

$\displaystyle\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\eta)+\mu{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle=\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\max_{\vec{\varphi}\in\mathds{R}^{m}}(\mathsf{A}\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\eta)\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}+\mu{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}-\operatorname{f}^{*}(\vec{\varphi})$
(31)
$\displaystyle=\max_{\vec{\varphi}\in\mathds{R}^{m}}\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}(\mathsf{A}\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\eta)\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}+\mu{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}-\operatorname{f}^{*}(\vec{\varphi})$
(32)
$\displaystyle=\max_{\vec{\varphi}\in\mathds{R}^{m}}\begin{cases}-\eta\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}-\operatorname{f}^{*}(\vec{\varphi})&\qquad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}\leq\mu\\\
-\infty&\qquad\text{else}\end{cases}$ (33)
$\displaystyle=-\min_{\vec{\varphi}\in\mathds{R}^{m}}\underbrace{\operatorname{f}^{*}(\vec{\varphi})+\eta\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}}_{\eqqcolon\operatorname{E}^{\dagger}(\vec{\varphi})}+\chi({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}\leq\mu).$
(34)

In particular,

$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])-\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}^{n}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])+\min_{\vec{\varphi}\in\mathds{R}^{m}\operatorname{\;s.t.\;}{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{*}\leq\mu}\operatorname{E}^{\dagger}(\vec{\varphi})\leq\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])+\operatorname{E}^{\dagger}(\vec{\varphi})$
(35)

for any feasible $\vec{\varphi}\in\mathds{R}^{m}$. We further derive the
criticality condition, if
$(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*},\vec{\varphi}^{*})$ is a saddle
point, then

$\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}-\eta\in\partial\operatorname{f}^{*}(\vec{\varphi}^{*}),\qquad\text{or
equivalently}\qquad\vec{\varphi}^{*}=\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}-\eta).$
(36)

We remark briefly that $\operatorname{E}^{\dagger}$ should be thought of as
the dual of $\operatorname{E}$ but without the constraint. We choose to omit
it here to highlight that it is only the constraint which changes between the
discrete and continuous cases; the value of $\operatorname{E}^{\dagger}$ will
remain the same.

Given $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}\in\mathds{U}^{n}$, the
optimality condition motivates a simple rule for choosing $\vec{\varphi}$:

$\vec{\varphi}_{n}\coloneqq\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta),\qquad\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])-\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{\prime}\in\mathds{U}^{n}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{\prime})\leq\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n})$
(37)

for some
$0\leq\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\leq\frac{\mu}{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}_{*}}$.
In the case
$\operatorname{f}(\cdot)=\frac{1}{2}{\left\lVert\cdot\right\rVert}_{\ell^{2}}^{2}$,
one can use the optimal choice

$\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]=\max\left(0,\min\left(\frac{-\eta\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}_{n}}{{\left\lVert\vec{\varphi}_{n}\right\rVert}_{\ell^{2}}^{2}},\frac{\mu}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}}\right)\right).$
(38)

To apply Algorithm 1, we are assuming that both
$\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)$
and
$\mathsf{\Pi}_{n}\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)$
are easily computable, therefore
$\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]$ and
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n})$
are also easy to compute.

####  6.3.2 Bounds for countable functionals

Extending the results of Section 6.3.1 to $\mathds{U}=\ell^{1}(\mathds{R})$ is
analytically very simple but computationally relies heavily on the specific
choice of $\mathsf{A}$. The computations of subdifferentials and gaps carry
straight over replacing $\mathsf{\Pi}_{n}$ with the identity and adding the
sets $J_{n}\subset\mathds{N}$ which define
$\mathds{U}^{n}=\\{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\ell^{1}\operatorname{\;s.t.\;}i\notin
J_{n}\implies\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{i}=0\\}$. Recall that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}\coloneqq\inf_{s\in\operatorname{sign}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}+\mu
s\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$ where the
$\operatorname{sign}$ function has the pointwise set-valued definition as
indicated in (27)-(28). Where $[u_{n}]_{i}=0$, the choice
$s_{i}=\min(1,\max(-1,-\mu^{-1}[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}))$
achieves the minimal value

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max_{i\in\mathds{N}}\begin{cases}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}+\mu|&[\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}]_{i}>0\\\
|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}-\mu|&[\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}]_{i}<0\\\
\max\left(|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|-\mu,0\right)&[\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}]_{i}=0\end{cases}$
(39)
$\displaystyle\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
$\displaystyle\leq\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}),\qquad\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\in\left[0,\frac{\mu}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}}\right]$
(40)

where
$\vec{\varphi}_{n}=\nabla\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}-\eta)\in\mathds{R}^{m}$
is always exactly computable.

In the countable case, the sets $J_{n}$ give a clear partition into
known/unknown values in these definitions. For $i\in J_{n}$ the computation is
the same as in Section 6.3.1, then for $i\notin J_{n}$ we know
$[\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}]_{i}=0$ which simplifies the
remaining computations. This leads to:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max\left(\max_{i\in
J_{n}}|[\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})]_{i}|,\
\sup_{i\notin
J_{n}}|[\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})]_{i}|\right)=\max\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*},\
\sup_{i\notin J_{n}}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|-\mu\right)$ (41)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max\left(\max_{i\in
J_{n}}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|,\ \sup_{i\notin
J_{n}}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|\right)\hskip
19.0pt=\max\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*},\
\sup_{i\notin J_{n}}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|\right).$ (42)

Both estimates only rely on an upper bound of $\max_{i\notin
J_{n}}|[\mathsf{A}^{*}\vec{\varphi}_{n}]_{i}|$. One example computing this
value is seen in Section 7.2.

####  6.3.3 Bounds for continuous functionals

Finally we extend the results of Section 6.3.1 to continuous problems. Similar
to the countable case (39)-(40), the exact formulae can be written down
immediately:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max_{\vec{x}\in\Omega}\begin{cases}|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})+\mu|&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})>0\\\
|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})-\mu|&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})<0\\\
\max\left(|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})|-\mu,0\right)&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}(\vec{x})=0\end{cases}$
(43)
$\displaystyle\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\inf_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{H}}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])$
$\displaystyle\leq\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}),\qquad\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\in\left[0,\frac{\mu}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}}\right]$
(44)

with $\operatorname{E}^{\dagger}$ as defined in (34). Recall that there is a
mesh $\mathds{M}^{n}$ corresponding to $\mathds{U}^{n}$ such that
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ is constant on each
$\omega\in\mathds{M}^{n}$, so we can rewrite these bounds:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max_{\omega\in\mathds{M}^{n}}\begin{cases}{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}+\mu\right\rVert}_{L^{\infty}(\omega)}&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}|_{\omega}>0\\\
{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}-\mu\right\rVert}_{L^{\infty}(\omega)}&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}|_{\omega}<0\\\
\max(0,{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{\infty}(\omega)}-\mu)&\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}|_{\omega}=0\end{cases}$
(45)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max_{\omega\in\mathds{M}^{n}}{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{\infty}(\omega)}.$
(46)

Now, both values can be estimated relying on pixel-wise supremum norms of
$\mathsf{A}^{*}\vec{\varphi}_{n}$ which we have assumed is sufficiently
smooth. We will therefore use a pixel-wise Taylor expansion to provide a
simple and accurate estimate. For instance, let $\vec{x}_{i}$ be the midpoint
of the pixel $\omega$, then

${\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{\infty}(\omega)}\leq|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})|+\frac{\operatorname{diam}(\omega)}{2}|[\nabla\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})|+\frac{\operatorname{diam}(\omega)^{2}}{8}|\mathsf{A}^{*}\vec{\varphi}_{n}|_{C^{2}}.$
(47)

In this work we chose a first order expansion because we are looking for
extrema of $\mathsf{A}^{*}\vec{\varphi}_{n}$, i.e. we are most interested in
the squares $\omega$ such that

$|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})|\approx\mu,\qquad|[\nabla\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})|\approx
0,\qquad[\nabla^{2}\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})\preceq 0.$
(48)

A zeroth order expansion would be optimally inefficient (approximating
$|[\nabla\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x}_{i})|$ with
$|\mathsf{A}^{*}\vec{\varphi}_{n}|_{C^{1}}$) and a second order expansion
would possibly be more elegant but harder to implement. We found that a first
order expansion was simple and efficient.

The bounds presented here for continuous problems emphasise the twinned
properties required for adaptive mesh optimisation. The mesh should be refined
greedily to the structures of $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$, but
also must be sufficiently uniform to provide a good estimate for
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$. This is a
classical exploitation/exploration trade-off; exploiting visible structure
whilst searching for other structures which are not yet visible.

###  6.4 Support detection

The main motivation for using L1 penalties in applications is because it
recovers sparse signals, in the case of compressed sensing the support of
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$ is also provably close to the
“true” support Duval2017a ; Poon2018 . If
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}\approx\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$
in the appropriate sense, then we should also be able to quantify the
statement
$\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\approx\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$.
Such methods are referred to as _safe screening_ rules ElGhaoui2010  which
gradually identify the support and allow the optimisation algorithm to
constrain parts of the reconstruction to 0. In this subsection we propose a
new simple screening rule which is capable of generalising to our continuous
subspace approximation setting. It is likely that more advanced methods
Bonnefoy2015 ; Ndiaye2017  can also be adapted, although that is beyond the
scope of this work. The key difference is the allowance of inexact
computations resulting from estimates such as (47).

The support of $\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$ has already been
characterised very precisely Duval2017a ; Poon2018 . In particular, the
support is at most $m$ distinct points and are a subset of
$\\{\vec{x}\in\Omega\operatorname{\;s.t.\;}|\mathsf{A}^{*}\vec{\varphi}^{*}|(\vec{x})=\mu\\}$
(an equivalent statement holds for the countable case). Less formally, this
can also be seen from the the subdifferential computations in Section 6.3, for
all $\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$
we have

$0\in\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})(\vec{x})=[\mathsf{A}^{*}\vec{\varphi}^{*}](\vec{x})+\mu\operatorname{sign}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\vec{x})).$
(49)

Heuristically, we will use strong convexity of $\operatorname{E}^{\dagger}$
from (34) and smoothness of $\mathsf{A}^{*}$ to quantify the statement:

$\text{if}\quad\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})\approx
0\quad\text{then}\quad\left\\{\vec{x}\operatorname{\;s.t.\;}|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})|\ll\mu\right\\}\subset\\{\vec{x}\operatorname{\;s.t.\;}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\vec{x})=0\\}.$

Recall that $\nabla\operatorname{f}$ is 1-Lipschitz if and only if
$\operatorname{f}^{*}$ is 1-strongly convex (Hiriart2013, , Chapter 10, Thm.
4.2.2). Therefore, if
$\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}$
and $\vec{\varphi}^{*}$ are both dual-feasible, then

$\tfrac{1}{2}{\left\lVert\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}-\vec{\varphi}^{*}\right\rVert}_{\ell^{2}}^{2}\leq\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})-\operatorname{E}^{\dagger}(\vec{\varphi}^{*})=\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})+\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\leq\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})+\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}),$
(50)

which gives an easily computable bound on
${\left\lVert\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}-\vec{\varphi}^{*}\right\rVert}_{\ell^{2}}$.
Now we estimate $\mathsf{A}^{*}\vec{\varphi}_{n}$ on the support of
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$:

$\displaystyle\min_{\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})}|[\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})|$
$\displaystyle\geq\min_{\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})}|[\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})|$
(51)
$\displaystyle=\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}^{-1}\min_{\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})}|[\mathsf{A}^{*}\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}](\vec{x})|$
(52)
$\displaystyle\geq\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}^{-1}\min_{\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})}|[\mathsf{A}^{*}\vec{\varphi}^{*}](\vec{x})|-|[\mathsf{A}^{*}\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}-\mathsf{A}^{*}\vec{\varphi}^{*}](\vec{x})|$
(53)
$\displaystyle=\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}^{-1}\min_{\vec{x}\in\operatorname{supp}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})}\mu-|[\mathsf{A}^{*}\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}-\mathsf{A}^{*}\vec{\varphi}^{*}](\vec{x})|$
(54)
$\displaystyle\geq\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}^{-1}\left(\mu-|\mathsf{A}^{*}|_{\ell^{2}\to
L^{\infty}}{\left\lVert\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}-\vec{\varphi}^{*}\right\rVert}_{\ell^{2}}\right).$
(55)

Therefore,

$|[\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}](\vec{x})|<\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}^{-1}\left(\mu-\sqrt{2(\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}))}|\mathsf{A}^{*}|_{\ell^{2}\to
L^{\infty}}\right)\qquad\implies\qquad\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}(\vec{x})=0.$
(56)

This equation is valid when $\vec{x}$ is either a continuous or countable
index, the only distinction is to switch to $\ell^{\infty}$ in the norm of
$\mathsf{A}^{*}$. To make the equivalent statement on the discretised problem,
simply replace
$\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}$
with $\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]$
and $\mathsf{A}^{*}$ with $\mathsf{\Pi}_{n}\mathsf{A}^{*}$. There are two
short observations on this formula:

  * •

The convergence guarantee from Theorem 4.2 is for the primal gap
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$,
rather than the primal-dual gap
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})$
used here. Although there is no guaranteed rate for the primal-dual gap, it is
much more easily computable than the primal gap.

  * •

In Section 6.1, $|\mathsf{A}^{*}|_{\ell^{2}\to C^{1}}<\infty$ was required to
compute a rate of convergence for
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$, but only
$|\mathsf{A}^{*}|_{\ell^{2}\to L^{\infty}}<\infty$ is needed to estimate the
support.

###  6.5 Operator norms

For numerical implementation of (18), we are required to accurately estimate
several operator norms of $\mathsf{A}$ of the form in (17). In particular,
there are kernels $\psi_{j}\in\mathds{H}$ such that
$(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])_{j}=\left\langle\psi_{j},\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rangle$
for each $j\in[m]$. Verifying that ${\left\lVert\mathsf{A}\right\rVert}\leq 1$
can be performed by computing
$|\mathsf{A}\mathsf{A}^{*}|_{\ell^{2}\to\ell^{2}}$, and the adaptivity
described in Sections 6.1, 6.3.3, and 6.4 requires the values of
$|\mathsf{A}^{*}|_{\ell^{2}\to L^{\infty}}$, $|\mathsf{A}^{*}|_{\ell^{2}\to
C^{1}}$, and $|\mathsf{A}^{*}|_{\ell^{2}\to C^{2}}$. The aim for this section
is to provide estimates of these norms and seminorms for the numerical
examples presented in Section 7.

The following lemma allows for exact computation of the operator norm of
$\mathsf{A}$.

###### Lemma 8

If $\mathsf{A}\colon\mathds{H}\to\mathds{R}^{m}$ has kernels
$\psi_{j}\in\mathds{H}$ for $j\in[m]$, then
$\mathsf{A}\mathsf{A}^{*}\in\mathds{R}^{m\times m}$ has entries
$(\mathsf{A}\mathsf{A}^{*})_{i,j}=\left\langle\psi_{i},\psi_{j}\right\rangle$,
so the spectral norm
${\left\lVert\mathsf{A}^{*}\mathsf{A}\right\rVert}={\left\lVert\mathsf{A}\mathsf{A}^{*}\right\rVert}$
can be computed efficiently.

###### Proof

To compute the entries of
$\mathsf{A}\mathsf{A}^{*}\colon\mathds{R}^{m}\to\mathds{R}^{m}$, observe that
for any $\vec{r}\in\mathds{R}^{m}$

$(\mathsf{A}\mathsf{A}^{*}\vec{r})_{i}=\left\langle\psi_{i},\mathsf{A}^{*}\vec{r}\right\rangle=\left\langle\psi_{i},\sum_{j=1}^{m}r_{j}\psi_{j}\right\rangle=\sum_{j=1}^{m}\left\langle\psi_{i},\psi_{j}\right\rangle
r_{j}$ (57)

as required. ∎

If ${\left\lVert\mathsf{A}^{*}\mathsf{A}\right\rVert}$ is not analytically
tractable, then Lemma 8 enables it to be computed using standard finite
dimensional methods. The operator $\mathsf{A}\mathsf{A}^{*}$ is always finite
dimensional, and can be computed without discretisation error.

In the continuous case, when $\mathds{H}=L^{2}(\Omega)$ we also need to
estimate the smoothness properties of $\mathsf{A}^{*}$. A generic result for
this is given in the following lemma.

###### Lemma 9

If $\mathsf{A}\colon L^{2}([0,1]^{d})\to\mathds{R}^{m}$ has kernels
$\psi_{j}\in L^{2}(\Omega)\cap C^{k}(\Omega)$ for $j\in[m]$, then for all
$\frac{1}{q}+\frac{1}{q^{*}}=1$, $q\in[1,\infty]$, we have

$\displaystyle|\mathsf{A}^{*}\vec{r}|_{C^{k}}$
$\displaystyle\coloneqq\sup_{\vec{x}\in\Omega}|\nabla^{k}[\mathsf{A}^{*}\vec{r}]|(\vec{x})\leq\sup_{\vec{x}\in\Omega}{\left\lVert(\nabla^{k}\psi_{j}(\vec{x}))_{j=1}^{m}\right\rVert}_{\ell^{q^{*}}}{\left\lVert\vec{r}\right\rVert}_{\ell^{q}},$
(58) $\displaystyle|\mathsf{A}^{*}|_{\ell^{2}\to C^{k}}$
$\displaystyle\coloneqq\sup_{{\left\lVert\vec{r}\right\rVert}_{\ell^{2}}\leq
1}|\mathsf{A}^{*}\vec{r}|_{C^{k}}\leq\sup_{\vec{x}\in\Omega}{\left\lVert(\nabla^{k}\psi_{j}(\vec{x}))_{j=1}^{m}\right\rVert}_{\ell^{q^{*}}}\times\begin{cases}1&q\geq
2\\\ \sqrt{m^{2-q}}&q<2\end{cases}.$ (59)

###### Proof

For the first inequality, we apply the Hölder inequality on $\mathds{R}^{m}$:

$|\nabla^{k}[\mathsf{A}^{*}\vec{r}]|(\vec{x})=\left|\sum_{j=1}^{m}\nabla^{k}\psi_{j}(\vec{x})r_{j}\right|\leq\left(\sum_{j=1}^{m}|\nabla^{k}\psi_{j}(\vec{x})|^{q^{*}}\right)^{\frac{1}{q^{*}}}{\left\lVert\vec{r}\right\rVert}_{\ell^{q}}={\left\lVert(\nabla^{k}\psi_{j}(\vec{x}))_{j}\right\rVert}_{\ell^{q^{*}}}{\left\lVert\vec{r}\right\rVert}_{\ell^{q}}\;.$

For the second inequality, if $q\geq 2$ and $\sum_{j=1}^{m}r_{j}^{2}\leq 1$,
then $|r_{j}|\leq 1$ for all $j$ and
${\left\lVert\vec{r}\right\rVert}_{\ell^{q}}^{q}\leq{\left\lVert\vec{r}\right\rVert}_{\ell^{2}}^{2}\leq
1$. If $q<2$ and ${\left\lVert\vec{r}\right\rVert}_{\ell^{2}}\leq 1$, then we
again use Hölder’s inequality:

$\sum_{j=1}^{m}r_{j}^{q}\leq\Big{(}\sum_{j=1}^{m}1^{Q^{*}}\Big{)}^{\frac{1}{Q^{*}}}\Big{(}\sum_{j=1}^{m}r_{j}^{qQ}\Big{)}^{\frac{1}{Q}}\leq
m^{\frac{2-q}{2}}$

for $Q=\frac{2}{q}$. ∎

The examples in Section 7 require explicit computations of the expressions in
Lemmas 8 and 9. These computations are provided in the appendix, Theorem C.1.

##  7 Numerical examples

We present four numerical examples. The first two are in 1D to demonstrate the
performance of different variants of Algorithm 1, both with and without
adaptivity. In particular, we explore sparse Gaussian deconvolution and sparse
signal recovery from Fourier data. We compare with the _continuous basis
pursuit_ (CBP) discretisation Ekanadham2011 ; Duval2017b  which is also
designed to achieve super-resolution accuracy within a convex framework. More
details of this method will be provided in Section 7.1.

The next example is 2D reconstruction from Radon or X-ray data with wavelet-
sparsity and a robust data fidelity. As the forward operator is not
sufficiently smooth, we must optimise in $\ell^{1}(\mathds{R})$, which
naturally leads to the choice of a wavelet basis.

Finally, we process a dataset which represents a realistic application in
biological microscopy, referred to as STORM microscopy. In essence, the task
is to perform 2D Gaussian de-blurring/super-resolution and denoising to find
the location of sparse spikes of signal.

In this section, the main aim is to minimise
$\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$,
and so this will be our main metric for the success of an algorithm, referred
to as the “continuous gap”. Lemma 6 only provides guarantees on the values of
$\min_{n\leq N}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})$
so it is this monotone estimate which is plotted. As
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})$ is not known
exactly, we always use the estimate $\min_{n\leq
N}\operatorname{E}_{0}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\approx\min_{n\leq
N}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\min_{n^{\prime}\leq
n}\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n^{\prime}})$.
Another quantity of interest is minimisation of the discrete energy
$\min_{n\leq
N}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\min_{n^{\prime}\leq
n}\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n^{\prime}})$
which will be referred to as the “discrete gap”. Note that for the adaptive
schemes the discrete gap may not be monotonic as the discrete dual problem
changes with $N$.

The code to reproduce these examples can be found
online111https://github.com/robtovey/2020SpatiallyAdaptiveFISTA.

###  7.1 1D continuous LASSO

In this example we choose $\mathds{U}=\mathcal{M}([0,1])$,
$\mathds{H}=L^{2}([0,1])$,
$\operatorname{f}(\cdot)=\frac{1}{2}{\left\lVert\cdot\right\rVert}_{\ell^{2}}^{2}$
and $\mathsf{A}\colon\mathds{U}\to\mathds{R}^{30}$ with either random Fourier
kernels:

$(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])_{j}=\int_{0}^{1}\cos(\IfEqCase{3}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k3]_{j}x)\mathop{}\\!\mathrm{d}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0](x),\qquad\IfEqCase{3}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k3]_{j}\sim\operatorname{Uniform}[-100,100],\
j=1,2,\ldots,30,\ \mu=0.02,$ (60)

or Gaussian kernels on a regular grid:

$(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0])_{j}=(2\pi\sigma^{2})^{-\frac{1}{2}}\int_{0}^{1}\exp\left(-\frac{(x-(j-1)\Delta)^{2}}{2\sigma^{2}}\right)\mathop{}\\!\mathrm{d}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0](x),\quad\sigma=0.12,\
\Delta=\tfrac{1}{29},\ j=1,2,\ldots,30,\ \mu=0.06.$ (61)

Several variants of FISTA are compared for these examples but the key
alternative shown here is the CBP discretisation. For this choice of
$\operatorname{f}$, we call (18) the continuous LASSO problem, for which there
are many numerical methods (c.f. Bredies2013 ; Castro2016 ; Boyd2017 ;
Catala2019 ) however, most require the solution of a non-convex problem. We
have focused on CBP because it approximates
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$ through a convex discrete
optimisation problem which is asymptotically exact in the limit $h\to 0$. It
can also be optimised with FISTA which allows for direct comparison with the
uniform and adaptive mesh approaches. The idea is that for a fixed mesh, the
kernels of $\mathsf{A}$ are expanded to first order on each pixel and a
particular first order basis is also chosen Ekanadham2011 ; Duval2017b . If
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$ has only one Dirac spike in each
pixel, then the zeroth order information should correspond to the mass of the
spike, and additional first order information should determine the location.

As shown in Section 6, in 1D we have
$a_{\operatorname{U}}=a_{\operatorname{E}}=2$. The estimates given in (25) and
(26) in dimension $d=1$ predict that the adaptive energy will decay at a rate
of
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\lesssim\frac{1}{n}$
so long as the pixel size also decreases at a rate of $h\sim\frac{1}{n}$. To
achieve these rates, we implement a refinement criterion from Lemma 6 with
guarantee of
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n_{k}-1})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\lesssim
2^{-k}$ using the estimates made in Section 6.3. We choose subspaces
$\mathds{U}^{n}$ to approximately enforce

$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})\leq
2(\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})+\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n})),$
(62)

i.e. the continuous gap is bounded by twice the discrete gap. In particular,
note that for
$\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\approx\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]$,

$\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n})=\tfrac{1}{2}{\left\lVert\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\vec{\varphi}_{n}\right\rVert}^{2}+\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}\eta\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}_{n}=\frac{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}}{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]}\left(\frac{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}}{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]}\tfrac{1}{2}{\left\lVert\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n}\right\rVert}^{2}+\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\eta\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\vec{\varphi}_{n}\right)\approx\frac{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}}{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]}\operatorname{E}^{\dagger}(\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]\vec{\varphi}_{n}).$
(63)

Converting this into a spatial refinement criteria, recall

$\frac{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]_{0}}{\IfEqCase{2}{{0}{\alpha}{1}{\beta}{2}{\gamma}{3}{a}{4}{b}{5}{c}}[k2]}\approx\frac{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}}=\frac{\max_{\omega\in\mathds{M}^{n}}{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{\infty}(\omega)}}{\max_{\omega\in\mathds{M}^{n}}|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}(\omega)|}\approx\max_{\omega\in\mathds{M}^{n}}\frac{{\left\lVert\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{\infty}(\omega)}}{|\mathsf{\Pi}_{n}\mathsf{A}^{*}\vec{\varphi}_{n}(\omega)|}$
(64)

is the maximum ratio of second vs. zeroth order Taylor approximations of
$\mathsf{A}^{*}\vec{\varphi}_{n}$ on pixel $\omega$. This was found to be an
efficient method of selecting pixels for refinement using quantities which had
already been computed. Note briefly that this greedy strategy directly targets
uncertainty, refinements also happen outside of the support of $u_{n}$ to
guarantee that this is representative of
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$. Such refinement is necessary to
avoid discrete minimisers of $\operatorname{E}$ which are not global
minimisers.

Figure 1: Rates of continuous/discrete gap convergence for different LASSO
algorithms with 128, 256, or 512 pixels. The “adaptive” method uses the
proposed algorithm. Both “fixed” and “CBP” use standard FISTA with a uniform
discretisation.

Figure 2: Convergence plots for solving 1D problems with different algorithms.
“Adaptive” methods use Algorithm 1 with fewer than 1024 pixels and the
remaining methods use a uniform discretisation of 1024 pixels.

Figure 3: Example reconstruction from the algorithms considered in Fig. 3.
Pixel boundaries are indicated on the $x$-axis and the filtering method of
Section 6.4 allows us to exclude the red shaded regions from
$\operatorname{supp}(u^{*})$. Values on the $y$-axis are normalised to units
of mass, i.e. a Dirac mass would have height 1.

##### Comparison of discretisation methods

In Fig. 3 we compare the three core approaches: fixed uniform discretisation,
adaptive discretisation, and CBP. In particular, we wish to observe their
convergence properties as the number of pixels is allowed to grow. In each
case we use a FISTA stepsize of $t_{n}=\frac{n+19}{20}$. The adaptive
discretisation is started with one pixel and limited to 128, 256, or 512
pixels while the fixed and CBP discretisations have uniform discretisations
with the maximum number of pixels. The main observations are:

  * •

The adaptive scheme is much more efficient, in both examples the adaptive
scheme with 128 pixels is at least as good as both fixed discretisations with
512 pixels. In fact, only a maximum of 214 pixels were needed by the adaptive
method in either example.

  * •

With Fourier kernels the uniform piecewise constant discretisation is more
efficient than CBP but in the Gaussian case this is reversed. This suggests
that the performance of CBP depends on the smoothness of $\mathsf{A}$.

  * •

The discrete gaps for non-adaptive optimisation behave as is common for FISTA,
initial convergence is polynomial until a locally linear regime activates
Tao2016 . CBP is always slower to converge than the piecewise constant
discretisation.

  * •

The adaptive refinement criterion succeeds in keeping the continuous/discrete
gaps close for all $n$, i.e. (62).

It is not completely fair to judge CBP with the continuous gap because,
although it generates a continuous representation, this continuous
representation is not necessarily consistent with the discrete gap being
optimised, unlike when discretised with finite element methods. On the other
hand, this is still the intended interpretation of the algorithm and we have
no more appropriate metric for success in this case.

##### Comparison of FISTA variants

Fig. 3 compares many methods with either fixed or adaptive discretisations.
Each adaptive scheme is allowed up to 1024 pixels and each uniform
discretisation uses exactly 1024. An example of each reconstruction method is
shown in Fig. 3. The adaptive method better identifies the support of $u^{*}$
and clearly localises pixels on that support. The reconstruction with uniform
grid fails to provably identify the support of
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}$, despite having found a
qualitatively accurate discrete minimiser. The “Greedy FISTA” implementation
was proposed by in Liang2018  and we include the adaptive variant despite a
lack of convergence proof. The remaining FISTA algorithms use a FISTA time
step of $t_{n}=\frac{n+a-1}{a}$ for the given value of $a$, as proposed in
Chambolle2015 . In this example CBP used the greedy FISTA implementation which
gave faster observed convergence. Fig. 3 compares the discrete gaps because it
is the accurate metric for fixed discretisations, and for the adaptive
discretisation it should also be an accurate predictor of the continuous gap.
The main observations are:

  * •

Each algorithm displays very similar convergence properties. The main
difference is that the reconstructions with fixed discretisations accelerate
after $10^{4}$-$10^{5}$ iterations.

  * •

During the initial “slow” phase, adaptive and fixed discretisations appear to
achieve very similar (discrete) convergence rates. The coarse-to-fine
adaptivity is not slower than fixed discretisations in this regime.

  * •

Lemma 6 accurately predicts the $\frac{1}{n}$ rate of the adaptive methods,
mirrored in the fixed discretisations. This suggests that high-resolution
discretisations are also initially limited by this $\frac{1}{n}$ rate before
entering the asymptotic regime, consistent with (14).

  * •

The fastest FISTA stepsize choice is consistently the greedy variant, although
$a=20$ is very comparable.

  * •

While each adaptive algorithm is allowed to use up to 1024 pixels, in Fig. 3
the most used was 235.

##### Comparison of fixed and adaptive discretisation

Motivated by the findings in Fig. 3, we now look more closely at the
performance of the $a=20$ and the greedy FISTA schemes. We have convergence
results for the former, but the latter typically performs the best for non-
adaptive optimisation and is never worse than $a=20$ in the adaptive setting.
The question is whether it is faster/more efficient to use the proposed
adaptive scheme, or to use a classical scheme at sufficiently high uniform
resolution. The fixed discretisations use 1024 pixels (i.e. constant pixel
size of $2^{-10}$ in Fig. 6) and the adaptive discretisation starts with two
pixels with an upper limit of 1024. As expected, the fixed discretisation
starts with a smaller continuous gap before plateauing to a sub-optimal gap
around
$\operatorname{E}=\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})+0.1$.

Fig. 6 shows convergence of pixel size and continuous gap with respect to
number of iterations. Fig. 6 shows the more practical attributes of continuous
gap and number of pixels against execution time. We see that the adaptive
discretisation is consistently capable of computing lower energies with fewer
pixels and in less time than the uniform discretisation. The convergence
behaviour is very consistent with respect to number of iterations.

Suppose that the numerical aim is to find a function
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n}$ with
$\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})-\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*})\leq
0.1$, all methods would converge after $O(10^{3})$ iterations, demonstrating
some equivalence between the two FISTA algorithms. For $n\in[10^{3},10^{4}]$,
in both problems, the adaptive schemes coincide with the fixed schemes in both
energy and minimum pixel size. On the other hand, we also see that the
adaptive scheme achieves this energy in almost an order of magnitude less time
and fewer pixels.

Figure 4: Continuous convergence of adaptive (coarse-to-fine pixel size)
compared with uniform discretisation (constant pixel size) with respect to
number of iterations.

Figure 5: Continuous convergence of adaptive compared with uniform
discretisation with respect to wall-clock time and total number of pixels
(memory requirement).

$\begin{aligned} J_{n}&=\\{(0,0),(0,1),(0,2),(1,2),(1,1)\\}\\\
\operatorname{leaf}(J_{n})&=\\{(0,2),(1,2),(1,1)\\}\\\
\mathds{M}^{n}&=\left\\{[0,\tfrac{1}{4}),[\tfrac{1}{4},\tfrac{1}{2}),[\tfrac{1}{2},1)\right\\}\end{aligned}$
$(0,0)$ $[0,1]$ $(0,1)$ $[0,\frac{1}{2}]$ $(0,2)$ $[0,\frac{1}{4}]$ $(1,2)$
$[\frac{1}{4},\frac{1}{2}]$ $(1,1)$ $[\frac{1}{2},1]$

Figure 6: Example tree representation of 1D wavelets. Left: nodes, leaves, and
mesh of discretisation. Right: arrangement into a tree with index $(j,k)$ and
corresponding support of wavelet $w_{j,k}$ underneath.

###  7.2 2D robust sparse wavelet reconstruction

In this example we consider $\mathsf{A}$ to be a 2D Radon transform. In
particular, the rows of $\mathsf{A}$ correspond to integrals over the sets
$\mathds{X}^{I}_{i}$ where

$\mathds{X}_{i}^{I}=\left\\{\vec{x}\in[-\tfrac{1}{2},\tfrac{1}{2}]^{2}\operatorname{\;s.t.\;}\vec{x}\vbox{\hbox{\leavevmode\resizebox{6.0pt}{}{$\cdot$}}}\begin{pmatrix}\cos\theta_{I}\\\
\sin\theta_{I}\end{pmatrix}\in\left[-\tfrac{1}{2}+\tfrac{i-1}{100},-\tfrac{1}{2}+\tfrac{i}{100}\right)\right\\},\quad\theta_{I}=\frac{180^{\circ}}{51}I$
(65)

for $i=\in[100]$, $I\in[50]$. This is not exactly in the form analysed by
Theorem C.1, only the sets
$\\{\mathds{X}^{I}_{i}\operatorname{\;s.t.\;}i\in[100]\\}$ for each $I$ are
disjoint, therefore we apply Theorem C.1 block-wise to estimate

${\left\lVert\mathsf{A}\right\rVert}_{L^{2}\to\ell^{2}}\leq\sqrt{\sum_{I\in[50]}\max_{i\in[100]}|\mathds{X}^{I}_{i}|}=\sqrt{\sum_{I\in[50]}\max_{i\in[100]}\int_{\mathds{X}^{I}_{i}}1\mathop{}\\!\mathrm{d}\vec{x}}=\sqrt{\sum_{I\in[50]}\max_{i\in[100]}\
(\mathsf{A}\mathds{1})_{i,I}}\;.$ (66)

$\mathsf{A}$ is not smooth, therefore we can’t bound
$|\mathsf{A}^{*}|_{C^{k}}$ for $k>0$, and so we must look to minimise over
$\ell^{1}$ rather than $L^{1}$. The natural choice is to promote sparsity in a
wavelet basis which can be rearranged into the form of (18):

$\min_{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\in\mathds{U}}\operatorname{f}(\mathsf{A}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]-\eta)+\mu{\left\lVert\mathsf{W}^{-1}\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]\right\rVert}_{\ell^{1}}=\min_{\widehat{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}\in\ell^{1}(\mathds{R})}\operatorname{f}(\mathsf{A}\mathsf{W}\widehat{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}-\eta)+\mu{\left\lVert\widehat{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}\right\rVert}_{\ell^{1}}.$
(67)

The minimisers are related by
$\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]^{*}=\mathsf{W}\widehat{\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]}^{*}$
and, for wavelet bases, $\mathsf{W}$ is orthonormal so
${\left\lVert\mathsf{A}\mathsf{W}\right\rVert}_{\ell^{2}\to\ell^{2}}={\left\lVert\mathsf{A}\right\rVert}_{L^{2}\to\ell^{2}}$.
In this example we consider the smoothed robust fidelity Rosset2007

$\operatorname{f}(\vec{\varphi})=\sum_{i=1}^{m}\begin{cases}10^{-4}|\varphi_{i}|&|\varphi_{i}|\geq
10^{-4}\\\
\tfrac{1}{2}|\varphi_{i}|^{2}+\tfrac{1}{2}10^{-8}&\text{else}\end{cases}\approx
10^{-4}{\left\lVert\vec{\varphi}\right\rVert}_{\ell^{1}}.$ (68)

From Section 6.3 we know that to track convergence and perform adaptive
refinement, it is sufficient to accurately bound
$|[\mathsf{W}^{\top}\mathsf{A}^{*}\vec{\varphi}_{n}]_{j}|$ for all $j\notin
J_{n}$. If $\mathsf{W}$ is a wavelet transformation then its columns,
$w_{j}\in L^{2}$, are simply the wavelets themselves and we can use the bound

$|\left\langle
w_{j},\mathsf{A}^{*}\vec{\varphi}_{n}\right\rangle|=\left|\left\langle
w_{j},\mathds{1}_{\operatorname{supp}(w_{j})}\mathsf{A}^{*}\vec{\varphi}_{n}\right\rangle\right|\leq{\left\lVert\mathds{1}_{\operatorname{supp}(w_{j})}\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{2}}\leq{\left\lVert\mathds{1}_{\mathds{X}}\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{2}}$
(69)

for all $\mathds{X}\supset\operatorname{supp}(w_{j})$. In the case of the
Radon transform, we can compute the left-hand side explicitly for the finitely
many $j\in J_{n}$, but we wish to use the right-hand side in a structured way
to avoid computing the infinitely many $j\notin J_{n}$. To do this, we will
take a geometrical perspective on the construction of wavelets to view them in
a tree format.

##### Tree structure of wavelets

Finite elements are constructed with a mesh which provided a useful tool for
adaptive refinement in Section 6.3.3. For wavelets, we will associate a tree
with every discretisation and the leaves of the tree correspond to a mesh.
This perspective comes from the multi-resolution interpretation of wavelets.
An example is seen in Fig. 6 for 1D Haar wavelets,
$w_{j,k}(x)=\sqrt{2}^{k}\psi(2^{k}x-j)$ where
$\psi=\mathds{1}_{[0,1)}-\mathds{1}_{[-1,0)}$.

In higher dimensions, the only two things which change are the number of
children ($2^{d}$ for non-leaves) and at each node you store the coefficients
of $2^{d}-1$ wavelets. The support on each node is still a disjoint partition
of unity consisting of regular cubes of side length $2^{-k}$ at level $k$. The
only change in our own implementation is to translate the support to
$[-\tfrac{1}{2},\tfrac{1}{2}]^{2}$. We briefly remark that the tree
structuring of wavelets is not novel and appears more frequently in the
Bayesian inverse problems literature Castillo2019 ; Kekkonen2021 .

##### Continuous gradient estimate

In Section 7.1 we used the continuous gap as a measure for convergence, for
wavelets we will use the continuous subdifferential. With the tree structure
we can easily adapt the results of Section 6.3 to estimate subdifferentials
(or function gaps). In particular,

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$
$\displaystyle=\max\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*},\max_{j\notin
J_{n}}|\left\langle
w_{j},\mathsf{A}^{*}\vec{\varphi}_{n}\right\rangle|-\mu\right)$ (70)
$\displaystyle\leq\max\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*},\max_{j\in\operatorname{leaf}(J_{n})}{\left\lVert\mathds{1}_{\operatorname{supp}(w_{j})}\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{2}}-\mu\right).$
(71)

##### Numerical results

We consider two phantoms where the ground-truth is either a binary disc or the
Shepp-Logan phantom. Both examples are corrupted with $2\text{\,}\mathrm{\char
37\relax}$ Laplace distributed noise. This is visualised in Fig. 9. All
optimisations shown are spatially adaptive using Haar wavelets and initialised
with $\mathds{U}^{0}=\\{x\mapsto c\operatorname{\;s.t.\;}c\in\mathds{R}\\}$.
The gradient metric shown throughout is the $\ell^{\infty}$ norm. Motivated by
(71), the spatial adaptivity is chosen to refine nodes
$j\in\operatorname{leaf}(J_{n})$ to ensure that

${\left\lVert\mathds{1}_{\operatorname{supp}(w_{j})}\mathsf{A}^{*}\vec{\varphi}_{n}\right\rVert}_{L^{2}}-\mu\leq
10{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\partial_{n}\operatorname{E}(\IfEqCase{0}{{0}{u}{1}{v}{2}{w}}[u0]_{n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$

for all $j$ and $n$ (i.e. so that the continuous gradient is less than 10
times the discrete gradient). We do not expect wavelet regularisation to have
state-of-the-art performance in the examples of Fig. 9. What they demonstrate
is the preference Haar wavelets have to align large discontinuities with a
coarse grid, even when the discretisation is allowed to be as fine as
necessary. There is an average of $2\cdot 10^{6}$ wavelet coefficients in each
discretised reconstruction, although the higher frequencies have much smaller
intensities. In limited data scenarios, wavelet regularisation automatically
selects a local “resolution” which reflects the quality of data. Particularly
in the Shepp-Logan reconstruction, we see that the outer ring is detected with
a finer precision than the dark interior ellipses.

The first numerical results shown in Fig. 9 compare the same adaptive FISTA
variants as shown in Fig. 3. In these examples we see that the greedy FISTA
and the $a=20$ algorithms achieve almost linear convergence while $a=2$ is
significantly slower. Interestingly, in both examples the $a=20$ variant uses
half as many wavelets as the Greedy variant, and therefore converges slightly
faster in time.

Figure 7: Phantoms, data and reconstructions for wavelet-sparse tomography
optimisation. Both examples are corrupted with $2\text{\,}\mathrm{\char
37\relax}$ Laplace distributed noise.

Figure 8: Convergence of different implementations of Algorithm 1 with an
unlimited number of pixels for sparse wavelet optimisation.

Figure 9: Example images from STORM dataset.

###  7.3 2D continuous LASSO

Our final application is a super-resolution/de-blurring inverse problem from
biological microscopy. In mathematical terms, the observed data is a large
number of sparse images which are corrupted by blurring and a large amount of
noise, examples are seen in Fig. 9. The task is to compute the centres of the
spikes of signal in each image and then re-combine into a single super-
resolved image, as in Fig. 11. This technique is referred to as _Single
Molecule Localisation Microscopy_ (SMLM), of which we consider the specific
example of _Stochastic Optical Reconstruction Microscopy_ (STORM). Readers are
directed to the references Sage2015 ; Sage2019 ; Schermelleh2019  for further
details. The LASSO formulation
($\operatorname{f}(\cdot)=\frac{1}{2}{\left\lVert\cdot\right\rVert}_{\ell^{2}}^{2}$)
has previously been shown to be effective in the context of STORM Huang2017 ;
Denoyelle2019 .

Here we use a simulated dataset provided as part of the 2016 SMLM
challenge222http://bigwww.epfl.ch/smlm/challenge2016/datasets/MT4.N2.HD/Data/data.html
for benchmarking software in this application. The corresponding LASSO
formulation is