# A Multiresolution approach to solve
large-scale optimization problems

Rosa Donat, Sergio López-Ureña∗ Universitat de València, 50th Doctor Moliner
street, Burjassot (Spain<EMAIL_ADDRESS>sergio.lopez-urena@uv.es). The authors
have been supported by grant MTM2017-83942 funded by Spanish MINEC and by
grant PID2020-117211GB-I00 funded by MCIN/AEI/10.13039/501100011033.

###### Abstract

General purpose optimization techniques can be used to solve many problems in
engineering computations, although their cost is often prohibitive when the
number of degrees of freedom is very large. We describe a multilevel approach
to speed up the computation of the solution of a large-scale optimization
problem by a given optimization technique. By embedding the problem within
Harten’s Multiresolution Framework (MRF), we set up a procedure that leads to
the desired solution, after the computation of a finite sequence of sub-
optimal solutions, which solve auxiliary optimization problems involving a
smaller number of variables. For convex optimization problems having smooth
solutions, we prove that the distance between the optimal solution and each
sub-optimal approximation is related to the accuracy of the interpolation
technique used within the MRF and analyze its relation with the performance of
the proposed algorithm. Several numerical experiments confirm that our
technique provides a computationally efficient strategy that allows the end
user to treat both the optimizer and the objective function as black boxes
throughout the optimization process.

##  1 Introduction

Optimization problems of the type

$\textrm{Find }z_{\min}\in\mathbb{R}^{N}\textrm{ such that
}F(z_{\min})=\min_{z\in\mathbb{R}^{N}}F(z),\qquad N>>1,$ (1)

arise frequently in applications, often in connection with the context of
calculus of variations, PDE constrained optimization, optimal control [8, 9,
25] or image processing [11, 10].

In many occasions, specially in engineering computations, it is known, a
priori, that the problem has a unique solution which can be computed (from a
given initial data) with an appropriate ‘off-the-shelf’ (black-box type)
optimization technique which, however, tends to be very time consuming.
Assuming that problem eq. 1 can be solved by using a given optimization
technique, ${\cal D}$, we describe a strategy to reduce the computational time
required by the direct application of $\mathcal{D}$ on a space with $N$
degrees of freedom. Our technique is based on embedding problem eq. 1 in
Harten’s Multiresolution Framework (MRF henceforth) [27], which can be loosely
described as a set of mathematical tools to obtain multilevel representations
of discrete data sets, similar to those obtained in the well known wavelet
framework.

Harten’s MRF has been successfully applied for different purposes and in
several scenarios (see e.g. [18, 1, 5, 14, 13, 7]) but, to the best of our
knowledge, the applications carried out in [19, 28] constitute the first
attempt to use Harten’s MRF in connection with constrained/unconstrained
optimization. The aim of this paper is to provide a complete mathematical
description of the technique used in [19, 28] to solve large-scale
optimization problems of the type eq. 1, together with some theoretical
results on its properties and performance, including several numerical
experiments that confirm our theoretical analysis. Our technique will be
referred throughout the paper as MR/OPT (for Multiresolution for
Optimization).

###  1.1 Related Work

The use of multilevel techniques to solve large-scale optimization problems
has been a common trend in the literature since the mid 80’s. Largely fueled
by the success of multigrid methods as efficient solvers for discretized
partial differential equations (PDE) [8, 9, 25]) and the fact that many large-
scale optimization problems arise from the discretization of a PDE [30, 15,
12, 22, 23, 24], there are nowadays many multilevel algorithms which exploit
the fact that, if eq. 1 comes from an underlying infinite dimensional problem,
it may be described at several discretization levels.

The idea behind many of these multilevel optimization techniques is to combine
an iterative (or descent) method with a correction which requires the solution
of a reduced, lower dimensional, auxiliary problem which is cheaper to solve.
For example, in the MG/OPT framework introduced by Nash in [30], the solutions
of low-dimensional problems (coarse resolution problems) are used to compute
search directions for higher dimensional problems in a multigrid fashion, an
idea that has been considered by Frandi and Papini in a series of papers [22,
23] where derivative-free multilevel optimization techniques are considered.
Other techniques, such as the trust-region multilevel optimization also make
use of auxiliary, lower dimensional, optimization problems (or surrogate
models [22, 23]) which make an explicit use of the family of discrete
versions, at different resolution levels, associated to the infinite-
dimensional/continuous problem.

The present paper describes a different strategy to solve eq. 1, which is
close, in spirit, to the so-called ‘cost-effective’ multiscale strategy used
in [13], whose objective is to reduce the cost of using a time-consuming
numerical scheme on a fine-mesh simulation involving hyperbolic PDEs. Our aim
is not to design a new optimization algorithm, but rather to speed up the
computation carried out by a convenient (but possibly very time-consuming)
optimization technique, starting at a given initial guess provided by the end-
user.

In the MR/OPT algorithm, the function $F$ in eq. 1 and the optimization tool
are treated as black boxes, with the only requirement that ${\cal D}$ is
appropriate to solve eq. 1 (except for the possibly very large computational
expense). Its design involves the computation of a sequence of iterates, which
solve auxiliary optimization problems in $\mathbb{R}^{m}$ ($m<N$) which do not
make use of any discretization of the underlying infinite-dimensional problem
at lower resolution levels (in contrast with [15, 22, 23, 24]). Somewhat
related multilevel optimization strategies can be found in e.g. [19, 26, 29]
(and references therein) in the context of aerodynamic shape optimization,
where it is now recognized that a shape optimization problem can be solved as
a sequence of optimization steps that evolve from a basic parametrization of
the shape to more refined parametrizations. These techniques are capable of
implementing geometry changes while maintaining the required constraints and
obtaining an improvement in the efficiency, accuracy, and robustness of widely
used optimization procedures. As in the MR/OPT strategy, they are not based on
discretizing the underlying infinite dimensional problem at several accuracy
levels.

The paper is organized as follows: In section 2 we briefly recall the main
ingredients of Harten’s MRF required for the description of the MR/OPT
strategy, which is carefully described in section 3. In section 2.1, we
describe the interpolatory MR setting, including the specific prediction
schemes used in the numerical examples of section 4, chosen to illustrate the
performance and limitations of our technique. We finish with some conclusions
and perspectives in section 5.

##  2 Harten’s Framework For Multiresolution

In this section we shall recall the essential ingredients of Harten’s MRF
which are required for the description of our algorithm (for a complete
description of Harten’s MRF see e.g. [27, 2, 3]).

A distinctive characteristic of Harten’s approach is that the discrete data at
different resolution levels in a specific MRF are assumed to come from a
particular discretization procedure, acting on an appropriate functional
space, which is ‘naturally’ linked to a hierarchy of nested meshes111$\\{{\cal
G}^{k}\\}_{k}$ is nested if ${\cal G}^{k}\subset{\cal G}^{k+1}$ $({\cal
G}^{k})_{k=0}^{L}$ covering a fixed spatial domain.

Assuming that $N_{k}$ is the dimension of the discrete data attached to ${\cal
G}^{k}$, where $N_{L}>N_{L-1}>\cdots>N_{0}$, a multiresolution (MR)
representation of a discrete data set, $z^{L}\in\mathbb{R}^{N_{L}}$, is an
equivalent representation (in $\mathbb{R}^{N_{L}}$) which consists in a coarse
representation of this data set, $z^{0}\in\mathbb{R}^{N_{0}}$, together with a
sequence of details, $d^{\ell}\in\mathbb{R}^{N_{\ell+1}-N_{\ell}}$,
$\ell=0,\cdots,L-1$, that can be understood as ‘essential’ (non-redundant)
difference in information between consecutive resolution levels.

In a nutshell, Harten’s MRF can be described as follows: For each pair of
consecutive resolution levels, we are given a decimation operator
$D_{k}^{k-1}$, which extracts $k-1$ level information from the discrete data
at level $k$, and a prediction operator $P_{k-1}^{k}$ which produces $k$ level
data from $k-1$ level information. $D_{k}^{k-1}$ is a linear operator, fully
determined by the underlying discretization procedure that defines the
framework. $P_{k-1}^{k}$ may be linear or nonlinear222This is a major
difference between Harten’s MRF and the classical wavelet theory.. These two
operators are required to satisfy a consistency relation at all resolution
levels (Iℓ is the identity operator on $\mathbb{R}^{N_{\ell}}$):

$D_{k}^{k-1}P_{k-1}^{k}=\textrm{I}_{k-1},$ (2)

(but $P_{k-1}^{k}D_{k}^{k-1}\neq\textrm{I}_{k}$). Notice that, for
$z^{k}\in\mathbb{R}^{N_{k}}$ we can write

$z^{k}=P_{k-1}^{k}D_{k}^{k-1}z^{k}+(\textrm{I}_{k}-P_{k-1}^{k}D_{k}^{k-1})z^{k}=P_{k-1}^{k}z^{k-1}+e^{k},$

hence each $z^{k}\in\mathbb{R}^{N_{k}}$ can be represented in terms of a
‘decimated’ version (at the next, coarser, resolution level) together with a
‘prediction error’, i.e.

$z^{k}\quad\equiv\quad\begin{cases}z^{k-1}&\\!\\!\\!:=D_{k}^{k-1}z^{k}\in\mathbb{R}^{N_{k-1}},\\\
e^{k}&\\!\\!\\!:=z^{k}-P_{k-1}^{k}z^{k-1}\in\mathbb{R}^{N_{k}}.\end{cases}$
(3)

The consistency relation eq. 2 implies that $D_{k}^{k-1}e^{k}=0$, so that the
representation of $z^{k}$ in terms of $(z^{k-1},e^{k})$ specified in eq. 3
(with $N_{k-1}+N_{k}$ components) is redundant. The details $d^{k-1}$
represent exactly the non-redundant information contained in $e^{k}$, so that
a one-to-one correspondence between $z^{k}$ and $(z^{k-1},d^{k-1})$ can be
established (hence $d^{k-1}\in\mathbb{R}^{N_{k}-N_{k-1}}$). By iterating this
procedure from level $L$ to level $1$, we obtain the full MR representation of
$z^{L}\in\mathbb{R}^{N_{L}}$:

$z^{L}\leftrightarrow(z^{L-1},d^{L-1})\leftrightarrow(z^{L-2},d^{L-2},d^{L-1})\leftrightarrow\cdots\leftrightarrow(z^{0},d^{0},\cdots,d^{L-1}).$
(4)

In eq. 4, $z^{k}=D_{k+1}^{k}\cdots D^{L-1}_{L}z^{L}\in\mathbb{R}^{N_{k}}$,
$L>k\geq 0$, i.e. $z^{k}$ are coarse versions of $z^{L}$ at the different
resolution levels, obtained by successive decimation. Moreover, the details
$d^{k-1}$ are related to the prediction error
$e^{k}=(I-P_{k-1}^{k}D_{k}^{k-1})z^{k}$, and provide information on the local
smoothness of the function which is associated to the discrete data $z^{L}$ by
the discretization operator which defines the framework (see also section
2.1).

If we denote by $M_{m,j}$ ($m<j$) the multiresolution transform operator that
acts between the resolution levels $m$ and $j$, then, for $0\leq\ell\leq k\leq
L$ and $\forall z^{L}\in\mathbb{R}^{N_{L}}$,

$M_{k}^{L}z^{L}=(z^{k},d^{k},\cdots d^{L-1})\quad\rightarrow\quad
M_{\ell}^{L}z^{L}=(M_{\ell}^{k}z^{k},d^{k},\cdots d^{L-1}),$ (5)

or, equivalently,

$M_{\ell,L}z^{L}=(z^{\ell},d^{\ell},\cdots,d^{k-1},d^{k},\cdots,d^{L-1})\,\equiv\,\begin{cases}M_{k,L}\,z^{L}=(z^{k},d^{k},\cdots
d^{L-1}),\\\ M_{\ell,k}z^{k}=(z^{\ell},d^{\ell},\cdots,d^{k-1}).\end{cases}$
(6)

We shall use these relations in various places throughout the paper.

An important property for the practical use of MR transformations is their
stability with respect to perturbations. In this paper, we shall assume that
the prediction operators
$P_{k-1}^{k}:\mathbb{R}^{N_{k-1}}\rightarrow\mathbb{R}^{N_{k}}$ are linear
(they can be represented as $\mathbb{R}^{N_{k}\times N_{k-1}}$ matrices) which
implies in turn that the multiscale transformations
$M_{m,j}:\mathbb{R}^{N_{j}}\rightarrow\mathbb{R}^{N_{j}}$ are also a linear
operators ($\mathbb{R}^{N_{j}\times N_{j}}$ matrices). In this case, the
stability of the associated MRF can be stated as follows

$||M_{m,j}||_{\infty}\leq C,\qquad||M_{m,j}^{-1}||_{\infty}\leq\tilde{C},$ (7)

with $C,\tilde{C}$ independent of $m,j$ or $L$, for $0\leq m\leq j\leq L$, and
$||\cdot||_{\infty}$ the usual operator norm.

Within Harten’s framework, stability is related to the convergence of the
process that results from the recursive application of the prediction
operators, which involves a refinement process akin to those found in the
theory of subdivision refinement333The convergence theory of these processes
is outside the scope of this paper. The reader is referred to [20] for the
concept of convergence in subdivision refinement.. The prediction operators
considered in this paper satisfy a property which guarantees the stability of
the associated MR transformations. This property, which we shall denote as
Property S, is defined below.

######  Definition 1 (Property S).

Let ${\cal P}=\\{P_{k}^{k+1}\\}_{k=0}^{\infty}$ be a sequence of linear
prediction operators. We will say that ${\cal P}$ satisfies property S if
there exist $d_{1},d_{2}>0$, such that for any $0\leq\ell<k$

$d_{1}||z^{\ell}||_{\infty}\leq||P_{\ell}^{k}z^{\ell}||_{\infty}\leq
d_{2}||z^{\ell}||_{\infty},\qquad\forall z^{\ell}\in\mathbb{R}^{N_{\ell}},$
(8)

where $P_{\ell}^{k}:=P_{k-1}^{k}\cdots P_{\ell+1}^{\ell+2}P_{\ell}^{\ell+1}$
is the operator that computes the recursive prediction process between levels
$l<k$.

######  Remark..

We shall assume that $P_{k}^{k}=I_{k}$. Notice that

$z^{L}=P_{k}^{L}z^{k}\quad\Leftrightarrow\quad
M_{k,L}z^{L}=(z^{k},0,\ldots,0).$ (9)

###  2.1 Harten’s Interpolatory MRF

In this paper, we shall use Harten’s interpolatory MR framework, which is
characterized by the assumption that the discrete data correspond to the
point-evaluations of a function, defined on the spatial domain underlying the
nested mesh structure.

We briefly describe next the one-dimensional case: Let us consider the
sequence of nested meshes on the interval $[0,1]$ obtained by the dyadic
refinement of a uniform mesh ${\cal G}^{0}$, with grid spacing
$h_{0}=1/J_{0}$. Thus, ${\cal G}^{k}=\\{ih_{k}\\}_{i=0}^{J_{k}}$ with
$h_{k}=h_{0}/2^{k}=1/J_{k}$ and the $k$-th level information consists in
$N_{k}=J_{k}+1$ discrete values, for each $k$, $0\leq k\leq L$.

The operators that define the transfer of information between two consecutive
resolution levels ($k$ fine, $k-1$ coarse) in the 1D interpolatory MRF are as
follows:
$\bullet$ Decimation by downsampling: For $z^{k}\in\mathbb{R}^{N_{k}}$,

$z^{k-1}_{i}=(D_{k}^{k-1}z^{k})_{i}=z^{k}_{2i},\quad 0\leq i\leq J_{k-1}.$
(10)

Notice that this decimation operator corresponds to the interpretation of the
discrete data as the point-values of an underlying function defined on the
unit interval $[0,1]$ at each resolution level. Obviously
$z^{k-1}\in\mathbb{R}^{N_{k-1}}$ is a coarse version of the data $z^{k}$,
which corresponds to the point-values of the function on the grid ${\cal
G}^{k-1}$.
$\bullet$ Prediction by interpolation. Because of eq. 10, the consistency
relation eq. 2 becomes

$z^{k-1}_{i}=(D_{k}^{k-1}P_{k-1}^{k}z^{k-1})_{i}=(P_{k-1}^{k}z^{k-1})_{2i}.$

This interpolation condition implies that the prediction errors at even points
are zero, thus, only the predicted values at odd grid points need to be
computed and stored. Then, the one-to-one correspondence between the discrete
sets $z^{k}$ and $(z^{k-1},d^{k-1})$, i.e. the two-level multiresolution
transformation $M_{k-1,k}$ and its inverse, are

$M_{k-1,k}z^{k}=(z^{k-1},d^{k-1})\
\leftrightarrow\begin{cases}z^{k-1}_{i}\\!\\!\\!\\!\\!\\!&=z^{k}_{2i},\,\,0\leq
i\leq J_{k-1},\\\
d_{i}^{k-1}\\!\\!\\!\\!\\!\\!&=z^{k}_{2i-1}-(P_{k-1}^{k}z^{k-1})_{2i-1},\,\,1\leq
i\leq J_{k-1},\end{cases}$ (11) $M_{k-1,k}^{-1}(z^{k-1},d^{k})=z^{k}\
\leftrightarrow\begin{cases}z^{k}_{2i}\\!\\!\\!\\!\\!\\!&=z^{k-1}_{i},\,\,0\leq
i\leq J_{k-1},\\\
z^{k}_{2i-1}\\!\\!\\!\\!\\!\\!&=d_{i}^{k-1}+(P_{k-1}^{k}z^{k-1})_{2i-1},\,\,1\leq
i\leq J_{k-1}.\end{cases}$ (12)

######  Remark..

From eq. 11, since $z^{k-1}=D_{k}^{k-1}z^{k}$, $z^{k}=D_{k+1}^{k}\cdots
D_{L}^{L-1}z^{L}$, we can also write $d_{i}^{k-1}=d_{i}^{k-1}(z^{L})$, which
clearly express that the detail coefficients are determined by the properties
of the discrete data at the finest resolution level.

In this paper, we shall consider prediction operators based on certain
piecewise polynomial interpolatory techniques of fixed polynomial degree, $n$,
odd. We refer the reader to [3, 16] for specific details, and simply provide
the prediction rules for $n=1,3,5$

$\begin{cases}(P_{k-1}^{k}z^{k-1})_{2i-1}\\!\\!\\!&=\sum_{\ell=1}^{(n+1)/2}\beta_{\ell}(z^{k-1}_{i+\ell}+z^{k-1}_{i+1-\ell}),\quad\frac{n+1}{2}\leq
i\leq J_{k-1}-\frac{n-1}{2},\\\
(P_{k-1}^{k}z^{k-1})_{2i}\\!\\!\\!&=z^{k-1}_{i},\quad 0\leq i\leq
J_{k-1},\end{cases}$ (13) $\begin{cases}n=1\Rightarrow&\beta_{1}=1/2,\\\
n=3\Rightarrow&\beta_{1}=9/16,\,\beta_{2}=-1/16,\\\
n=5\Rightarrow&\beta_{1}=150/256,\,\beta_{2}=-25/256,\,\beta_{3}=3/256,\end{cases}$
(14)

with the following special boundary rules at the left boundary of $[0,1]$ for
$n=3$,

$(P_{k-1}^{k}z^{k-1})_{1}=\frac{5}{16}z^{k-1}_{0}+\frac{15}{16}z^{k-1}_{1}-\frac{5}{16}z^{k-1}_{2}+\frac{1}{16}z^{k-1}_{3},$
(15)

and $n=5$,

$\\!\begin{cases}\\!\\!(P_{k-1}^{k}z^{k-1})_{1}\\!\\!\\!\\!\\!\\!&=\frac{63}{256}z^{k-1}_{0}\\!+\frac{315}{256}z^{k-1}_{1}\\!-\\!\frac{105}{128}z^{k-1}_{2}\\!+\\!\frac{63}{128}z^{k-1}_{3}\\!-\\!\frac{45}{256}z^{k-1}_{4}\\!+\\!\frac{7}{256}z^{k-1}_{5},\\\\[6.0pt]
\\!\\!(P_{k-1}^{k}z^{k-1})_{3}\\!\\!\\!\\!\\!\\!&=-\frac{7}{256}z^{k-1}_{0}\\!+\\!\frac{105}{256}z^{k-1}_{1}\\!+\\!\frac{105}{128}z^{k-1}_{2}\\!-\\!\frac{35}{128}z^{k-1}_{3}\\!+\\!\frac{21}{256}z^{k-1}_{4}\\!-\\!\frac{3}{256}z^{k-1}_{5}.\end{cases}$
(16)

The corresponding rules at the right end of $[0,1]$ can be obtained by
symmetry. We prove in the Appendix that these prediction operators satisfy
property S with $d_{1}=1$, hence, the associated MR transformations are
stable, i.e. satisfy eq. 7.

######  Remark..

The prediction operators defined by eq. 13 are related to the well known
(stationary444Stationary subdivision rules are independent of the refinement
level.) subdivision rules of the interpolatory Deslauries-Dubuc subdivision
schemes [20, 16] for infinite sequences. To deal with bounded domains, we have
to take care of the, possibly special, prediction rules close to the
boundaries, where the required centered stencils might leave the domain. For
$n=1$, i.e. linear interpolation, only the boundary values
$z^{k-1}_{0},z^{k-1}_{N_{k-1}}$ are required to apply the prediction operator
specified in eq. 13 at each resolution level, hence no special rules are
required close to the boundaries of the domain. However, for $n=3$ (or $n=5$)
interpolation polynomials of degree 3 (or 5) are involved, hence special rules
are required to predict $k$-data from data from $(k-1)$-th grid values at the
positions $i=1$ and $i=J_{k-1}-1$ (or $i=1,2$, $i=J_{k-1}-1,J_{k-1}-2$ for
$n=5$). To obtain eqs. 15 and 16, we have used the stencil of the $n+1$ points
closest to the boundary, in order to obtain the required predicted value [3].
Hence, all polynomial pieces in the interpolatory reconstruction process have
the same degree (and the same approximation order).

The decay of the scale coefficients with respect to the resolution level under
smoothness assumptions, is a well known fact that is very simple to explain in
the interpolatory MRF. Since the details are interpolation errors, their size
(and behavior across scales) is determined by the approximation order of the
interpolatory technique used to define the prediction operators. The rules in
eqs. 13, 14, 15, and 16, involve polynomials of degree $n$, hence

$d^{k-1}_{i}={\cal O}(h_{k}^{n+1}).$ (17)

when the data around the spatial location corresponding to the indices
$(2i-1,k)$ correspond to the point values of a sufficiently smooth function.

In this paper, we shall make use of the 1D interpolatory MRF just described
and its 2D version obtained by using a tensor-product approach [6]. When using
the 1D prediction operators specified above in the 2D tensor-product context,
the prediction operators have the same order of approximation, with respect to
the uniform mesh spacing of the (tensor product) hierarchical mesh structure,
as in the 1D setting. Hence, the decay of the scale coefficients in the 2D
tensor-product interpolatory MRF is also given by eq. 17.

##  3  MR/OPT: A multiresolution approach for the efficient solution of large-
scale optimization problems

We shall describe next a multilevel strategy to compute the solution of eq. 1,
assuming that the end-user provides an initial guess $\bar{z}$ and a preferred
optimization tool ${\cal D}$ so that

$\quad z_{\min}=\mathcal{D}(\bar{z},F,N)=\operatorname*{arg\
min}\\{F(z),\,\,\,z\in\mathbb{R}^{N}\\}.$ (18)

The notation above is meant to emphasize the fact that the optimization
problem eq. 1, with $F$ as the objective function on a discrete space with $N$
degrees of freedom, can be solved using $\cal D$ as optimization tool and
$\bar{z}$ as the initial guess.

Considering $N=N_{L}$ as the finest resolution level in a sequence of $L+1$
levels of refinement corresponding to discrete spaces $\mathbb{R}^{N_{k}}$ (
$N_{0}<N_{1}<\cdots<N_{L}=N$, $L>0$) within a MRF which is appropriate for the
given problem, the MR/OPT strategy produces a sequence of ‘sub-optimal’
solutions, $\\{z^{L,k}\\}_{k=0}^{L+1}$, satisfying

$z^{L,0}:=\bar{z},\quad F(\bar{z})=F(z^{L,0})\geq F(z^{L,1})\geq\cdots\geq
F(z^{L,L+1}),\quad z^{L,L+1}=z_{\min}.$ (19)

The $k$-th level sub-optimal solution satisfies

$F(z^{L,k+1})=\min\\{F(z),\,z\in\Xi_{k}\\}$

where $\Xi_{k}$ is a subset of $\mathbb{R}^{N}$, whose construction depends on
the MRF being used and on the sequence of prediction operators ${\cal P}$.

In this paper we always assume that the prediction operators in ${\cal P}$ are
linear and satisfy property S. Hence the MR transformations used are also
linear operators and eq. 7 is satisfied.

At the coarsest resolution level ($k=0$) we set $z^{L,0}:=\bar{z}$, the
initial guess provided by the user. Its MR representation,
$M_{0,L}z^{L,0}=M_{0,L}\bar{z}=:(\bar{z}^{0},\bar{d}^{0},\ldots,\bar{d}^{L-1})$,
serves to define the set $\Xi_{0}\subset\mathbb{R}^{N_{L}}$ and the associated
$0$-th-level objective function as follows

$\displaystyle\Xi_{0}:=$
$\displaystyle\\{M_{0,L}^{-1}(\bar{z}^{0}+\varepsilon^{0},\bar{d}^{0},\ldots,\bar{d}^{L-1}),\,\,\varepsilon^{0}\in
R^{N_{0}}\\},$ $\displaystyle F_{0}(\varepsilon^{0}):=$ $\displaystyle
F(M_{0,L}^{-1}(\bar{z}^{0}+\varepsilon^{0},\bar{d}^{0},\bar{d}^{1},\ldots,\bar{d}^{L-1})),$
$\displaystyle\varepsilon^{0}\in\mathbb{R}^{N_{0}}.$

Since the MR operators are linear, taking into account eq. 9, we can write

$M_{0,L}^{-1}(\bar{z}^{0}+\varepsilon^{0},\bar{d}^{0},\bar{d}^{1},\ldots,\bar{d}^{L-1})=\bar{z}+P_{0}^{L}\varepsilon^{0}.$

Hence $\Xi_{0}$ is an affine space with $N_{0}$ degrees of freedom, which is
defined by considering only perturbations (of the MR representation) of the
initial data at the coarsest resolution level. These perturbations are moved
to $\mathbb{R}^{N_{L}}$ by successive prediction.

We then compute (assuming it is feasible555We shall see that this is indeed
the case for convex optimization problems.)

$\varepsilon^{0}_{*}:=\operatorname*{arg\
min}\\{F_{0}(\varepsilon^{0}),\quad\varepsilon^{0}\in\mathbb{R}^{N_{0}}\\}={\cal
D}(0,F_{0},{N_{0}}),$ (20)

and define the sub-optimal solution at level $0$ as

$z^{L,1}:=\bar{z}+P_{0}^{L}\varepsilon^{0}_{*}=z^{L,0}+P_{0}^{L}\varepsilon^{0}_{*}.$
(21)

Obviously

$z^{L,1}=\operatorname*{arg\ min}\\{F(z),\,\,z\in\Xi_{0}\\},\quad
F(z^{L,1})=F_{0}(\varepsilon^{0}_{*})\leq F_{0}(0)=F(z^{L,0})=F(\bar{z}).$

Notice that the computation of $z^{L,1}$ can be carried out involving only
direct calls to the objective function $F$, since
$F_{0}(\varepsilon^{0})=F(\bar{z}+P_{0}^{L}\varepsilon^{0})$. Notice that the
user only needs to provide $F$, and the MRF takes care of the definition of
the objective function $F_{0}$. In addition, considering
$\varepsilon^{0}=0\in\mathbb{R}^{N_{0}}$ as the initial guess to solve eq. 20
(the $0$-th level auxiliary optimization problem) corresponds to considering
$z^{L,0}=\bar{z}$ as the ‘corresponding’ initial guess in the affine space
$\Xi_{0}\subset\mathbb{R}^{N}$. Notice also that

$M_{0,L}z^{L,1}=M_{0,L}\bar{z}+M_{0,L}P_{0}^{L}\varepsilon^{0}_{*}=(\bar{z}^{0}+\varepsilon^{0}_{*},\bar{d}^{0},\cdots,\bar{d}^{L-1}),$

i.e. the difference between $z^{L,0}$ and $z^{L,1}$ occurs only at the
coarsest resolution level, while the ‘details’ in their MR representations
coincide at all higher resolution levels. Taking into account eq. 5, we can
write

$M_{1,L}z^{L,1}=(z^{1}_{*},\bar{d}^{1},\cdots,\bar{d}^{L}),\quad
M_{0,1}z^{1}_{*}=(\bar{z}_{0}+\varepsilon^{*}_{0},\bar{d}^{0})\,\equiv\,z^{1}_{*}=M_{0,1}^{-1}(\bar{z}_{0},\bar{d}^{0})+P_{0}^{1}\varepsilon^{0}_{*}$

so that we can repeat the process as if the level 1 was the coarsest
resolution level.

The general $k$-th step of the algorithm, $1\leq k\leq L$, which will provide
the sub-optimal solution $z^{L,k+1}$, can be described as follows: We assume
that we have obtained $\varepsilon_{*}^{m}\in\mathbb{R}^{N_{m}}$, $0\leq m\leq
k-1$ (the solution of the $m$-th level auxiliary optimization problem) and the
associated sub-optimal solutions, $\\{z^{L,m}\\}_{m=1}^{k}$, which satisfy the
following properties ($\bar{z}=z^{L,0}$)

$z^{L,m}=z^{L,m-1}+P_{m-1}^{L}\varepsilon_{*}^{m-1};\quad
M_{m,L}z^{L,m}=(z^{m}_{*},\bar{d}^{m},\ldots,\bar{d}^{L-1}),\quad 1\leq m\leq
k,\,\,$

where $z^{0}_{*}=\bar{z}^{0}$ and

$\begin{cases}z^{m}_{*}\\!\\!&=M_{m-1,m}^{-1}(z_{*}^{m-1},\bar{d}^{m-1})+P_{m-1}^{m}\varepsilon^{m-1}_{*},\\\
\bar{d}^{j}\\!\\!&=d^{j}(\bar{z}),\,\,m-1\leq j\leq L-1.\end{cases}$

Then, we proceed as follows:
1. Consider the $k$-th level space of approximation

$\Xi_{k}:=\\{M_{k,L}^{-1}(z^{k}_{*}+\varepsilon^{k},\bar{d}^{k},\cdots,\bar{d}^{L-1}),\,\,\varepsilon^{k}\in\mathbb{R}^{N_{k}}\\}=z^{L,k}+\\{P_{k}^{L}\varepsilon^{k},\,\,\varepsilon^{k}\in\mathbb{R}^{N_{k}}\\}$
(22)

which is an affine space in $\mathbb{R}^{N}$ with ${N_{k}}$ degrees of
freedom.

2. Define the $k$-th level objective function $F_{k}:\mathbb{R}^{N_{k}}\rightarrow\mathbb{R}$,

$F_{k}(\varepsilon^{k}):=F(M_{k,L}^{-1}(z^{k}_{*}+\varepsilon^{k},\bar{d}^{k},\ldots,\bar{d}^{L-1}))=F(z^{L,k}+P_{k}^{L}\varepsilon^{k}).$
(23)

Notice that only $F$ and ${\cal P}=\\{P_{k}^{k+1}\\}_{k}$ is required to
define $F_{k}$.

3. Assuming that it is feasible, compute the solution to the $k$-th level auxiliary optimization problem

$\varepsilon^{k}_{*}:=\operatorname*{arg\
min}\\{F_{k}(\varepsilon^{k}),\,\,\,\varepsilon^{k}\in\mathbb{R}^{N_{k}}\\}={\cal
D}(0,F_{k},N_{k}).$ (24)

4\.  Define

$\quad z^{L,k+1}:=z^{L,k}+P_{k}^{L}\varepsilon^{k}_{*}.$ (25)

Then, $z^{L,k+1}=\operatorname*{arg\ min}\\{F(z),\,\,z\in\Xi_{k}\\}$,
$F(z^{L,k+1})=F_{k}(\varepsilon^{k}_{*})\leq F_{k}(0)=F(z^{L,k})$ and

$M_{k,L}z^{L,k+1}=M_{k,L}z^{L,k}+M_{k,L}P_{k}^{L}\varepsilon^{k}_{*}=(z^{k}_{*}+\varepsilon^{k}_{*},\bar{d}^{k},\ldots,\bar{d}^{L-1})$
(26)

hence

$\displaystyle d^{i}(z^{L,k+1})=\bar{d}^{i}=d^{i}(\bar{z}),\ k\leq i\leq L-1,$
(27) $\displaystyle
M_{k+1,L}z^{L,k+1}=(z^{k+1}_{*},\bar{d}^{k+1},\cdots,\bar{d}^{L-1}),\quad
z^{k+1}_{*}:=M_{k,k+1}^{-1}(z^{k}_{*},\bar{d}^{k})+P_{k}^{k+1}\varepsilon^{k}_{*}.$
(28)

We can easily prove the following properties of the affine spaces $\Xi_{k}$.

######  Lemma 1.

The spaces $\Xi_{k}$ in eq. 22 satisfy

  1. 1.

$\Xi_{k}\subset\Xi_{k+1},\quad 0\leq k<L.$

  2. 2.

$\Xi_{k}=\\{z^{L,\ell}+P_{k}^{L}\varepsilon^{k}\ :\
\varepsilon^{k}\in\mathbb{R}^{N_{k}}\\},\quad\forall\ell,0\leq\ell\leq k\leq
L.$

###### Proof.

The $\Xi_{k}$ are affine spaces, hence to prove item 1 it is enough to check
that $P_{k}^{L}(\mathbb{R}^{N_{k}})\subset P_{k+1}^{L}(\mathbb{R}^{N_{k+1}})$
and $z^{L,k}\in\Xi_{k+1}$. For this, we simply notice that
$\forall\varepsilon_{k}\in\mathbb{R}^{N_{k}}$

$P_{k}^{L}\varepsilon^{k}=P_{k+1}^{L}P_{k}^{k+1}\varepsilon^{k},\quad
z^{L,k}=z^{L,k+1}-P_{k+1}^{L}P_{k}^{k+1}\varepsilon^{k}_{*},\quad
P_{k}^{k+1}\varepsilon^{k}\in\mathbb{R}^{N_{k+1}}.$

Since $z^{L,\ell}\in\Xi_{\ell}\subset\Xi_{k}$, item 2 follows immediately. ∎

######  Remark..

By construction, $\Xi_{L}=\mathbb{R}^{N_{L}}=\mathbb{R}^{N}$, hence
$z^{L,L+1}=z_{\min}$ and the relations in eq. 19 are obviously satisfied.

The MR/OPT strategy substitutes the direct computation of the original (large-
scale) optimization problem, by the computation of the solutions of a sequence
of auxiliary optimization problems, each one of them associated to a level of
refinement in a multiresolution ladder. It is reasonable to assume, when these
auxiliary optimization problems admit a solution, it will be obtained quite
fast at low resolution levels, due to the reduced number of degrees of freedom
in the spaces where they are defined. In addition, it is expected that if the
distance between the sub-optimal solutions and the true solution gets smaller
while climbing up the MR ladder, then, even though each auxiliary
minimization problem in the MR ladder involves an increasing number of
variables, it is also expected that the improved initial guess (the previous
sub-optimal solution) chosen to carry out the optimization process will
improve the performance of the optimizer. In this case, the MR/OPT strategy
leads to an overall gain in functional evaluations.

In the following section, we shall see that, in some cases, it is possible to
ensure that the auxiliary optimization problems are of the same type as the
original one, hence their solution can be computed by the same optimization
technique, ${\cal D}$. In addition, we shall see that the previous
observations on the performance of the MR/OPT strategy can be theoretically
justified in some cases.

###  3.1 Some theoretical results on MR/OPT

We have assumed that the auxiliary optimization problems can be solved with
the same optimization tool as the full problem eq. 18. It is simple to see
that this is the case for quadratic minimization problems of the type

$\text{Find }z_{\min}\in\mathbb{R}^{N}\text{ such that
}F(z_{\min})=\min_{z\in\mathbb{R}^{N}}F(z),\quad
F(z)=\frac{1}{2}z^{T}Az-b^{T}z+c.$ (29)

with $A$ a symmetric and positive definite matrix.

######  Proposition 1.

Let us consider the quadratic minimization problem eq. 29, with $A$ symmetric
and positive definite. If the prediction operators in ${\cal P}$ satisfy eq.
8, then the auxiliary minimization problem at level $k$ eq. 24 is also a
quadratic minimization problem of the same type.

###### Proof.

Using the definition of $F$ in eq. 29, we have

$\displaystyle F_{k}(\varepsilon^{k})$
$\displaystyle=F(z^{L,k}+P_{k}^{L}\varepsilon^{k})=\frac{1}{2}(z^{L,k}+P_{k}^{L}\varepsilon^{k})^{T}A(z^{L,k}+P_{k}^{L}\varepsilon^{k})-b^{T}(z^{L,k}+P_{k}^{L}\varepsilon^{k})+c$
$\displaystyle=F(z^{L,k})+\frac{1}{2}(P_{k}^{L}\varepsilon^{k})^{T}AP_{k}^{L}\varepsilon^{k}+(z^{L,k})^{T}AP_{k}^{L}\varepsilon^{k}-b^{T}P_{k}^{L}\varepsilon^{k}$
$\displaystyle=\frac{1}{2}(\varepsilon^{k})^{T}(P_{k}^{L})^{T}AP_{k}^{L}\varepsilon^{k}-((P_{k}^{L})^{T}(b-Az^{L,k}))^{T}\varepsilon^{k}+F(z^{L,k}).$

Hence,

$F_{k}(\varepsilon^{k})=\frac{1}{2}(\varepsilon^{k})^{T}A_{k}\varepsilon^{k}-b_{k}^{T}\varepsilon^{k}+c_{k},\quad\begin{cases}A_{k}\\!\\!\\!\\!&=(P_{k}^{L})^{T}AP_{k}^{L}\in\mathbb{R}^{N_{k}\times
N_{k}},\\\
b_{k}\\!\\!\\!\\!&=(P_{k}^{L})^{T}(b-Az^{L,k})\in\mathbb{R}^{N_{k}},\\\
c_{k}&=F(z^{L,k})\in\mathbb{R}.\end{cases}$ (30)

Obviously $A_{k}$ is symmetric. If $A$ is positive definite, then $A_{k}$ is
also positive definite, since $P_{k}^{L}$ is injective, by eq. 8:

$\varepsilon^{k}\neq 0\quad\longleftrightarrow\quad
P_{k}^{L}\varepsilon^{k}\neq 0\quad\longleftrightarrow\quad
0<(P_{k}^{L}\varepsilon^{k})^{T}A(P_{k}^{L}\varepsilon^{k})=(\varepsilon^{k})^{T}A_{k}\varepsilon^{k}.$

∎

These quadratic optimization problems are special cases of a larger class of
objective functions for which the same property holds: The auxiliary
optimization problems are of the same type as the original one.

######  Proposition 2.

Let $F:\mathbb{R}^{N_{L}}\rightarrow\mathbb{R}$ be. If the prediction
operators in ${\cal P}$ satisfy property S, the auxiliary functions
$F_{k}:\mathbb{R}^{N_{k}}\rightarrow\mathbb{R}$ defined in eq. 23 satisfy:

  1. 1.

If $F$ is convex and/or $F\in{\cal C}^{2}(\mathbb{R}^{N_{L}},\mathbb{R})$,
then $F_{k}$ is also convex and/or $F_{k}\in{\cal
C}^{2}(\mathbb{R}^{N_{k}},\mathbb{R})$.

  2. 2.

If the hessian matrix $\nabla^{2}F(\xi^{L})$ is a positive definite matrix
$\forall\xi^{L}\in\mathbb{R}^{N_{L}}$, then $\nabla^{2}F_{k}(\xi^{k})$ is a
positive definite matrix $\forall\xi^{k}\in\mathbb{R}^{N_{k}}$.

  3. 3.

If $F$ is coercive, i.e.
$\lim_{||z^{L}||_{\infty}\rightarrow\infty}F(z^{L})=+\infty$, then $F_{k}$ is
coercive.

###### Proof.

Since $P_{k}^{L}$ is a linear operator, the proof of item 1 follows easily
from the expression of $F_{k}$ in eq. 23. To prove item 2, let
$\xi^{k}\in\mathbb{R}^{N_{k}}$ be and notice that $\forall\varepsilon^{k}\neq
0$,

$(\varepsilon^{k})^{T}\nabla^{2}F_{k}(\xi^{k})\varepsilon^{k}=(P_{k}^{L}\varepsilon^{k})^{T}\nabla^{2}F(z^{L,k}+P_{k}^{L}\xi^{k})(P_{k}^{L}\varepsilon^{k})>0,$

since $\varepsilon^{k}\neq 0\,\rightarrow\,P_{k}^{L}\varepsilon^{k}\neq 0$, by
eq. 8.

To prove item 3, let $0<\eta\in\mathbb{R}$. Since $F$ is coercive,
$\exists\delta=\delta(F,\eta)$ such that $F(z^{L})>\eta$, $\forall
z^{L}\in\mathbb{R}^{N_{L}}$: $||z^{L}||_{\infty}>\delta$. Define
$\delta^{\prime}:=d_{1}^{-1}(\delta+||z^{L,k}||_{\infty})$, with $d_{1}$ in
eq. 8. Then $\forall\varepsilon^{k}:$
$||\varepsilon^{k}||_{\infty}>\delta^{\prime}$ we have

$\|z^{L,k}+P_{k}^{L}\varepsilon^{k}\|_{\infty}\geq\|P_{k}^{L}\varepsilon^{k}\|_{\infty}-\|z^{L,k}\|_{\infty}\geq
d_{1}\|\varepsilon^{k}\|_{\infty}-\|z^{L,k}\|_{\infty}>d_{1}\delta^{\prime}-\|z^{L,k}\|=\delta,$

hence,

$F_{k}(\varepsilon^{k})=F(z^{L,k}+P_{k}^{L}\varepsilon^{k})>\eta,\quad\forall\varepsilon^{k}\,:\,||\varepsilon^{k}||_{\infty}>\delta^{\prime}.$

∎

The following results aim to estimate the difference between sub-optimal
solutions. For the sake of completeness and ease of reference, we include the
following technical lemma.

######  Lemma 2.

Let $F\in\mathcal{C}^{2}(\mathbb{R}^{N},\mathbb{R})$ be a convex function,
such that its Hessian matrix is always positive definite. Let
$A\subset\mathbb{R}^{N}$ be a compact convex set and let $x_{A}\in A$ be such
that $\nabla F(x_{A})=0$. Let $B\subset A$ be a compact convex set and
$x_{B}=\operatorname*{arg\ min}\\{F(x),\,x\in B\\}$. Then there exists
$c=c(A,F,\|\cdot\|)\geq 1$ such that

$\|x_{A}-x_{B}\|\leq c\ \text{dist}_{||\cdot||}(x_{A},B),$

where $\text{dist}_{||\cdot||}(x_{A},B)=\min\\{||x-x_{A}||,\,x\in B\\}$, for
any norm in $\mathbb{R}^{N}$.

###### Proof.

Since $\nabla F(x_{A})=0$, $B\subset A$ and $A$ is convex, for any $x\in B$,
$\exists\xi_{x}\in A$ such that

$F(x)-F(x_{A})=\frac{1}{2}(x-x_{A})^{T}\nabla^{2}F(\xi_{x})(x-x_{A}).$ (31)

The properties of $F$ and $A$ allow us to define $\rho_{\min}$, $\rho_{\max}$
as follows

$\displaystyle\rho_{\min}:=\min_{\xi\in A}\min_{x\neq
0}\frac{x^{T}\nabla^{2}F(\xi)x}{\|x\|^{2}}>0,\qquad\rho_{\max}:=\max_{\xi\in
A}\max_{x\neq 0}\frac{x^{T}\nabla^{2}F(\xi)x}{\|x\|^{2}}>0,$

so that

$\rho_{\min}\|x\|^{2}\leq
x^{T}\nabla^{2}F(\xi)x\leq\rho_{\max}\|x\|^{2},\qquad\forall\xi\in A,\ \forall
x\in\mathbb{R}^{N}.$

Let $w$ be an arbitrary but fixed point in $B$. Since $F(x_{B})\leq F(w)$,
applying eq. 31 with $x=w$ and $x=x_{B}$ and using the inequalities above, we
get

$\frac{\rho_{\min}}{2}\|x_{B}-x_{A}\|^{2}\leq F(x_{B})-F(x_{A})\leq
F(w)-F(x_{A})\leq\frac{\rho_{\max}}{2}\|w-x_{A}\|^{2}.$

Then,

$\|x_{A}-x_{B}\|\leq c\|x_{A}-w\|,\qquad c:=\sqrt{\rho_{\max}/\rho_{\min}}.$
(32)

The result follows from observing that the last inequality is true for every
$w\in B$. ∎

######  Proposition 3.

Let $F\in\mathcal{C}^{2}(\mathbb{R}^{N_{L}},\mathbb{R})$ be a convex coercive
function such that $\nabla^{2}F(\xi^{L})$ is a positive definite matrix
$\forall\xi^{L}\in\mathbb{R}^{N_{L}}$. If ${\cal
P}=\\{P_{k}^{k+1}\\}_{k=0}^{L-1}$ satisfies property S, then the sub-optimal
solutions $\\{z^{L,k}\\}_{k=1}^{L+1}$ in eq. 25 are well defined and unique.
In addition, for any compact and convex set $D$ which contains in its interior
the sequence $\\{z^{L,k}\\}_{k=0}^{L+1}$, there exists $C=C(F,D,{\cal P})>0$
such that for $\ell,k$: $0\leq\ell\leq k<L$,

$\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}\leq C\
\text{dist}_{\infty}(z^{L,k+1},\Xi_{\ell}\cap D).$ (33)

###### Proof.

By proposition 2, $\varepsilon^{k}_{*}$ in eq. 24 is uniquely defined, and,
thus, so is $z^{L,k+1}$, for each $k$, $0\leq k\leq L$.

Now, let $\ell,k$ be a pair of (fixed) indices, $0\leq\ell\leq k\leq L$. From
eq. 25, we can write

$\displaystyle z^{L,k+1}$
$\displaystyle=z^{L,k}+P_{k}^{L}\varepsilon^{k}_{*},\quad\varepsilon^{k}_{*}=\operatorname*{arg\
min}\\{F_{k}(\varepsilon^{k}),\,\varepsilon^{k}\in\mathbb{R}^{N_{k}}\\}.$
$\displaystyle z^{L,\ell+1}$
$\displaystyle=z^{L,\ell}+P_{\ell}^{L}\varepsilon^{\ell}_{*},\quad\varepsilon^{\ell}_{*}=\operatorname*{arg\
min}\\{F_{\ell}(\varepsilon^{\ell}),\,\varepsilon^{\ell}\in\mathbb{R}^{N_{\ell}}\\}.$

Thus, since $P_{\ell}^{L}=P_{k}^{L}P_{\ell}^{k}$,

$z^{L,k+1}-z^{L,\ell+1}=z^{L,k}-z^{L,\ell}+P_{k}^{L}(\varepsilon_{*}^{k}-P_{\ell}^{k}\varepsilon_{*}^{\ell}).$

Since $\Xi_{\ell}\subset\Xi_{k}$, $\exists!\,\varepsilon^{k}_{\ell}$ :
$z^{L,\ell}=z^{L,k}+P_{k}^{L}\varepsilon^{k}_{\ell}$666It is unique because of
eq. 8.. Hence we can write

$z^{L,k+1}-z^{L,\ell+1}=P_{k}^{L}(\varepsilon^{k}_{A}-\varepsilon^{k}_{B}),\quad\begin{cases}\varepsilon^{k}_{A}&=\varepsilon^{k}_{*}-\varepsilon^{k}_{\ell},\\\
\varepsilon^{k}_{B}&=P_{\ell}^{k}\varepsilon^{\ell}_{*}.\end{cases}$ (34)

Thus, by eq. 8,

$\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}\leq
d_{2}\|\varepsilon^{k}_{A}-\varepsilon^{k}_{B}\|_{\infty}.$

To estimate the right-hand-side of the above relation, we seek to apply lemma
2. For this, we consider $D$ to be a compact and convex set which contains in
its interior the sequence of iterates $\\{z^{L,k}\\}_{k=0}^{L+1}$, and define
the sets

$A_{\ell,k}:=\\{\varepsilon^{k}\in\mathbb{R}^{N_{k}}:z^{L,\ell}+P_{k}^{L}\varepsilon^{k}\in
D\\},\quad B_{\ell,k}:=\\{\varepsilon^{k}\in
A_{\ell,k}:\varepsilon^{k}=P_{\ell}^{k}\varepsilon^{\ell},\,\varepsilon^{\ell}\in\mathbb{R}^{N_{\ell}}\\}.$

Since $D$ is compact and convex, and the prediction operators are linear,
$A_{\ell,k}$ and $B_{\ell,k}$ are compact and convex sets in
$\mathbb{R}^{N_{k}}$. Notice that lemma 1 and eq. 8 imply that there are one-
to-one correspondences between $A_{\ell,k}\,\Leftrightarrow\,\Xi_{k}\cap D$
and $B_{\ell,k}\,\Leftrightarrow\,\Xi_{\ell}\cap D$. The sets $\Xi_{k}\cap D$
and $\Xi_{\ell}\cap D$ are compact and convex in $\mathbb{R}^{N_{L}}$.

Taking into account the relations above, we easily get

$z^{L,\ell+1}=z^{L,\ell}+P_{k}^{L}P_{\ell}^{k}\varepsilon^{\ell}_{*}\,\,\rightarrow\,\,\varepsilon^{k}_{B}\in
B_{\ell,k},\,\,z^{L,k+1}=z^{L,\ell}+P_{k}^{L}\varepsilon^{k}_{A}\,\,\rightarrow\,\,\varepsilon^{k}_{A}\in
A_{\ell,k}.$ (35)

Hence, considering the function

$F_{\ell,k}:\mathbb{R}^{N_{k}}\rightarrow\mathbb{R},\quad
F_{\ell,k}(\varepsilon^{k}):=F(z^{L,\ell}+P_{k}^{L}\varepsilon^{k})$

we have

$F_{\ell,k}(\varepsilon^{k}_{A})=F(z^{L,k+1})=\min_{z^{L}\in\Xi_{k}}F(z^{L}),\quad
F_{\ell,k}(\varepsilon^{k}_{B})=F(z^{L,\ell+1})=\min_{z^{L}\in\Xi_{\ell}}F(z^{L}),$

therefore

$F_{\ell,k}(\varepsilon^{k}_{A})=\min_{\varepsilon^{k}\in
A_{\ell,k}}F_{\ell,k}(\varepsilon^{k}),\quad
F_{\ell,k}(\varepsilon^{k}_{B})=\min_{\varepsilon^{k}\in
B_{\ell,k}}F_{\ell,k}(\varepsilon^{k}).$

Since
$F_{\ell,k}(\varepsilon^{k})=F_{k}(\varepsilon^{k}_{\ell}+\varepsilon^{k})$,
we have that $\nabla F_{\ell,k}(\varepsilon^{k}_{A})=\nabla
F_{k}(\varepsilon_{*}^{k})=0$. Hence, applying lemma 2, there exists
$c_{\ell,k}=c_{\ell,k}(A_{\ell,k},F_{\ell,k})$ such that

$\|\varepsilon_{A}^{k}-\varepsilon_{B}^{k}\|_{\infty}\leq c_{\ell,k}\
\text{dist}_{\infty}(\varepsilon_{A}^{k},B_{\ell,k}).$

Let $\varepsilon^{k}\in B_{l,k}$ and
$\varepsilon^{\ell}\in\mathbb{R}^{N_{\ell}}$:
$\varepsilon^{k}=P_{\ell}^{k}\varepsilon^{\ell}$. Then, using eq. 8 and eq. 35
we have

$||P_{\ell}^{k}\varepsilon^{l}-\varepsilon^{k}_{A}||_{\infty}\leq
d_{1}^{-1}||P_{k}^{L}(P_{\ell}^{k}\varepsilon^{l}-\varepsilon^{k}_{A})||_{\infty}=d^{-1}_{1}||z^{L,\ell}+P_{\ell}^{L}\varepsilon^{\ell}-z^{L,k+1}\|_{\infty},$

thus

$\displaystyle\text{dist}_{\infty}(\varepsilon_{A}^{k},B_{\ell,k})$
$\displaystyle\leq
d_{1}^{-1}\inf\left\\{\|z^{L,\ell}+P_{\ell}^{L}\varepsilon^{\ell}-z^{L,k+1}\|_{\infty}\
:\ P_{\ell}^{k}\varepsilon^{\ell}\in B_{k,\ell}\right\\}$
$\displaystyle=d_{1}^{-1}\text{dist}(z^{L,k+1},\Xi_{\ell}\cap D).$

Therefore

$\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}\leq
d_{2}||\varepsilon^{k}_{A}-\varepsilon^{k}_{B}||_{\infty}\leq
c_{\ell,k}\frac{d_{2}}{d_{1}}\text{dist}(z^{L,k+1},\Xi_{\ell}\cap D).$ (36)

To see that $c_{\ell,k}$ on the right hand side can be substituted by a
universal constant, independent of $\ell,k$, we proceed as follows: Recall
that, from lemma 2,
$c_{\ell,k}=c(A_{\ell,k},F_{\ell,k})=\sqrt{\rho_{\max}^{l,k}/\rho_{\min}^{l,k}}$
with

$\displaystyle\rho_{\min}^{l,k}=\min_{\xi^{k}\in
A_{\ell,k}}\min_{\varepsilon^{k}\neq
0}\frac{(\varepsilon^{k})^{T}\nabla^{2}F_{\ell,k}(\xi^{k})\varepsilon^{k}}{\|\varepsilon^{k}\|_{\infty}^{2}},\qquad\rho_{\max}^{l,k}=\max_{\xi^{k}\in
A_{\ell,k}}\max_{\varepsilon^{k}\neq
0}\frac{(\varepsilon^{k})^{T}\nabla^{2}F_{\ell,k}(\xi^{k})\varepsilon^{k}}{\|\varepsilon^{k}\|_{\infty}^{2}}.$

Taking into account that $\forall\,\xi^{k}\in\mathbb{R}^{N_{k}}$

$\quad\nabla^{2}F_{\ell,k}(\xi^{k})=(P_{k}^{L})^{T}\nabla^{2}F(z^{k})P_{k}^{L},\quad
z^{k}=z^{L,\ell}+P_{k}^{L}\xi^{k},$

and recalling that eq. 8 implies $P_{k}^{L}\varepsilon^{k}\neq 0$ when
$\varepsilon_{k}\neq 0$, we can write

$\displaystyle\frac{(\varepsilon^{k})^{T}\nabla^{2}F_{\ell,k}(\xi^{k})\varepsilon^{k}}{\|\varepsilon^{k}\|_{\infty}^{2}}$
$\displaystyle=\frac{(P_{k}^{L}\varepsilon^{k})^{T}\nabla^{2}F(z^{k})(P_{k}^{L}\varepsilon^{k})}{\|P_{k}^{L}\varepsilon^{k}\|_{\infty}^{2}}\frac{\|P_{k}^{L}\varepsilon^{k}\|_{\infty}^{2}}{\|\varepsilon^{k}\|_{\infty}^{2}}$
$\displaystyle\geq\frac{(P_{k}^{L}\varepsilon^{k})^{T}\nabla^{2}F(z^{k})(P_{k}^{L}\varepsilon^{k})}{\|P_{k}^{L}\varepsilon^{k}\|_{\infty}^{2}}d_{1}^{2}\geq\min_{\varepsilon^{L}\neq
0}\frac{(\varepsilon^{L})^{T}\nabla^{2}F(z^{k})\varepsilon^{L}}{\|\varepsilon^{L}\|_{\infty}^{2}}d_{1}^{2}.$

Since $\xi^{k}\in
A_{\ell,k}\,\Rightarrow\,z^{k}=z^{L,\ell}+P_{k}^{L}\xi^{k}\in\Xi_{k}\cap D$,
we have that

$\rho_{\min}^{l,k}\geq\rho_{\min}(F,D)d_{1}^{2},\quad\text{ with
}\quad\rho_{\min}(F,D)\\!\\!:=\min_{z\in D}\min_{\varepsilon^{L}\neq
0}\frac{(\varepsilon^{L})^{T}\nabla^{2}F(z)\varepsilon^{L}}{\|\varepsilon^{L}\|_{\infty}^{2}},$
(37)

proceeding analogously with the upper bounds, we have that

$\rho_{\max}^{l,k}\leq\rho_{\max}(F,D)d_{2}^{2},\quad\text{ with
}\quad\rho_{\max}(F,D)\\!\\!:=\max_{z\in D}\max_{\varepsilon^{L}\neq
0}\frac{(\varepsilon^{L})^{T}\nabla^{2}F(z)\varepsilon^{L}}{\|\varepsilon^{L}\|_{\infty}^{2}}.$
(38)

Thus, we obtain eq. 33 with

$C=C(F,D,{\cal
P}):=\sqrt{\frac{\rho_{\max}(F,D)}{\rho_{\min}(F,D)}}{\frac{d_{2}^{2}}{d_{1}^{2}}}.$

∎

This result allows us to relate the difference between sub-optimal solutions
to the size of certain detail coefficients.

######  Theorem 1.

Let $F\in{\cal C}^{2}(\mathbb{R}^{N_{L}},\mathbb{R})$ and ${\cal P}$ as in
proposition 3. Then, for any pair of indices $\ell,k$, $0\leq\ell\leq k<L$
there exists $C>0$, independent of $\ell,k$, such that

$\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}\leq
C\|(d^{\ell}(z^{L,k+1}-z^{L,0}),\ldots,d^{k-1}(z^{L,k+1}-z^{L,0}))\|_{\infty}.$
(39)

###### Proof.

From eq. 26 and eq. 27, we can write

$M_{\ell,L}z^{L,\ell+1}=(z^{\ell},\bar{d}^{\ell},\ldots,\bar{d}^{L-1}),\qquad
M_{k,L}z^{L,k+1}=(z^{k},\bar{d}^{k},\ldots,\bar{d}^{L-1}),$

with $z^{m}:=z^{m}_{*}+\varepsilon^{m}_{*}\in\mathbb{R}^{N_{m}}$, $m=k,\ell$
and $\bar{d}^{j}=d^{j}(z^{L,0})$. Since $\ell\leq k$, by eq. 5,

$\displaystyle M_{\ell,L}z^{L,k+1}$
$\displaystyle=(M_{\ell,k}z^{k},\bar{d}^{k},\ldots,\bar{d}^{L-1})=(x^{\ell},\hat{d}^{\ell},\ldots,\hat{d}^{k-1},\bar{d}^{k},\ldots,\bar{d}^{L-1}),\quad$

with $x^{\ell}\in\mathbb{R}^{N_{\ell}}$, $\hat{d}^{j}=d^{j}(z^{L,k+1})$ (see
remark Remark.). Then,

$M_{\ell}^{L}(z^{L,k+1}-z^{L,\ell+1})=(x^{\ell}-z^{\ell},\hat{d}^{\ell}-\bar{d}^{\ell},\cdots,\hat{d}^{k-1}-\bar{d}^{k-1},0,\cdots,0).$

Thus, $w^{\ell,k}:=z^{L,\ell+1}+P_{\ell}^{L}(x^{\ell}-z^{\ell})\in\Xi_{\ell}$
satisfies

$M_{\ell,L}(z^{L,k+1}-w^{\ell,k})=\left(0,\hat{d}^{\ell}-\bar{d}^{\ell},\ldots,\hat{d}^{k-1}-\bar{d}^{k-1},0,\ldots,0\right).$
(40)

Let $D$ be a compact set containing all the iterates
$\\{z^{L,m}\\}_{m=0}^{L+1}$ and all the points $w^{\ell,k}$,
$\forall\ell,k:0\leq\ell\leq k<L$. Applying proposition 3
$\exists\,C=C(F,D,{\cal P})$ such that

$\displaystyle\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}$ $\displaystyle\leq C\
\text{dist}(z^{L,k+1},\Xi_{\ell}\cap D)\leq
C\|z^{L,k+1}-w^{\ell,k}\|_{\infty}.$

Since
$\|z^{L,k+1}-w^{\ell,k}\|_{\infty}=\|M_{\ell,k}^{-1}M_{\ell,k}(z^{L,k+1}-w^{\ell,k})\|_{\infty}$,
by eq. 40 and eq. 7, we have

$\displaystyle\|z^{L,k+1}-z^{L,\ell+1}\|_{\infty}$ $\displaystyle\leq
C\left\|M_{\ell,L}^{-1}\left(0,\hat{d}^{\ell}-\bar{d}^{\ell},\ldots,\hat{d}^{k-1}-\bar{d}^{k-1},0,\ldots,0\right)\right\|_{\infty}$
$\displaystyle\leq
C\tilde{C}\left\|(d^{\ell}(z^{L,k+1}-z^{L,0}),\ldots,d^{k-1}(z^{L,k+1}-z^{L,0}))\right\|_{\infty},$

since $\hat{d}^{m}-\bar{d}^{m}=d^{m}(z^{L,k+1}-z^{L,0})$, by the linearity of
the MR transformations. ∎

Under smoothness assumptions, the decay of the detail coefficients is related
to the approximation properties of the prediction scheme, which allows to
obtain quantitative estimates for the distance between sub-optimal solutions.
.

######  Corollary 1.

Let $F\in{\cal C}^{2}(\mathbb{R}^{N_{L}},\mathbb{R})$ as in proposition 3, and
${\cal P}=\\{P_{k}^{k+1}\\}_{k\geq 0}$ the sequence of interpolatory
prediction operators given in section 2.1. If the initial guess
$\bar{z}=z^{L,0}$ and $z_{\min}=\operatorname*{arg\
min}\\{F(z),\,z\in\mathbb{R}^{N_{L}}\\}$ can be associated to the point
evaluations on ${\cal G}^{L}$ of sufficiently smooth functions, then for
$0\leq\ell<L$,

  1. 1.

$\|z_{\min}-z^{L,\ell+1}\|_{\infty}=O(h_{\ell+1}^{n+1})$,

  2. 2.

$\|z^{L,\ell+1}-z^{L,\ell}\|_{\infty}=O(h_{\ell}^{n+1})$,

  3. 3.

$\|\varepsilon^{\ell}_{*}\|_{\infty}=O(h_{\ell}^{n+1})$.

###### Proof.

If $z_{\min}-\bar{z}$ can be considered as point evaluations of a sufficiently
smooth function, we have that
$\|d^{j}(z_{\min}-\bar{z})\|_{\infty}=O(h_{j+1}^{n+1})$, $0\leq j<L$. To prove
item 1 we take $k=L$ in theorem 1. Since $z^{L,L+1}=z_{\min}$,
$\bar{z}=z^{L,0}$, we have

$\|z_{\min}-z^{L,\ell+1}\|_{\infty}\leq
C\|(d^{\ell}(z_{\min}-\bar{z}),\ldots,d^{L-1}(z_{\min}-\bar{z}))\|_{\infty}=O(h_{\ell+1}^{n+1}).$

item 2 and item 3 follow from the previous result, since

$\|z^{L,\ell}-z^{L,\ell+1}\|_{\infty}\leq\|z_{\min}-z^{L,\ell}\|_{\infty}+\|z_{\min}-z^{L,\ell+1}\|_{\infty}$

and $z^{L,\ell+1}-z^{L,\ell}=P_{\ell}^{L}\varepsilon^{\ell}_{*}$, hence, by
eq. 8,

$\|\varepsilon^{\ell}_{*}\|_{\infty}\leq
d_{1}^{-1}||P_{\ell}^{L}\varepsilon^{\ell}_{*}||_{\infty}=d_{1}^{-1}\
\|z^{L,\ell+1}-z^{L,\ell}\|_{\infty}.$

∎

We summarize the results shown in this section as follows: For (suitable)
convex optimization problems, the distance between sub-optimal solutions is
related to the size of certain detail coefficients. Provided that both the
solution and the initial guess are smooth, the $\infty$-distance between
consecutive sub-optimal solutions decreases when climbing up the MR ladder (at
a rate which depends on the properties of the prediction operators). In this
case, even though the auxiliary minimization problems eq. 24 involve an
increasing number of degrees of freedom when the resolution level increases,
it seems reasonable to assume that they will be efficiently solved by the
given optimization technique, due to the fact that the initial guess and the
solution of each auxiliary problem are increasingly closer.

Our theoretical results are numerically validated in the next section.

##  4 Numerical Experiments

The common setting for our numerical experiments is as follows: We seek to
solve eq. 1 by embedding the problem into an appropriate interpolatory MR
Framework and using the MR/OPT strategy with a given minimization tool (${\cal
D}$) and an initial guess, $\bar{z}$, which we assume supplied by the end-
user. We recall the following two important points in our approach:
1\. The solution of the auxiliary minimization problems eq. 24 can be carried
out using only the objective function $F$ (as a black box routine). Moreover,
in practice only the prediction rules of the MR transformation are really
required for the application of our algorithm.
2\.  The chosen optimizer, ${\cal D}$, is also a black-box routine in the
algorithm.

In this paper we shall consider two different optimization tools from MATLAB:
fminunc, a _steepest descent method_ which requires a gradient estimation, and
the general purpose optimization tool patternserach, which is a _coordinate
search method_. When using them as black-box tool, a series of parameters need
to be specified. We set the running parameters seeking to stop ${\cal D}$ when
the max-norm of the difference between two consecutive iterates within the
call falls below a specified tolerance ($tol_{\mathcal{D}}$)777In our
implementation the running parameters (consult its documentation for further
information) are TolX $=tol_{\mathcal{D}}$, TolFun=0 and
MaxFunEvals=MaxIter=Inf. For fminunc we additionally established
OptimalityTolerance=0, Algorithm=quasi-newton, FinDiffRelStep=1e-13 and
FiniteDifferenceType=central..

In addition, it seems reasonable to consider

$\|z^{L,k+1}-z^{L,k}\|\leq tol_{M}$ (41)

as a suitable stopping criteria for the MR/OPT strategy, in order to avoid
unnecessary calls to the optimizer. In our numerical experiments we always
choose $tol_{M}=tol_{\mathcal{D}}$.

The tests considered in the following subsections are taken from the
specialized literature. In all the cases considered, the MR/OPT strategy leads
to a more efficient computation than the direct use of ${\cal D}$ to solve eq.
1. To measure the performance of the MR/OPT strategy, we have considered, as
in [23], the total number of evaluations of the objective function $F$, since
often the remaining operations may have a negligible cost compared to that. In
addition, we also examine the ratio between the number of functional
evaluations required to solve each $k$-th auxiliary problem versus its number
of degrees of freedom, because it also provides very valuable information on
the behavior of the strategy.

All the computations have been done with MATLAB R2020b on a MacBook Air with a
1,8 GHz Intel Core i5 double core processor and a 8 GB 1600 MHz DDR3 memory.

###  4.1  Quadratic Optimization Problems

Boundary Value Problems (BVP) in ordinary and partial differential equations
are a common source of large scale optimization problems. In many cases, the
discrete setting corresponds to a quadratic minimization problem of the type
considered in proposition 1. Here, we shall consider 1D and 2D BVP with smooth
solutions in order to check numerically the theoretical results stated in the
previous section.

We consider first the 1D BVP (from [30]),

$\begin{cases}-u^{\prime\prime}(t)+2u(t)=f(t),&t\in(0,1)\\\
u(0)=u(1)=0.\end{cases}$ (42)

where $f(t):=10^{6}t(1-t)(t-1/2)(t-1/4)(3/4-t)$. Using the standard centered
second order discretization for $u^{\prime\prime}$ on the uniform grid
$\mathcal{G}:=(i/J)_{i=0}^{J}$ leads to the linear system

$(-z_{i-1}+2z_{i}-z_{i+1})J^{2}+2z_{i}=f(i/J),\qquad i=1,2,\ldots,J-1,$

for the unknown values $z_{i}=u(i/J)+O(h^{2})$, $h=1/J$, $1\leq i\leq J-1$
($z_{0}=z_{J}=0$ because of the boundary conditions). In practice, $J$ is
sufficiently large so that the (discrete) solution of the linear system is
sufficiently close to the true solution of the ODE, according to the end-user.

The coefficient matrix for the system $A\in\mathbb{R}^{(J-1)\times(J-1)}$, is
banded, symmetric and definite positive, hence the solution can be found by
solving a quadratic minimization problem of the type stated in proposition 1.
We compute the solution of this minimization problem using the MR/OPT strategy
in combination with the 1D Interpolatory MR setting described in section 2.1.
For this, we consider $J=J_{L}$, ${\cal G}={\cal G}_{L}$ as the finest mesh in
a ladder of nested uniform grids on $[0,1]$ of the form ${\cal
G}_{k}=(i/J_{k})_{i=0}^{J_{k}}$, $J_{k}=J_{k+1}/2$, $0\leq k<L$. Hence, the
discrete spaces in the corresponding 2D-interpolatory MR setting have
dimension $(J_{k}+1)$. In our computations, we take $J_{L}=128=2^{7}$, $L=5$
(i.e. $J_{0}=4$), $tol_{\mathcal{D}}=10^{-6}=tol_{M}$.

Considering the zero vector as initial data, we comply with the boundary
conditions, i.e $\bar{z}=z^{L,0}=\vec{0}\in\mathbb{R}^{N_{L}}$. At each
resolution level we enforce $\varepsilon^{k}_{0}=0=\varepsilon^{k}_{J_{k}}$
when applying the prediction operators and, as a result, the sub-optimal
solutions associated to the auxiliary optimization problems naturally satisfy
the boundary conditions. In fig. 1, we display $z^{L,k+1}$ for $k=0,1,2$ and
the different prediction operators in section 2.1. For $n=1$, the sub-optimal
solution is a piecewise linear function, while for $n=3,5$ we can clearly
observe that the sub-optimal solutions are smoother888These facts are related
to the properties of the prediction scheme.. The plots shown correspond to
${\cal D}=$fminunc, but are indistinguishable from those obtained with ${\cal
D}=$pattersearch.

|  |
---|---|---
Figure 1: 1D BVP ($\mathcal{D}=$ fminunc). The sub-optimal solutions
$z^{L,1},z^{L,2},z^{L,3}$, from left to right, for $n=1,3,5$.

|
---|---
|
$n\backslash k$ | 1 | 2 | 3 | 4 | 5 | Theoretical rate
---|---|---|---|---|---|---
1 | -0.51 | 1.84 | 1.80 | 2.00 | 1.98 | 2
3 | 0.48 | 3.54 | 2.47 | 3.37 | 3.80 | 4
5 | 0.45 | 4.24 | 5.33 | 5.78 | 6.09 | 6

Figure 2: 1D BVP. $tol_{\cal D}=tol_{M}=10^{-6}$. Horizontal axis, $k$
(resolution level). Top row, $\mathcal{D}=\texttt{fminunc}$; Bottom row,
$\mathcal{D}=\texttt{patternsearch}$. Left column, $\infty$-distance between
$z^{L,k+1}$ and $z_{\min}$. Right column, $\infty$-distance between
consecutive sub-optimal solutions. Table, numerical estimation of $r$ in eq.
43, for $\mathcal{D}=\texttt{fminunc}$.

In fig. 2, we clearly observe that the decay rates of (the max-norm of) the
difference between the sub-optimal solutions and $z_{\min}$999Computed with
$A\backslash b$ of MATLAB., and the difference between two consecutive sub-
optimal solutions, depend on the prediction scheme, but not on the
optimization tool, just as predicted by our theoretical results. The table
below the plots, shows a numerical estimation of the decay rate corresponding
to the values $||z^{L,k+1}-z^{L,k}||_{\infty}$ obtained when using fminunc
(the results for patternsearch are similar and are not shown), computed as

${r=\log_{2}\frac{\|z^{L,k}-z^{L,k-1}\|_{\infty}}{\|z^{L,k+1}-z^{L,k}\|_{\infty}}.}$
(43)

Since the optimization problem fulfills all the hypothesis stated in Theorem 1
and Corollary 1, the numerical decay rates coincide with the expected decay
rates of the detail coefficients eq. 17 for each choice of the prediction
scheme.

|
---|---
# $F$-evaluations by | $n=1$ | $n=3$ | $n=5$ | Direct-$\mathcal{D}$
---|---|---|---|---
fminunc | 38 678 | 17 089 | 8 910 | 49 980
patternsearch | 101 307 742 | 7 938 578 | 1 063 433 | 176 168 800

Figure 3: 1D BVP, as in fig. 2. Ratio between the number of functional
evaluations and the number of degrees of freedom involved in the solution of
the $k$-th auxiliary problem, versus $k$. Left:
$\mathcal{D}=\texttt{fminunc}$; Right: $\mathcal{D}=\texttt{patternsearch}$.
Table: Total number of functional evaluations required to find $z^{L,L+1}$ in
each case.

In fig. 3 we examine the efficiency of the MR/OPT strategy. The table below
the plots clearly shows that the MR/OPT strategy leads to a significant
improvement in efficiency compared to the direct application of the
optimization tool. Since fminunc is a much more efficient optimizer for this
type of problems, its use leads to a smaller number of functional evaluations,
both when applied directly or when used in the MR/OPT strategy. It is worth
observing that the MR/OPT strategy leads to large improvements in efficiency:
the total number of evaluations carried out with the MR/OPT for $n=5$, is less
than 18% for fminunc and 0.7% for patternsearch.

We can see in the plots of fig. 3 that the number of functional evaluations
over the number of degrees of freedom increases in a similar way for $n=1,3,5$
for the lower resolution levels, but this ratio decays at the higher
resolution levels for $n=3,5$. A plausible explanation is that the distance
between the initial guess and the solution of the $k$-th auxiliary problem
decreases with $n$ and $k$, as specified in corollary 1, hence the
corresponding optimization problems can be solved in less iterations. Thus,
the efficiency improves with the approximation order of the prediction scheme
used.

Notice, however, that even in the ‘least efficient’ case ($n=1$), the MR/OPT
reaches the desired solution more efficiently than with the direct optimizer,
without any additional effort by the end-user.

Our second test case is the P2D in [22],

$\begin{array}[]{rcl}-(u_{xx}(x,y)+u_{yy}(x,y))&=&f(x,y),\\\\[2.0pt]
u(x,y)&=&0,\end{array}\qquad\begin{array}[]{l}(x,y)\in\text{int}(\Omega),\\\\[2.0pt]
(x,y)\in\ \partial\Omega,\end{array}$ (44)

$\Omega=[0,1]\times[0,1]$, $\partial\Omega$ and $\text{int}(\Omega)$ denoting
the boundary and the interior of $\Omega$, but with a smooth, non-polynomial
right hand side, $f(x,y)=\sin(4\pi x(1-x)y(1-y))$.

A standard discretization of eq. 44 (using the classical 5-point Laplacian,
see [22]), on a uniform grid $\mathcal{G}:=(i/J,j/J)_{i,j=0}^{J}$, leads to a
system of equations that can be written in matrix form as $Az=b$,
$z(i,j)=u(i/J,j/J)+O(h^{2})$, where $h=1/J$ is the uniform mesh spacing of
${\cal G}$ and $A$ is symmetric and definite positive. As in the 1D case,
finding $z$ may be formulated as a quadratic, convex, minimization problem,
but in the 2D case there are $N=(J-1)^{2}$ degrees of freedom, so that solving
the associated minimization problem is a much more demanding task, for the
same accuracy in the resulting approximation.

We solve the minimization problem using the MR/OPT strategy in combination
with the 2D-tensor product interpolatory MR described in section 2.1. For
this, we consider $J=J_{L}$, ${\cal G}={\cal G}_{L}$ as the finest mesh in a
ladder of nested uniform grids on $\Omega$ of the form ${\cal
G}_{k}=(i/J_{k},j/J_{k})_{i,j=0}^{J_{k}}$, $J_{k}=J_{k+1}/2$, $0\leq k<L$. The
discrete spaces in the corresponding 2D-interpolatory MR setting have
dimension $(J_{k}+1)^{2}$. In our computations, we take zero initial data,
$J_{L}=128$, $L=5$, $tol_{M}=tol_{\mathcal{D}}=10^{-7}$ and ${\cal
D}=$fminunc.

fig. 4 shows several plots analogous to those in fig. 2-fig. 3, where we
observe the same behavior as in the 1D case: The actual errors between the
auxiliary solutions and the true solution decrease, and also the efficiency of
the MR/OPT strategy increases, with the approximation order of the prediction
scheme. The decay rate of the max-norm errors between the sub-optimal
solutions and the true solution, and also between successive sub-optimal
solutions, coincides with the order of approximation of the prediction
operators (see table 1), which is consistent with corollary 1.

It is worth noticing that in this test case the max-norm of the difference
between sub-optimal solutions falls below $tol_{M}=tol_{\cal D}$ at $k=3$ for
$n=5$, and $k=4$ for $n=3$, which makes the MR/OPT strategy stop101010For
larger values of $tol_{\cal D}=tol_{M}$, the strategy stop even at lower
resolution levels.. We also remark the very large increment in efficiency
obtained for $n=3,5$, compared to $n=1$, which is probably due to the
smoothness of the solution. For $n=5$, the solution has been obtained with
only $0.1\%$ of the functional evaluations required by the direct use of
fminunc.

|  |
---|---|---

n | 1 | 3 | 5 | Direct-$\mathcal{D}$
---|---|---|---|---
# of $F$-evaluations | 2 742 108 | 41 206 | 11 136 | 12 581 010

Figure 4:  2D Poisson problem eq. 44, $L=5$, $tol_{M}=tol_{\mathcal{D}}=10^{-7}$. ${\cal D}$ is fminunc.  $n\backslash k$ | 1 | 2 | 3 | 4 | 5 | Theoretical rate
---|---|---|---|---|---|---
1 | 1.84 | 1.91 | 1.89 | 1.95 | 1.98 | 2
3 | 5.27 | 3.36 | 4.14 | 4.27 | - | 4
5 | 6.72 | 7.59 | 3.62 | - | - | 6
Table 1: Test-case: 2D Poisson problem. Decay rate eq. 43 for the data in fig.
4.

###  4.2 Non-quadratic problems

In this subsection we shall consider: a convex, but non-quadratic,
minimization problem and a non-convex problem, in order to evaluate the
performance of our MR/OPT strategy, in more general situations

We consider first the MINS problem in [22], which arises from a discretization
of the following minimal surface problem

$\min_{u}\int_{\Omega}\sqrt{1+\|\nabla u(x,y)\|_{2}^{2}}\ d(x,y),$ (45)

with the boundary conditions

$u_{0}(x,y)=\begin{cases}x(1-x),&\text{if }y\in\\{0,1\\},\\\
0,&\text{otherwise}.\end{cases}$

Its solution is approximated in [22] by the solution of the minimization
problem that results from considering as objective function

$F(z):=\frac{1}{2J^{2}}\sum_{i,j=0}^{J-1}\sqrt{1+a^{2}+b^{2}}+\sqrt{1+c^{2}+d^{2}},$

with

$\displaystyle a=$ $\displaystyle J(z_{i,j+1}-z_{i,j}),\qquad$ $\displaystyle
b=J(z_{i+1,j+1}-z_{i,j+1}),$ $\displaystyle c=$ $\displaystyle
J(z_{i+1,j+1}-z_{i+1,j}),\quad$ $\displaystyle d=J(z_{i+1,j}-z_{i,j}).$

We use the interpolatory 2D framework and the MR/OPT strategy to solve this
nonlinear, convex, optimization problem, with $J_{L}=128$, $L=5$,
$tol_{M}=tol_{\mathcal{D}}=10^{-6}$ and
$z^{L,0}_{i,j}=\frac{i}{J}\left(1-\frac{i}{J}\right)$,
$i,j\in\\{0,1,\ldots,J\\}$, which complies with the specified boundary
conditions. The solution obtained with MR/OPT and $\cal D$=fminunc with the
above parameters is represented in fig. 7-left.

In this case, we do not know the exact solution and in fig. 5-left we only
plot the max-norm of the difference between sub-optimal solutions. The
observed decay rate is $r\approx 2$, for all prediction schemes. To explain
this apparent contradiction with the results in corollary 1, we carry out a
numerical inspection of the derivatives of the solution shown in fig. 7-left.
We observe large spikes in third order derivatives at the four corners of the
unit square, which is an indication that the obtained discrete data can only
be associated to a $\mathcal{C}^{2}$ function. Hence, the observed decay rate
is, in fact, the one predicted by the theoretical result in corollary 1 since,
in this case, the interpolation error corresponding to an $n$-th degree
polynomial $(n\geq 1)$ is only $O(h^{2})$, with $h$ the distance between the
points in the grid.

The plot in fig. 5-right shows, again, the increase in efficiency of the
MR/OPT strategy for $n=3,5$, which we assume to be related to the fact that
the errors between consecutive sub-optimal solutions are smaller than for
$n=1$, hence the corresponding auxiliary optimization problems require less
iterations. For $n=5$, the solution of the original optimization problem,
shown in fig. 7-right, is obtained with less than 2 % of the total number of
evaluations of $F$, with respect to the direct aplication of ${\cal D}$.

|
---|---
interpolation degree | n=1 | n=3 | n=5 | Direct-$\mathcal{D}$
---|---|---|---|---
total # evaluations | 2 417 132 | 235 771 | 180 990 | 10 226 103

Figure 5: Minimal surface problem. $L=5$, $tol_{M}=tol_{\mathcal{D}}=10^{-6}$
and ${\cal D}$ is fminunc.

Our last test concerns a nonlinear non-convex optimization problem from [22,
section 5.1.4]: the MOREBV problem, which is obtained from a discretization on
a uniform grid on $\Omega=[0,1]\times[0,1]$ of the following BVP

$\begin{array}[]{rcl}-(u_{xx}(x,y)+u_{yy}(x,y))+\frac{1}{2}(u(x,y)+x+y+1)^{3}&=&0,\quad(x,y)\in\text{int}(\Omega),\\\\[2.0pt]
u(x,y)&=&0,\quad(x,y)\in{\partial\Omega}.\end{array}$ (46)

Using the classical 5-point discretization of the Laplacian, the resulting
system of nonlinear equations can be rewritten as a nonlinear least-squares
problem with a non-convex objective function given by the expression (see [22]
and references therein). The theoretical results of section 3.1 do not apply
in this case.

$F(z):=\sum_{i,j=1}^{J-1}\Big{(}(4z_{i,j}-z_{i-1,j}-z_{i+1,j}-z_{i,j-1}-z_{i,j+1})\\\
+\frac{1}{2J^{2}}\left(z_{i,j}+i/J+j/J+1\right)^{3}\Big{)}^{2}.$ (47)

We solve this problem by using the MR/OPT strategy within the 2D interpolatory
MR framework used in the previous two test cases. As before, setting
$\bar{z}=\vec{0}$ leads to a sequence of sub-optimal solutions satisfying the
boundary conditions. The solution obtained with the same parameters as in the
previous test case ($n=5$) is shown in fig. 7-right.

In fig. 6 we show the results obtained for this test case. The left plot in
this figure shows that when taking $n=1$ the distance between consecutive sub-
optimal solutions does not decay111111The set of sub-optimal solutions can be
found in the Appendix., but it does decrease for $n=3,5$ (at a rate $r\approx
2$). As a consequence, the number of calls to the objective function increases
dramatically for $n=1$, while it remains moderate for $n=3,5$. More research
is required to explain the behavior of the MR/OPT strategy for this non-convex
optimization problem, for which no theoretical results are provided in this
paper. We have observed that, as in the MINS problem, certain third order
finite differences of $z^{L,L+1}$ also exhibit peaks in the corners of
$\Omega$. We do remark that the solution corresponding at $n=5$ has been
obtained with less than $0.3\%$ of the total number of functional evaluations
than for $n=1$ (which is, at least, less than 69% of the total number of
evaluations for the direct use of ${\cal D}$).

|
---|---
interpolation degree | n=1 | n=3 | n=5 | Direct-$\mathcal{D}$
---|---|---|---|---
total # evaluations | 201 583 677 | 1 168 621 | 495 258 | 294 879 519

Figure 6:  MOREBV problem. $L=7$, $tol_{M}=tol_{\mathcal{D}}=10^{-6}$ and
${\cal D}$ is fminunc.

Figure 7: The computed solution, $z^{L,L+1}$, of the minimal surface (left)
and the MOREBV (right) problems taking $n=5$.

##  5 Conclusions And Perspectives

We have described and analyzed a multilevel strategy (MR/OPT) which can be
used to reduce the cost of solving large scale optimization problems. MR/OPT
uses Harten’s MRF to define a hierarchy of auxiliary problems that provide a
sequence of sub-optimal solutions, one at each resolution level, so that at
the last step the full (large-scale) optimization problem is solved, but with
an initial guess that is much closer to the desired solution than the original
initial guess (which is assumed to be provided by the end-user).

The solutions of the auxiliary problems depend on the prediction operators of
the chosen MRF. The results in this paper clearly demonstrate that the
accuracy of the prediction scheme plays a significant role in the overall
performance of the MR/OPT strategy, when the initial guess and the solution
can be interpreted as point-value discretizations of sufficiently smooth
functions.

Harten’s MRF provides a flexible and versatile multilevel structure that
allows to consider both the optimization tool and the objective function eq. 1
as black boxes, which is an important asset in engineering applications. If an
optimization problem such a eq. 1 can be solved with a specific optimization
tool, ${\cal D}$, MR/OPT will reduce the cost of finding its solution. For
convex optimization problems, we have presented several theoretical results on
the properties and the performance of the MR/OPT strategy, which have been
numerically validated in a series of examples from the specialized literature.
More analytic work is necessary in order to determine the properties of the
MR/OPT strategy in more general situations, as well as the influence of the
underlying minimization code or the specification of constraints.

The strategy can be improved in various ways, that we intend to investigate in
future papers. In particular, it might be advantageous to consider a level-
dependent tolerance criterion for the optimizer $\mathcal{D}$: large
tolerances in coarse levels and small ones in fine levels. Following this
idea, to consider different optimizers might also be profitable: e.g. global
searchers in the first levels and, later, fast/local optimizers. Another line
of research concerns the use of nonlinear reconstruction techniques instead of
the linear ones used in this paper, since they might be more efficient in
problems where the solution presents discontinuities or it is important to
preserve specific features of the solution, such as monotonicity or convexity.

## References

  * [1] Remi Abgrall, Pietro Marco Congedo, and Gianluca Geraci.  A one-time truncate and encode multiresolution stochastic framework.  Journal of Computational Physics, 257:19–56, 2014.
  * [2] Francesc Aràndiga and Rosa Donat.  Nonlinear multiscale decompositions: the approach of A. Harten.  Numerical Algorithms, 23(2-3):175–216, 2000.
  * [3] Francesc Aràndiga, Rosa Donat, and Ami Harten.  Multiresolution based on weighted averages of the hat function i: Linear reconstruction techniques.  SIAM Journal on Numerical Analysis, 36(1):160–203, 1998.
  * [4] Francesc Aràndiga, Rosa Donat, and Ami Harten.  Multiresolution based on weighted averages of the hat function i: Linear reconstruction techniques.  SIAM Journal on Numerical Analysis, 36(1):160–203, 1998.
  * [5] Francesc Aràndiga, Pep Mulet, and Vicent Renau.  Lossless and near-lossless image compression based on multiresolution analysis.  Journal of Computational and Applied Mathematics, 242:70–81, 2013\.
  * [6] Barna L. Bihari and Ami Harten.  Multiresolution schemes for the numerical solution of 2-d conservation laws i.  SIAM J. Sci. Comput., 18:315–354, 1997.
  * [7] Barna L Bihari and Donald Schwendeman.  Multiresolution schemes for the reactive euler equations.  Journal of Computational Physics, 154(1):197–230, 1999.
  * [8] Achi Brandt.  Multi-level adaptive solutions to boundary-value problems.  Mathematics of computation, 31(138):333–390, 1977.
  * [9] WL Briggs.  A Multigrid Tutorial.  1987\.
  * [10] Raymond H Chan and Ke Chen.  A multilevel algorithm for simultaneously denoising and deblurring images.  SIAM Journal on Scientific Computing, 32(2):1043–1063, 2010.
  * [11] Tony F Chan and Ke Chen.  An optimization-based multilevel algorithm for total variation image denoising.  Multiscale Modeling & Simulation, 5(2):615–645, 2006.
  * [12] Cheng Chen, Zaiwen Wen, and Ya-xiang Yuan.  A general two-level subspace method for nonlinear optimization.  Journal of Computational Mathematics, 36(6), 2018.
  * [13] Guillaume Chiavassa and Rosa Donat.  Point value multiscale algorithms for 2d compressible flows.  SIAM Journal on Scientific Computing, 23(3):805–823, 2001.
  * [14] Albert Cohen, Sidi Kaber, Siegfried Müller, and Marie Postel.  Fully adaptive multiresolution finite volume schemes for conservation laws.  Mathematics of Computation, 72(241):183–225, 2003.
  * [15] B Colson, M Porcelli, and Ph L Toint.  Aircraft fuselage sizing with multilevel optimization.  2013\.
  * [16] Johan M De Villiers, Karin M Goosen, and Ben M Herbst.  Dubuc–deslauriers subdivision for finite sequences and interpolation wavelets on an interval.  SIAM Journal on Mathematical analysis, 35(2):423–452, 2003.
  * [17] Johan M De Villiers, Karin M Goosen, and Ben M Herbst.  Dubuc–deslauriers subdivision for finite sequences and interpolation wavelets on an interval.  SIAM Journal on Mathematical analysis, 35(2):423–452, 2003.
  * [18] Rosa Donat and Sergio López-Ureña.  High-accuracy approximation of piecewise smooth functions using the truncation and encode approach.  Applied Mathematics and Nonlinear Sciences, 2(2):367–384, 2017\.
  * [19] Rosa Donat, Sergio López-Ureña, and Marc Menec.  A novel multi-scale strategy for multi-parametric optimization.  In European Consortium for Mathematics in Industry, pages 593–600. Springer, 2016.
  * [20] Nira Dyn.  Subdivision schemes in cagd.  Advances in numerical analysis, 2:36–104, 1992.
  * [21] Nira Dyn.  Subdivision schemes in cagd.  Advances in numerical analysis, 2:36–104, 1992.
  * [22] Emanuele Frandi and Alessandra Papini.  Coordinate search algorithms in multilevel optimization.  Optimization Methods and Software, 29(5):1020–1041, 2014.
  * [23] Emanuele Frandi and Alessandra Papini.  Improving direct search algorithms by multilevel optimization techniques.  Optimization Methods and Software, 30(5):1077–1094, 2015.
  * [24] Serge Gratton, Annick Sartenaer, and Philippe L Toint.  Recursive trust-region methods for multiscale nonlinear optimization.  SIAM Journal on Optimization, 19(1):414–444, 2008.
  * [25] Wolfgang Hackbusch.  Convergence of multigrid iterations applied to difference equations.  Mathematics of Computation, 34(150):425–440, 1980.
  * [26] Xiaocong Han and David W Zingg.  An adaptive geometry parametrization for aerodynamic shape optimization.  Optimization and Engineering, 15(1):69–91, 2014.
  * [27] Ami Harten.  Multiresolution representation of data: A general framework.  SIAM Journal on Numerical Analysis, 33(3):1205–1256, 1996.
  * [28] S López-Ureña, JR Torres-Lapasió, R Donat, and MC García-Alvarez-Coque.  Gradient design for liquid chromatography using multi-scale optimization.  Journal of Chromatography A, 1534:32–42, 2018.
  * [29] DA Masters, NJ Taylor, TCS Rendall, and Christian B Allen.  Multilevel subdivision parameterization scheme for aerodynamic shape optimization.  AIAA Journal, 55(10):3288–3303, 2017.
  * [30] Stephen G Nash.  A multigrid approach to discretized optimization problems.  Optimization Methods and Software, 14(1-2):99–116, 2000.

##  Appendix A On the Property S for subdivision schemes

For completeness, in this section we provide a proof, without entering into
too many details, which is well-known in subdivision schemes theory.

######  Lemma 3.

Let ${\cal P}=(P_{k}^{k+1})_{k=0}^{\infty}$ be the sequence of linear,
interpolatory, prediction operators defined in Section 2.1. Then there exists
$C>0$ such that

$\|\varepsilon^{\ell}\|_{\infty}\leq\|P_{\ell}^{k}\varepsilon^{\ell}\|_{\infty}\leq
C\|\varepsilon^{\ell}\|_{\infty},\qquad\forall 0\leq\ell\leq k.$

###### Proof.

Any interpolatory prediction operator fulfills

$\|\varepsilon^{\ell}\|_{\infty}=\sup_{0\leq i\leq
J_{\ell}}|\varepsilon_{i}^{\ell}|=\sup_{0\leq i\leq
J_{\ell}}|(P_{\ell}^{\ell+1}\varepsilon^{\ell})_{2i}|\leq\sup_{0\leq i\leq
J_{\ell+1}}|(P_{\ell}^{\ell+1}\varepsilon^{\ell})_{i}|=\|P_{\ell}^{\ell+1}\varepsilon^{\ell}\|_{\infty}.$

Applying this inductively, we obtain the first inequality

$\|\varepsilon^{\ell}\|_{\infty}\leq\|P_{\ell}^{\ell+1}\varepsilon^{\ell}\|_{\infty}\leq\|P_{k-1}^{k}\cdots
P_{\ell}^{\ell+1}\varepsilon^{\ell}\|_{\infty}=\|P_{\ell}^{k}\varepsilon^{\ell}\|_{\infty}.$

The second inequality is a consequence of the fact that the subdivision
process associated to the recursive application of the predictions operators
is convergent (see e.g. [21, 4, 17]), which means that for each vector
$\varepsilon^{\ell}$ there exists a continuous function, which we denote by
$P_{\ell}^{\infty}\varepsilon^{\ell}\in\mathcal{C}([0,1],\mathbb{R})$, such
that

$\lim_{\ell<m\to\infty}\sup_{0\leq i\leq
J_{m}}|(P_{\ell}^{m}\varepsilon^{\ell})_{i}-(P_{\ell}^{\infty}\varepsilon^{\ell})(h_{m}i)|=0,\quad$

with $h_{m}=J_{m}^{-1}$ being the grid spacing at $m$-level. For interpolatory
prediction operators
$(P_{\ell}^{k}\varepsilon^{\ell})_{i}=(P_{\ell}^{\infty}\varepsilon^{\ell})(h_{k}i)$,
$0\leq i\leq J_{k}$, $0\leq k$, by the interpolation condition, hence

$\|P_{\ell}^{k}\varepsilon^{\ell}\|_{\infty}=\sup_{0\leq i\leq
J_{k}}|(P_{\ell}^{k}\varepsilon^{\ell})_{i}|=\sup_{0\leq i\leq
J_{k}}|(P_{\ell}^{\infty}\varepsilon^{\ell})(h_{k}i)|\leq\sup_{t\in[0,1]}|(P_{\ell}^{\infty}\varepsilon^{\ell})(t)|.$

Since
$\varepsilon^{\ell}=\sum_{i=0}^{J_{\ell}}\varepsilon_{i}\delta_{i}^{\ell}$,
with $\delta_{i}^{\ell}$ the canonical $i$-th vector in
$\mathbb{R}^{N_{\ell}}$, and the prediction operators are linear, the result
will follow from the properties of the limit functions
$\phi_{i}^{\ell}=P_{\ell}^{\infty}\delta_{i}^{\ell}$. It can be proven (see
[4, 17]) that these limit functions have compact support (narrower than
$2nh_{\ell}$, with $n$ as in Section 2.1) and satisfy a refinability property
which ensures that there exists $\tilde{C}>0$ such that
$||\phi^{\ell}_{i}||_{\infty}\leq\tilde{C}$, $\forall i,0\leq i\leq J_{\ell}$,
$\forall\ell\geq 0$. Given $\varepsilon^{\ell}\in\mathbb{R}^{N_{\ell}}$,

$P_{\ell}^{\infty}\varepsilon^{\ell}(t)=\sum_{i=0}^{J_{\ell}}\varepsilon^{\ell}_{i}\phi_{i}^{\ell}(t),\quad
t\in[0,1]\quad\rightarrow\quad|(P_{\ell}^{\infty}\varepsilon^{\ell})(t)|\leq\|\varepsilon^{\ell}\|_{\infty}\sum_{i=0}^{J_{\ell}}|\phi_{i}^{\ell}(t)|$

the result follows from observing that for any $t\in[0,1]$,
$\phi^{\ell}_{i}(t)\neq 0$ only for at most $2n$ such functions (see [4, 17]).
Hence for some $C>0$, which does not depend on $\ell$, the second inequality
is obtained. ∎

##  Appendix B Limit function basis of the prediction operators

The prediction operators described in section 2.1 are based on $n$-th degree
polynomial interpolation, with $n=1,3,5$. Since they are linear operators, the
refined data can be written in terms of a basis. In figs. 8 and 9, some limit
function basis for $n=3,5$, respectively, are shown.

|
---|---
|
|
Figure 8: Some of the basis functions for the space of generated functions by the prediction schemes in section 2.1 for $n=3$. The starting level has 17 equal-spaced nodes in [0,1], and the prediction operators are applied 10 times, so the final level has grid points. |
---|---
|
|
Figure 9: Some of the basis functions for the space of generated functions by
the prediction schemes in section 2.1 for $n=5$. The starting level has 17
equal-spaced nodes in [0,1], and the prediction operators are applied 10
times, so the final level has 65537 grid points.

##  Appendix C Smoothness of the MINS and MOREBV solutions

In this section we discuss about the regularity of the solutions to the
problems ‘MINS’ eq. 45 and ‘MOREBV’ eq. 46.

According the boundary conditions, the derivatives
$\frac{\partial^{3}}{\partial x^{3}}$ $\frac{\partial^{3}}{\partial y^{3}}$ of
the solution $u$ to the continuous problem must be zero in the four corners
$(0,0),(1,0),(0,1),(1,1)$. In the pictures in fig. 10 there is a sudden change
in the value of the derivatives in the corners. On the one direction ($x$ or
$y$, depending on the picture), the derivative is zero, while in the other
direction, the absolute value grows when approaching to the corner. If there
exists some continuous function which point evaluations are the discrete
solution $z_{\min}$, then its third derivatives seem to be not continuous in
the four corners. Hence, the solutions are not $\mathcal{C}^{3}$ globally. For
comparison, we show the graphics of $\frac{\partial^{2}}{\partial x^{2}}$
$\frac{\partial^{2}}{\partial y^{2}}$ in Figure fig. 11, which seem to be
continuous, but not differentiable.

|
---|---
|
Figure 10: The derivatives $\frac{\partial^{3}}{\partial x^{3}}$ $\frac{\partial^{3}}{\partial y^{3}}$ (left and right columns, respectively) in the interval $[0,1]\times[0,1]$ of the solutions to the problems ‘MINS’ eq. 45 and ‘MOREBV’ eq. 46 (top and bottom rows). The solutions were computed on a $65\times 65$ uniform grid up to machine precision, solving the minimization problem with the direct use of fminunc.  |
---|---
|
Figure 11: The derivatives $\frac{\partial^{2}}{\partial x^{2}}$
$\frac{\partial^{2}}{\partial y^{2}}$ (left and right columns, respectively)
in the interval $[0,1]\times[0,1]$ of the solutions to the problems ‘MINS’ eq.
45 and ‘MOREBV’ eq. 46 (top and bottom rows). The solutions were computed on a
$65\times 65$ uniform grid up to machine precision, solving the minimization
problem with the direct use of fminunc.

##  Appendix D Sub-optimal solutions of the MOREBV problem

In this section we show the sub-optimal solutions found by MR/OPT when solving
the MOREBV problem as in section 4.2.

Observe in fig. 12 that, for the linear prediction operators in section 2.1
with $n=1$, how slowly the sub-optimal solutions approaches the solution (pay
attention to the z-axis). The low regularity of the data generated by this
prediction operators may explain this behavior, since the continuous
optimization problem requires the computation of a Laplacian.

However, we can see in figs. 13 and 14 that for $n=3,5$ the sub-optimal
solution converges quickly to $z_{\min}$.

$z^{L,1}$ | $z^{L,2}$
---|---
|
$z^{L,3}$ | $z^{L,4}$
|
$z^{L,5}$ | $z^{L,6}$
|
Figure 12: The suboptimal solutions to the problem ‘MOREBV’ eq. 46 for $n=1$.  $z^{L,1}$ | $z^{L,2}$
---|---
|
$z^{L,3}$ | $z^{L,4}$
|
Figure 13: Some suboptimal solutions to the problem ‘MOREBV’ eq. 46 for $n=3$.  $z^{L,1}$ | $z^{L,2}$
---|---
|
$z^{L,3}$ | $z^{L,4}$
|
Figure 14: Some suboptimal solutions to the problem ‘MOREBV’ eq. 46 for $n=5$.