On high-order multilevel optimization strategies

TL;DR

This paper introduces a multilevel framework for high-order unconstrained optimization, reducing computational costs by leveraging problem hierarchies, with theoretical convergence guarantees and practical efficiency demonstrated through numerical experiments.

Contribution

It extends high-order optimization methods with a multilevel approach, providing convergence analysis and demonstrating computational savings over traditional methods.

Findings

The multilevel methods achieve convergence to stationary points.

Numerical results show significant computational savings.

Theoretical analysis confirms convergence and complexity bounds.

Abstract

We propose a new family of multilevel methods for unconstrained minimization. The resulting strategies are multilevel extensions of high-order optimization methods based on q-order Taylor models (with q >= 1) that have been recently proposed in the literature. The use of high-order models, while decreasing the worst-case complexity bound, makes these methods computationally more expensive. Hence, to counteract this effect, we propose a multilevel strategy that exploits a hierarchy of problems of decreasing dimension, still approximating the original one, to reduce the global cost of the step computation. A theoretical analysis of the family of methods is proposed. Specifically, local and global convergence results are proved and a complexity bound to reach first order stationary points is also derived. A multilevel version of the well known adaptive method based on cubic regularization (ARC, corresponding to q = 2 in our setting) has been implemented. Numerical experiments clearly highlight the relevance of the new multilevel approach leading to considerable computational savings in terms of floating point operations compared to the classical one-level strategy.

Tables

Table 1. Table 1: Solution of the minimization problem ( 33 ) with the one level ARC method and a four level ARC (MARC) (case of n h = 4096 subscript 𝑛 ℎ 4096 n_{h}=4096 ) with u ¯ 1 = 1 ​ rand ​ ( n h , 1 ) subscript ¯ 𝑢 1 1 rand subscript 𝑛 ℎ 1 \bar{u}_{1}=1\,\mathrm{rand}(n_{h},1) , u ¯ 2 = 3 ​ rand ​ ( n h , 1 ) subscript ¯ 𝑢 2 3 rand subscript 𝑛 ℎ 1 \bar{u}_{2}=3\,\mathrm{rand}(n_{h},1) . i ​ t T 𝑖 subscript 𝑡 𝑇 it_{T} denotes the average number of iterations over ten simulations, i ​ t f 𝑖 subscript 𝑡 𝑓 it_{f} the average number of iterations in which the fine level model has been used, 𝚁𝙼𝚂𝙴 𝚁𝙼𝚂𝙴 {\tt RMSE} the root-mean square error with respect to the true solution and save the ratio between the total number of floating point operations required for the Cholesky factorizations in ARC and MARC, respectively.

			$n_h=4096$
$u_0$	Method	$i{t}_T/i{t}_f$	RMSE	save
${\overline{u}}_1$	ARC	6/6	${10}^{-4}$
	MARC	9/4	${10}^{-4}$	1.7-2.0-2.3
${\overline{u}}_2$	ARC	17/17	${10}^{-4}$
	MARC	10/3	${10}^{-4}$	1.9-5.8-8.3

Table 2. Table 2: Solution of the minimization problem ( 33 ) with the one level ARC method and a four level ARC (MARC) (case of n h = 16384 subscript 𝑛 ℎ 16384 n_{h}=16384 ) with u ¯ 1 = 1 ​ rand ​ ( n h , 1 ) subscript ¯ 𝑢 1 1 rand subscript 𝑛 ℎ 1 \bar{u}_{1}=1\,\mathrm{rand}(n_{h},1) , u ¯ 3 = 6 ​ rand ​ ( n h , 1 ) subscript ¯ 𝑢 3 6 rand subscript 𝑛 ℎ 1 \bar{u}_{3}=6\,\mathrm{rand}(n_{h},1) . i ​ t T 𝑖 subscript 𝑡 𝑇 it_{T} denotes the average number of iterations over ten simulations, i ​ t f 𝑖 subscript 𝑡 𝑓 it_{f} the average number of iterations in which the fine level model has been used, 𝚁𝙼𝚂𝙴 𝚁𝙼𝚂𝙴 {\tt RMSE} the root-mean square error with respect to the true solution and save the ratio between the total number of floating point operations required for the Cholesky factorizations in ARC and MARC, respectively.

			$n_h=16384$
$u_0$	Method	$i{t}_T/i{t}_f$	RMSE	save
${\overline{u}}_1$	ARC	6/6	${10}^{-5}$
	MARC	12/3	${10}^{-5}$	1.5-2.0-2.5
${\overline{u}}_3$	ARC	FAIL	FAIL
	MARC	18/5	${10}^{-5}$	-

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x)$$

$$m_{q,k}(x_k,s)=T_{q,k}(x_k,s).$$

$$T_{q,k}(x_k,s)+\frac{{\lambda }_k}{q+1}\Vert s{\Vert }^{q+1},$$

$${\rho }_k=\frac{ared}{pred}=\frac{f(x_k)-f(x_k+s_k)}{m_{q,k}(x_k)-m_{q,k}(x_k,s_k)}.$$

$$\lambda_{k+1}=\bigg{\{}\begin{array}[]{ll}\max\{\lambda_{\min},\gamma_{2}\lambda_{k}\},&\text{ if }\rho_{k}\geq\eta_{2},\\ \max\{\lambda_{\min},\gamma_{1}\lambda_{k}\},&\text{ if }\rho_{k}<\eta_{2},\\ \end{array}$$

$$\Vert {\nabla }_s{m}_{q,k}(x_k,s_k)+{\lambda }_k\Vert s_k{\Vert }^{q-1}{s}_k\Vert \le \theta \Vert s_k{\Vert }^q,$$

$$T(u,v,w)$$

$$T(v,w)(i)$$

$$T(u_1,\dots ,u_i)$$

$$T(u_1,\dots ,u_{i-1})(j_1)$$

$${\nabla }_s^i{m}_{q,k}^H(x_{0,k}^H,\underset{i\mathrm{times}}{\underbrace{s^H,\dots ,s^H}})=[\mathcal{R}({\nabla }_x^i{f}^h(x_k^h))](\underset{i\mathrm{times}}{\underbrace{s^H,\dots ,s^H}}),i=1,\dots ,q,$$

$$[\mathcal{R}({\nabla }_x^i{f}^h(x_k^h))](s_1^H,\dots ,s_i^H)$$

$$⟨[\mathcal{R}({\nabla }_x^i{f}^h(x_k^h))](s_1^H,\dots ,s_{i-1}^H),s_i^H⟩$$

$${\nabla }_s{m}_{q,k}^H(x_{0,k}^H{)}^T{s}^H=(R{\nabla }_x{f}^h(x_k^h){)}^T{s}^H,(s^H{)}^T{\nabla }_x^2{m}_{q,k}^H(x_{0,k}^H)s^H=(s^H{)}^{T}R{\nabla }_x^2{f}^h(x_k^h)P{s}^H.$$

$$m_{q,k}^H(x_{0,k}^H,s^H)=f^H(x_{0,k}^H+s^H)+\sum _{i=1}^q\frac{1}{i!}[\mathcal{R}({\nabla }_x^i{f}^h(x_k^h))-{\nabla }_x^i{f}^H(x_{0,k}^H)](\underset{i\mathrm{times}}{\underbrace{s^H,\dots ,s^H}}),$$

$$m_{2,k}^H(x_{0,k}^H,s^H)=$$

$$+$$

$$\Vert R{\nabla }_x{f}^h(x_k^h)\Vert \ge {\kappa }_H\Vert {\nabla }_x{f}^h(x_k^h)\Vert \text{ and }\Vert R{\nabla }_x{f}^h(x_k^h)\Vert >{ϵ}_H,$$

$$m_{q,k}^H(x_{0,k}^H,s^H)+\frac{{\lambda }_k}{q+1}\Vert s^H{\Vert }^{q+1}$$

$$m_{q,k}^h(x_k^h,s_k^h)=\{\begin{array}{ll}{T}_{q,k}^h(x_k^h,s_k^h)& \text{(Taylor model)},\\ m_{q,k}^H(R{x}_k^h,s_k^H),s_k^h=P{s}_k^H& \text{(lower level model)}.\end{array}$$

$${\rho }_k=\frac{f^h(x_k^h)-f^h(x_k^h+s_k^h)}{m_{q,k}^h(x_k^h)-m_{q,k}^h(x_k^h,s_k^h)},$$

$$m_{q,k}^h(x_k^h)-m_{q,k}^h(x_k^h,s_k^h)=\{\begin{array}{ll}{T}_{q,k}^h(x_k^h)-T_{q,k}^h(x_k^h+s_k^h),& \text{(Taylor model)},\\ m_{q,k}^H(R{x}_k^h)-m_{q,k}^H(R{x}_k^h,s_k^H),& \text{(lower level model)}.\end{array}$$

$$m_{q,k}^{l-1}(R_l{x}_k^l,s^{l-1})$$

$$+\sum _{i=1}^q\frac{1}{i!}[\mathcal{R}({\nabla }_x^i{f}^l(x_k^l))-{\nabla }_x^i{f}^{l-1}(R_l{x}_k^l)](\underset{i\mathrm{times}}{\underbrace{s^{l-1},\dots ,s^{l-1}}}).$$

$$\lambda_{k+1}^{l}=\bigg{\{}\begin{array}[]{ll}\max\{\lambda_{\min},\gamma_{2}\lambda_{k}^{l}\},&\text{ if }\rho_{k}^{l}\geq\eta_{2},\\ \max\{\lambda_{\min},\gamma_{1}\lambda_{k}^{l}\},&\text{ if }\rho_{k}^{l}<\eta_{2}\\ \end{array}$$

$$\Vert {\nabla }^q{f}^h(x)-{\nabla }^q{f}^h(y){\Vert }_T$$

$$\Vert {\nabla }^q{f}^H(x)-{\nabla }^q{f}^H(y){\Vert }_T$$

$$\Vert H{\Vert }_T\stackrel{\mathrm{def}}{=}\underset{\Vert u_1\Vert =\cdots =\Vert u_q\Vert =1}{\max }|H(u_1,\dots ,u_q)|$$

$$g(x+s)=$$

$$|g(x+s)-T_q(x,s)|\le \frac{L}{q}\Vert s{\Vert }^{q+1},$$

$$\Vert \nabla g(x+s)-{\nabla }_s{T}_q(x,s)\Vert \le L\Vert s{\Vert }^q.$$

$$T_{q,k}^h(x_k^h)-T_{q,k}^h(x_k^h,s_k^h)$$

$$m_{q,k}^H(x_{0,k}^H)-m_{q,k}^H(x_{0,k}^H,s_k^H)$$

$$\Vert {\nabla }_s{T}_{q,k}^h(x_k^h,s_k^h)+{\lambda }_k\Vert s_k^h{\Vert }^{q-1}{s}_k^h\Vert$$

$$\Vert {\nabla }_s{m}_{q,k}^H(x_{0,k}^H,s_k^H)+{\lambda }_k\Vert s_k^H{\Vert }^{q-1}{s}_k^H\Vert$$

Full text

1 Introduction

We propose a new family of high-order multilevel optimization methods for unconstrained minimization. Exploiting ideas stemming from multilevel methods allows us to reduce the cost of the step computation, which represents the major cost per iteration of the standard single level procedures. We have been mainly inspired by two driving ideas: the use of high-order models in optimization as introduced in [2], and the multilevel recursive strategy proposed in [13].

When solving unconstrained minimization problems, quadratic models are widely used. These are usually regularized by a quadratic term. For example, trust-region methods have been widely studied and used to globalize Newton-like iterations [12, 22]. Lately in the literature, a different option has received a growing attention: the use of a cubic overestimator of the objective function as a regularization technique for the computation of the step from one iterate to the next, giving rise to quadratic models with cubic regularization. This idea first appeared in [14] and then was reconsidered in [21], where the authors proved that the method has a better worst-case complexity bound compared to standard trust-region methods. Later, in [8, 10], an adaptive variant of the method has been proposed, based on a dynamical choice of the regularization parameters and on an approximate solution of the subproblems. The resulting method is known as adaptive method based on cubic regularization (ARC) and is shown to preserve the attractive global complexity bound established in [21]. In recent years the method has attracted further interest, see for example [9, 23, 26].

In recent publications also methods of higher order start to gain interest, see for example [2, 24]. In [2] in particular, it has been observed that the good complexity bound of ARC can be made even lower, if one is willing to use higher-order derivatives. In specific applications this computation is indeed feasible, for example when considering partially separable functions [11]. The authors in [2] present a family of methods that generalizes ARC, and that uses high-order regularized models. Specifically, they are based on models of order $q\geq 1$, regularized by a term of order $q+1$. ARC belongs to this family and corresponds to the choice $q=2$. The authors in [2] propose a unifying framework to describe the theoretical properties of the methods in this class. It is proved that the method based on the $q$-th order model requires at most $O\left(\epsilon^{-\frac{q+1}{q}}\right)$ function evaluations to find a first-order critical point, where $\epsilon$ denotes the absolute accuracy level.

However, the use of higher-models come along with higher computational costs. The main cost per iteration of the methods described in [2] is represented by the step computation through the model minimization. This cost is proportional to the dimension of the problem, it can therefore be significant for large-scale problems. For second-order models this issue has been faced for example in [13], where the authors exploit ideas coming from multigrid [15] to reduce the cost of the minimization. Indeed, the idea of making use of more grids to solve a large-scale problem has been extended also to optimization, see for example [13, 16, 17, 18, 19, 20, 25]. These methods share with classical multigrid methods the idea of exploiting a hierarchy of problems (in this case a sequence of nonlinear functions) defined on lower dimensional spaces, approximating the original objective function $f$. The simplified expressions of the objective function are used to build models that are cheaper to minimize, and are used to define the step. Specifically, in [13], the authors present an extension of classical multigrid methods for nonlinear optimization problems, [4, 5] or [6, Ch. 3], to a class of multilevel trust-region based optimization algorithms.

Our contributions

Inspired by the ideas presented in [2, 13], we propose a family of multilevel optimization methods using high-order regularized models that generalizes the methods proposed in both papers. The aim is to decrease the computational cost of the methods in [2], extending the ideas in [13] to higher-order models. We also develop a theoretical analysis for the resulting family of methods. The main theoretical results are provided in Theorems 1, 2 and 3, respectively. In these theorems we successively prove the global convergence property of the methods, then evaluate a worst-case complexity bound to reach a first-order critical point and finally provide local convergence rates. The global convergence analysis generalizes the results in [2] and appears as much simpler than that in [13]. Moreover we establish local convergence results towards second-order stationary points, that are not present neither in [2] nor in [13]. These results not only generalize those in [27], that are valid only for $q=2$, but also apply to the one level methods in [2]. From a practical point of view, we implemented the method of the family corresponding to $q=2$. This represents a multilevel version of ARC method.

To the best of our knowledge, this is the first time that multilevel optimization strategies, based on models of generic order $q\geq 1$, are proposed, and that a unifying framework is introduced to study their convergence. Moreover, multilevel versions of ARC have never been analysed or tested numerically before.

The manuscript is organized as follows. In Section 2, we briefly introduce the family of optimization methods using high-order regularized models considered in [2]. Section 3 and Section 4 represent our main contribution. We introduce in Section 3 the multilevel extensions of the methods presented in Section 2, and we provide a theoretical analysis in Section 4. Specifically, we focus on global convergence in Section 4.1, worst-case complexity in Section 4.2 and local convergence in Section 4.3. In Section 5 we then present results related to numerical experiments performed with the multilevel method corresponding to $q=2$. Finally, conclusions are drawn in Section 6.

2 High-order iterative optimization methods

Let $q\geq 1$ be an integer. Let us consider a minimization problem of the form:

[TABLE]

with $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ a bounded below and $q$-times continuously differentiable function, called the objective function.

Classical iterative optimization methods for unconstrained minimization are based on the use of a model to approximate the objective function at each iteration. In this section, we describe the iterative optimization methods using high-order models presented in [2].

2.1 Model definition and step acceptance

At each iteration $k$, given the current iterate $x_{k}$, the objective function is approximated by the Taylor series $T_{q,k}$ of $f(x_{k}+s)$ (with $s\in\mathbb{R}^{n}$) truncated at order $q$. The Taylor model of order $q$ denoted as $m_{q,k}$ is then defined as:

[TABLE]

A step $s_{k}$ is then found minimizing (possibly approximately) the regularized model

[TABLE]

where $\lambda_{k}$ is a positive value called regularization parameter. The step $s_{k}$ is used to define a trial point i.e. $x_{k+1}=x_{k}+s_{k}$. At each iteration, it has to be decided whether to accept the step or not. This decision is based on the accordance between the decrease in the function and in the model. More precisely, at each iteration both the decrease achieved in the model, that we call predicted reduction, $pred=m_{q,k}(x_{k})-m_{q,k}(x_{k},s_{k})$, and that achieved in the objective function, that we call actual reduction, $ared=f(x_{k})-f(x_{k}+s_{k})$, are computed. The step acceptance is then based on the ratio:

[TABLE]

If the model is accurate, $\rho_{k}$ will be close to one. Then, the step $s_{k}$ is accepted if $\rho_{k}$ is larger than or equal to a chosen threshold $\eta_{1}\in(0,1)$ and is rejected otherwise. In the first case, the step is said to be successful, and otherwise the step is unsuccessful.

After the step acceptance, the regularization parameter is updated for the next iteration. The update is still based on the ratio (4). If the step is successful, the regularization parameter is decreased, otherwise it is increased. The whole procedure is stopped when a minimizer of $f$ is reached. Usually, the stopping criterion is based on the norm of the gradient, i.e. given an absolute accuracy level $\epsilon>0$ the iterations are stopped as soon as $\|\nabla_{x}f(x_{k})\|<\epsilon$. The whole procedure is sketched in Algorithm 1.

2.2 Minimization of the model

The main computational work per iteration in this kind of methods is represented by the minimization of the regularized model (3). This is the most expensive task, and the cost naturally depends on the dimension of the problem. However, from the convergence theory of such methods, it is well known that it is not necessary to minimize the model exactly to get a globally convergent method.

A well-known possibility is to minimize the model until the Cauchy decrease is achieved, i.e. until a fraction of the decrease provided by the Cauchy step (the step that minimizes the model in the direction of the negative gradient) is obtained. In [2] the authors consider a different stopping criterion for the inner iterations, the one originally proposed in [8, 10], which has the advantage of allowing for simpler convergence proofs. The inner iterations are stopped as soon as the norm of the gradient of the regularized model becomes lower or equal than a multiple of the power $q$ of the norm of the step $s_{k}$:

[TABLE]

for a chosen constant $\theta>0$.

For very large-scale problems however, even an approximate minimization of (3) may be really costly. Then, in the next section we propose multilevel variants of the procedures, that rely on simplified models of the objective function, cheaper to optimize, allowing to reduce the global cost of the optimization procedure.

3 Multilevel optimization methods

We describe the multilevel extension of the family of methods presented in Section 2. The procedures are inspired by the multilevel trust-region approach presented in [13], where only second-order models with quadratic regularization have been considered. Here, we generalize this approach by allowing also higher-order models, i.e. $q>2$.

3.1 Preliminaries and notations

In standard optimization methods the minimization of (3) represents the major cost per iteration, which crucially depends on the dimension $n$ of the problem. When $n$ is large, the solution cost is therefore often significant. We want to reduce this cost by exploiting the knowledge of alternative simplified expressions of the objective function. More specifically, we assume that we know a collection of functions $\{f_{l}\}_{l=1}^{l_{\max}}$ such that each $f_{l}$ is a $q$-times continuously differentiable function from $\mathbb{R}^{n_{l}}\rightarrow\mathbb{R}$ and $f^{l_{\max}}(x)=f(x)$ for all $x\in\mathbb{R}^{n}$. We will also assume that, for each $l=2,\dots,{l_{\max}}$, $f_{l}$ is more costly to minimize than $f_{l-1}$. This is the typical scenario when the problem arises from the discretization of an infinite dimensional problem and $f_{l}$ represent increasingly finer discretizations. In this case, $n_{l}\geq n_{l-1}$ for all $l$, but of course this is not the only possible application. As we do not assume the hierarchy to come from a discretization process, we do not use the terminology typically used in the field of multigrid methods. We will then use ’levels’ rather than ’grids’.

The methods we propose are recursive procedures, so it suffices to describe the two-level case. Then, for sake of simplicity, from now on, we will assume that we have just two approximations to our objective $f$ at disposal. This amounts to consider $l_{\max}=2$.

For ease of notation, we will denote by $f^{h}:\mathcal{D}^{h}\subseteq\mathbb{R}^{n_{h}}\rightarrow\mathbb{R}$ the approximation at the highest level ($f^{h}(x)=f^{l_{\max}}(x)$ in the notation previously used) and by $f^{H}:\mathcal{D}^{H}\subseteq\mathbb{R}^{n_{H}}\rightarrow\mathbb{R}$ the other approximation available, that is cheaper to optimize. The quantities on the highest level will be denoted by a superscript $h$, whereas the quantities on the lower level will be denoted by a superscript $H$. Let $x_{k}^{h}$ denote the $k$-th iteration at the highest level. In the following $\|\cdot\|$ will denote the Euclidean norm. We will use the same notation for all the spaces we will consider, the space on which the norm is defined will be clear by the context.

Remark 1**.**

To deal with high order derivatives, and to properly handle the concept of coherence between lower and higher level model, we will need to use a tensor notation, that we introduce here for convenience of the reader, see [3]. We first consider a tensor of order three, and then extend the definition to a tensor of order $i\in\mathbb{N}$.

Definition 1**.**

Let $T\in\mathbb{R}^{n\times n\times n}$, and $u,v,w\in\mathbb{R}^{n}$. Then $T(u,v,w)\in\mathbb{R}$, $T(u,v)\in\mathbb{R}^{n}$ and

[TABLE]

Definition 2**.**

Let $i\in\mathbb{N}$ and $T\in\mathbb{R}^{n^{i}}$, and $u_{1},\dots,u_{i}\in\mathbb{R}^{n}$. Then $T(u_{1},\dots,u_{i})\in\mathbb{R}$, $T(u_{1},\dots,u_{i-1})\in\mathbb{R}^{n}$ and

[TABLE]

3.2 Construction of the lower level model

The main idea is to use $f^{H}$ to construct, in the neighbourhood of the current iterate, an alternative model $m_{q,k}^{H}$ to the Taylor model $m_{q,k}^{h}$ in (2) for $f^{h}=f$ [13]. The alternative model $m_{q,k}^{H}$ should be cheaper to optimize than $m_{q,k}^{h}$, and will be used, whenever suitable, to define the step. Of course, for $f^{H}$ to be useful at all in minimizing $f^{h}$, there should be some relation between the variables of these two functions. We henceforth assume the following.

Assumption 1**.**

Let us assume that there exist two full-rank linear operators $R:\mathbb{R}^{n_{h}}\rightarrow\mathbb{R}^{n_{H}}$ and $P:\mathbb{R}^{n_{H}}\rightarrow\mathbb{R}^{n_{h}}$ such that $P=\alpha R^{T}$, for a fixed scalar $\alpha>0$. Let us assume also that it exists $\kappa_{R}>0$ such that $\max\{\|R\|,\|P\|\}\leq\kappa_{R}$, where $\|\cdot\|$ denotes the matrix norm induced by the Euclidean norm at the fine level.

In the following, we can assume $\alpha=1$, without loss of generality, as the problem can be easily scaled to handle the case $\alpha\neq 1$.

At each iteration $k$ at highest level we set $x_{0,k}^{H}=R~{}x_{k}^{h}$, i.e. the initial iterate at the lower level is set as the projection of the current iterate, and we define the lower level model $m_{q,k}^{H}$ as a modification of the coarse function $f^{H}$. Given $q$, $f^{H}$ is modified adding $q$ correction terms, to enforce the following relation:

[TABLE]

where $\mathcal{R}(\nabla^{i}_{x}f^{h}(x_{k}^{h}))$ is such that for all $i=1,\dots,q$ and $s_{1}^{H},\dots,s_{i}^{H}\in\mathbb{R}^{n_{H}}$

[TABLE]

where $\nabla^{i}_{x}f^{h},\nabla^{i}_{s}m_{q,k}^{H}$ denote the $i$-th order tensor of $f^{h}$ and $m_{q,k}^{H}$ respectively, $\langle\cdot,\cdot\rangle$ denotes the scalar product, and for generic $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $s_{1},\dots,s_{i}\in\mathbb{R}^{n}$, $\nabla^{i}g(x,s_{1},\dots,s_{i})$ is the same as $\nabla^{i}g(x)(s_{1},\dots,s_{i})$, which is given in Definition 2.

For instance, if $q=2$, relation (6) simply becomes:

[TABLE]

Relation (6) crucially ensures that the behaviours of $f^{h}$ and $m_{q,k}^{H}$ are coherent up to order $q$ in a neighbourhood of $x_{k}^{h}$ and $x_{0,k}^{H}$. To achieve (6), we define $m_{q,k}^{H}$ as

[TABLE]

with $\mathcal{R}(\nabla^{i}_{x}f^{h}(x_{k}^{h}))$ defined in (7)-(8). When $q=2$ this is simply:

[TABLE]

3.3 Step computation and step acceptance

At each generic iteration $k$ of our method, a step $s^{h}_{k}$ has to be computed to define the new iterate. Then, one has the choice between the Taylor model (2) and a lower level model (9).

Obviously, it is not always possible to use the lower level model. For example, it may happen that $\nabla_{x}f^{h}(x_{k}^{h})$ lies in the nullspace of $R$ and thus that $R\nabla_{x}f^{h}(x_{k}^{h})$ is zero while $\nabla_{x}f^{h}(x_{k}^{h})$ is not. In this case, the current iterate appears to be first-order critical for $m_{q,k}^{H}$ while it is not for $f^{h}$. Using the model $m_{q,k}^{H}$ is hence potentially useful only if $\|\nabla_{s}m_{q,k}^{H}(x_{0,k}^{H})\|=\|R\nabla_{x}f^{h}(x_{k}^{h})\|$ is large enough compared to $\|\nabla_{x}f^{h}(x_{k}^{h})\|$ [13]. We therefore restrict the use of the model $m_{q,k}^{H}$ to iterations where

[TABLE]

for some constant $\kappa_{H}\in(0,\min\{1,\|R\|\})$ and where $\epsilon_{H}\in(0,1)$ is a measure of the first-order criticality for $m_{q,k}^{H}$ that is judged sufficient at level $H$ [13]. Note that, given $\nabla_{x}f^{h}(x_{k}^{h})$ and $R$, this condition is easy to check before even attempting to compute a step at a lower level.

If the Taylor model is chosen, then we just compute a step as in standard methods, minimizing (possibly approximately) the corresponding regularized model (3). If the lower level model is chosen, we then minimize the following regularized model:

[TABLE]

(possibly approximately) and obtain a point $x_{*,k}^{H}$ such that (if the minimization is successful) the value of the regularized model has been reduced, and a step $s_{k}^{H}=x_{*,k}^{H}-x_{0,k}^{H}$ (note that the iteration indices always refer to the highest level, we are not indexing the iterations on the lower level for the minimization of the lower level model). This step has to be prolongated back on the fine level, i.e. we define $s_{k}^{h}=Ps_{k}^{H}$.

Then, $m_{q,k}^{h}$ will be defined as:

[TABLE]

In both cases, after the step is found, we have to decide whether to accept it or not. The step acceptance is based on the ratio:

[TABLE]

where we remind that, from (12), the denominator is defined as:

[TABLE]

As in the standard form of the methods, the step is accepted if it provides a sufficient decrease in the function, i.e. if given $\eta_{1}>0$, $\rho_{k}\geq\eta_{1}$. The regularization parameter is also updated as in Algorithm 1. We sketch the whole procedure in Algorithm 2.

Some comments are necessary to explain Step 6 in Algorithm 2. The generic framework sketched in Algorithm 2 comprises different possible methods. Specifically, one of the flexible features (inherited by the method in [13]) is that, to ensure convergence, the minimization at lower levels can be stopped after the first successful iteration, as we will see in the next section. This therefore opens the possibility to consider both fixed form recursion patterns and free form ones. A free form pattern is obtained when Algorithm 2 is run carrying the minimization at each level out, until the norm of the gradient becomes small enough. The actual recursion pattern is then uniquely determined by the progress of minimization at each level and may be difficult to forecast. By contrast, the fixed form recursion patterns are obtained by specifying a maximum number of successful iterations at each level, a technique directly inspired from the definitions of V- and W-cycles in multigrid algorithms [15].

4 Convergence theory

In this section, we provide a theoretical analysis of the proposed family of multilevel methods. Inspired by the convergence theory reported in [2], we prove global convergence of the proposed methods to first-order critical points and we provide a worst-case complexity bound to reach such a point, generalizing the theory proposed in [2, 13]. At the same time the proposed analysis also appears as simpler than that in [13], since the regularization parameter $\lambda_{k}$ is directly updated, rather than the trust-region radius, and since we use the stopping criterion (5) as in [2]. The use of this criterion allows for simpler convergence proofs and enables us to concentrate on the multilevel algorithm, that is the main contribution of the paper. Moreover, we also propose local convergence results, which also apply to the methods in [2], and that extend those in [27] to higher-order models.

Note that, as the methods are recursive, we can restrict the analysis to the two-level case. For the analysis we need the following regularity assumptions as in [2].

Assumption 2**.**

Let $f^{h}$ and $f^{H}$ be $q$-times continuously differentiable and bounded below functions. Let us assume that the $q$-th derivative tensors of $f^{h}$ and $f^{H}$ are Lipschitz continuous, i.e. that there exist constants $L_{h},L_{H}$ such that

[TABLE]

where $\|\cdot\|_{T}$ is the tensor norm recursively induced by the Euclidean norm on the space of $q$-th order tensors, which for a tensor $H$ of order $q$ is given by

[TABLE]

where the action of $H$ on $(u_{1},\dots,u_{q})$ is given in Definition (2).

We remind three useful relations, following from Taylor’s theorem, see for example relations (2.3) and (2.4) in [2].

Lemma 1**.**

Let $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a $q$-times continuously differentiable function with Lipschitz continuous $q$-th order tensor, with $L$ the corresponding Lipschitz constant. Given its Taylor series $T_{q}(x,s)$ truncated at order $q$, it holds:

[TABLE]

4.1 Global convergence

In this section we prove the global convergence property of the method. Our analysis proceeds in three steps. First, we bound the quantity $\lvert 1-\rho_{k}\rvert$ to prove that $\lambda_{k}$ must be bounded above. Then, we relate the norm of the step and the norm of the gradient. Finally, we use these two ingredients to conclude proving that the norm of the gradient goes to zero.

4.1.1 Upper bound for the regularization parameter $\lambda_{k}$

At iteration $k$ we either minimize (decrease) the regularized Taylor model (3), or the regularized lower level model (11). Consequently, it respectively holds:

[TABLE]

In both cases, the minimization process is stopped as soon as the stopping condition

[TABLE]

for $\theta>0$ is satisfied, respectively. In both cases we are sure that it will exist a point that satisfies (18), as when the level is selected, a standard one-level optimization method is used, and the analysis in [2] applies.

Let us consider the quantity

[TABLE]

with the denominator defined in (13). If at step $k$ the Taylor model is chosen, from relation (15) applied to $f^{h}$ and (17a) we obtain the inequality:

[TABLE]

If the lower level model is used, we have

[TABLE]

Let us consider the numerator in this expression. From relations (9) and (14) applied to $f^{H}$, using its Taylor series $T_{q,k}^{H}$, it follows

[TABLE]

Similarly the relation (14) applied to $f^{h}$ yields

[TABLE]

From relation (6) we can rewrite $f^{h}(x_{k}^{h})-T_{q,k}^{h}(x_{k}^{h},s_{k}^{h})$ as:

[TABLE]

Then, subtracting (21) from (20), we obtain

[TABLE]

Using Assumption 2, we obtain:

[TABLE]

From relation (17b) we finally obtain:

[TABLE]

Then, in both cases (when either a Taylor model or a lower level model is used), it exists a strictly positive constant $K$ such that the following relation holds:

[TABLE]

Using this last relation and the updating rule of the regularization parameter, we deduce that $\lambda_{k}$ must be bounded above. Indeed, in case of unsuccessful iterations, $\lambda_{k}$ is increased. If $\lambda_{k}$ is increased, the ratio appearing in the right hand side of (22) is progressively decreased, until it becomes smaller than $1-\eta_{1}$. In this case, $\rho_{k}>\eta_{1}$, so a successful step is taken and $\lambda_{k}$ is decreased. Hence $\lambda_{k}$ cannot be greater than

[TABLE]

4.1.2 Relating the steplength to the norm of the gradient

Our next step is to show that the steplength cannot be arbitrarily small, compared to the norm of the gradient of the objective function. If the Taylor model is used, from [2, Lemma 2.3] it follows:

[TABLE]

If the lower level model is chosen, we have:

[TABLE]

By (16), the first term can be bounded by $\kappa_{R}L_{h}\|s_{k}^{h}\|^{q}$. Considering that $s_{k}^{h}=Ps_{k}^{H}$ and $\|P\|\leq\kappa_{R}$, we obtain the upper bound $\kappa_{R}^{2}L_{h}\|s_{k}^{H}\|^{q}$. Regarding the second term, taking into account that from relations $s_{k}^{h}=Ps_{k}^{H}$, $R=P^{T}$, and (8), for all $p^{H}\in\mathbb{R}^{n_{H}}$ it holds:

[TABLE]

we can write

[TABLE]

Then, from

[TABLE]

we obtain

[TABLE]

which represents the Taylor remainder for the approximation of $\nabla_{x}f^{H}$ by $\nabla_{s}T_{q,k}^{H}$. Therefore, by relation (16), this quantity can be bounded above by $L_{H}\|s_{k}^{H}\|^{q}$. The third term, from (18), is less than $\theta\|s_{k}^{H}\|^{q}$. Then, since $\lambda_{k}\leq\lambda_{\max}$, we finally obtain

[TABLE]

4.1.3 Proof of global convergence

Let us consider the sequence of successful iterations ($\rho_{k}\geq\eta_{1}$). They are divided into two groups, $K_{s,f}$ the successful iterations at which the fine model has been employed and $K_{s,l}$ the ones at which the lower level model has been employed. Let us define $k_{1}$ the index of the first successful iteration. We remind that at successful iterations $\rho_{k}\geq\eta_{1}$. Due to the updating rule of the regularization parameter in Algorithm 2 we have $\lambda_{k}\geq\lambda_{\min}$. Hence from relations (13), (17), (24) and (26), (10) it follows that:

[TABLE]

Hence we conclude that $\sum_{K_{s,f}\cup K_{s,l}}\|\nabla_{x}f^{h}(x_{k}^{h}+s_{k}^{h})\|$ is a bounded series and therefore has a convergent subsequence. Then, $\|\nabla_{x}f^{h}(x_{k}^{h}+s_{k}^{h})\|$ converges to zero on the subsequence of successful iterations.

We can then state the global convergence property towards first-order critical points in the following theorem.

Theorem 1**.**

Let Assumptions 1 and 2 hold. Let $\{x_{k}^{h}\}$ be the sequence of fine level iterates generated by Algorithm 2. Then, $\{\|\nabla_{x}f^{h}(x_{k}^{h})\|\}$ converges to zero on the subsequence of successful iterations.

4.2 Worst-case complexity

We now want to evaluate the worst-case complexity of our methods, to reach a first order stationary point. We assume then that the procedure is stopped as soon as $\|\nabla_{x}f^{h}(x_{k}^{h})\|\leq\epsilon$ for $\epsilon>0$. The proof is similar to that of Theorem 2.5 in [2].

To evaluate the complexity of the proposed methods, we have to bound the number of successful and unsuccessful iterations performed before the stopping condition is met. Let us then define $k_{f}$ the index of the last iterate for which $\|\nabla_{x}f^{h}(x_{k}^{h})\|>\epsilon$, $K_{s}=\{0<j\leq k_{f}\,|\,\rho_{j}\geq\eta_{1}\}$ the set of successful iterations before iteration $k_{f}$, and $K_{u}$ its complementary in $\{1,\dots,k_{f}\}$. We can use the same reasoning as that used to derive (27), but considering in the sum just the successful iterates in $K_{s}$. Remind that before termination $\|\nabla_{x}f^{h}(x_{k}^{h})\|>\epsilon$ and, in case the lower level model is used, $\|R\nabla f^{H}(x_{k}^{h})\|>\kappa_{H}\|\nabla f^{H}(x_{k}^{h})\|>\kappa_{H}\epsilon$ (otherwise at that iteration the Taylor model would have been used). It then follows:

[TABLE]

from which we get the desired bound on the total number of successful iterations. We can then bound the cardinality of $K_{u}$, with respect to the cardinality of $K_{s}$. From the updating rule of the regularization parameter, it holds:

[TABLE]

Then, proceeding inductively, we conclude that:

[TABLE]

Then,

[TABLE]

and, given that $\gamma_{1}<1$, we obtain:

[TABLE]

We can then state the following result.

Theorem 2**.**

Let Assumptions 1 and 2. Let $f_{low}$ denote a lower bound on $f$ and let $k_{1}$ denote the index of the first successful iteration in Algorithm 2. Then, given an absolute accuracy level $\epsilon>0$, Algorithm 2 needs at most

[TABLE]

iterations in total to produce an iterate $x_{k}^{h}$ such that $\|\nabla_{x}f(x_{k})\|\leq\epsilon$, where

[TABLE]

with $K_{1}$ and $K_{2}$ defined in (24), (26), $\gamma_{1},\gamma_{3},\lambda_{0},\lambda_{\min}$ defined in Algorithm 2 and $\lambda_{\max}$ defined in (23).

Theorem 2 reveals that the use of lower level steps does not deteriorate the complexity of the method, and that the complexity bound $O(\epsilon^{-\frac{q+1}{q}})$ is preserved. This is a very satisfactory result, because each iteration of the multilevel methods will be less expensive than one iteration of the corresponding one-level method, thanks to the use of the cheaper lower level models. Consequently, if the number of iterations in the multilevel strategy is not increased, we can expect global computational savings.

4.3 Local convergence

In this section we study the local convergence of the proposed methods towards second-order stationary points. We assume $q\geq 2$ in this section, otherwise the problem is not well-defined. Thanks to the use of high order models, our methods are expected to attain a fast local convergence rate, especially for growing $q$. The results reported here are inspired by [27] and extend the analysis proposed therein.

We denote by $\mathcal{X}$ the set of second-order critical points of $f$, i.e. of points $x^{*}$ satisfying the second-order necessary conditions:

[TABLE]

i.e. $\nabla_{x}^{2}f(x^{*})$ is a symmetric positive semidefinite matrix. We denote by $\mathcal{B}(x,\rho)=\{y\,s.t.\,\|y-x\|\leq\rho\}$ and for all $x\in\mathbb{R}^{n}$,  $\mathcal{L}(f(x))=\{y\in\mathbb{R}^{n}\,|\,f(y)\leq f(x)\}$ for $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$.

Remark 2**.**

From the assumption that $f$ is a $q$ times continuously differentiable function, it follows that its $i$-th derivative tensor is locally Lipschitz continuous for all $i\leq q-1$.

Following [27], we first prove an intermediate lemma that allows us to relate, at generic iteration $k$, the norm of the step and the distance of the current iterate from the space of second-order stationary points. This lemma holds without need of assuming a stringent non-degeneracy condition, but rather under a local error bound condition, which is a much weaker requirement as it can be satisfied also when $f$ has non isolated second-order critical points.

Assumption 3**.**

There exist strictly positive scalars $\kappa_{EB},\rho>0$ such that

[TABLE]

where $\mathcal{X}$ is the set of second-order critical points of $f$, $\mathrm{dist}(x,\mathcal{X})$ denotes the distance of $x$ to $\mathcal{X}$ and $\mathcal{N}(\mathcal{X},\rho)=\{x\;|\;\mathrm{dist}(x,\mathcal{X})\leq\rho\}$.

This condition has been proposed for the first time in [27]. It is different from other error bound conditions in the literature as, in contrast to them, $\mathcal{X}$ is not the set of first-order critical points, but of second-order-critical points. In addition to being useful for proving convergence, it is also interesting on its own, as it is shown to be equivalent to a quadratic growth condition ([27, Theorem 1]) under mild assumptions on $f$.

Lemma 2**.**

Let Assumptions 1 and 2 hold. Let $\{x_{k}^{h}\}$ be the sequence generated by Algorithm 2 and $x_{k}^{*}$ be a projection point of $x_{k}^{h}$ onto $\mathcal{X}$. Assume that it exists a strictly positive constant $\underline{\rho}$ such that $\{x_{k}^{h}\}\in\mathcal{B}(x_{k}^{*},\underline{\rho})$ and that $\nabla_{x}^{2}f$ is Lipschitz continuous in $\mathcal{B}(x_{k}^{*},\underline{\rho})$ with Lipschitz constant $L_{2}$. Then, it holds:

[TABLE]

with

[TABLE]

with $\lambda_{\max}$ and $\kappa_{R}$ defined respectively in (23) and Assumption 1.

The proof of the lemma is reported in the appendix. This lemma can be used to prove that if it exists an accumulation point of $\{x_{k}^{h}\}$ that belongs to $\mathcal{X}$, then the full sequence converges to that point and that the rate of convergence depends on $q$. First, we can prove that the set of accumulation points is not empty.

Lemma 3**.**

Let Assumptions 1 and 2 hold. Let $\{x_{k}^{h}\}$ be the sequence of fine level iterates generated by Algorithm 2. If $\mathcal{L}(f(x_{k}^{h}))$ is bounded for some $k\geq 0$, then the sequence has an accumulation point that is a first-order stationary point.

Proof.

As $\{f(x_{k}^{h})\}$ is a decreasing sequence, and $\mathcal{L}(f(x_{k}^{h}))$ is bounded for some $k\geq 0$, $\{x_{k}^{h}\}$ is a bounded sequence and it has an accumulation point. From Theorem 1, all the accumulation points are first-order stationary points. ∎

Theorem 3**.**

Let Assumptions 1 and 2 hold. Let $\{x_{k}^{h}\}$ be the sequence of fine level iterates generated by Algorithm 2. Assume that $\mathcal{L}(f(x_{k}^{h}))$ is bounded for some $k\geq 0$ and that it exists an accumulation point $x^{*}$ such that $x^{*}\in\mathcal{X}$. Then, the whole sequence $\{x_{k}^{h}\}$ converges to $x^{*}$ and it exist strictly positive constants $c\in\mathbb{R}$ and $\bar{k}\in\mathbb{N}$ such that:

[TABLE]

Proof.

As $x^{*}$ is an accumulation point, we have that $\lim\limits_{k\rightarrow\infty}\mathrm{dist}(x_{k}^{h},\mathcal{X})=0$. Then, it exist $\rho$ and $k_{1}$ such that $x_{k}^{h}\in\mathcal{N}(\mathcal{X},\rho)$ for all $k\geq k_{1}$. Therefore, from Assumption 3 it holds

[TABLE]

Moreover, from Remark 2, $\nabla_{x}^{2}f$ is locally Lipschitz continuous, so Lemma 2 applies to all $k\geq k_{1}$.

Let us first consider the case in which the Taylor model is employed. It follows from (31), (24) and (29) that for all $k\geq k_{1}$

[TABLE]

If the lower level model is employed, from (31), (10), (26) and (36) it follows that for all $k\geq k_{1}$

[TABLE]

Then in both cases, it exists $\bar{C}$ such that

[TABLE]

where

[TABLE]

With this result, we can prove the convergence of $\{x_{k}^{h}\}$ with standard arguments. We repeat for example the arguments of the proof of [27, Theorem2] for convenience. Let $\eta>0$ be an arbitrary value. As $\lim_{k\rightarrow\infty}\mathrm{dist}(x_{k}^{h},\mathcal{X})=0$, it exists $k_{2}\geq 0$ such that

[TABLE]

Then,

[TABLE]

From (29), it then holds for all $k\geq\bar{k}$ and $j\geq 0$:

[TABLE]

i.e. that $\{x_{k}^{h}\}_{k\geq\bar{k}}$ is a Cauchy sequence and so the whole sequence is convergent. Finally we establish the $q$-th order rate of convergence of the sequence. For any $k\geq\bar{k}$,

[TABLE]

Combining this with $\mathrm{dist}(x_{k}^{h},\mathcal{X})\leq\|x_{k}^{h}-x^{*}\|$, and setting $c=2C\bar{C}$ we obtain the thesis (30).

Therefore $\{x_{k}^{h}\}$ converges at least with order $q$ to $x^{*}$. ∎

5 Numerical results

In this section, we report on the practical performance of a method in the family.

We have implemented the method corresponding to $q=2$ in Algorithm 2 in Julia [1] (version 0.6.1). This is a multilevel extension of the method AR$2$ in Algorithm 1, which is better known as ARC [8, 10]. We will therefore denote the implemented multilevel method as MARC (multilevel adaptive method based on cubic regularization), rather than MAR$2$.

We consider the following two-dimensional nonlinear problem in the unit square domain $S_{2}$:

[TABLE]

where $g$ is obtained such that the analytical solution to this problem is given by

[TABLE]

The negative Laplacian operator is discretized using finite difference, giving a symmetric positive definite matrix $A$, that also takes into account the boundary conditions. The discretized version of the problem is then a system of the form $Au+e^{u}=g$, where $u$, $g$, $e^{u}$ are vectors in $\mathbb{R}^{n_{h}}$, in which the columns of matrices $U$, $G$, $E$ are stacked, with $U_{i,j}=u(x_{i},y_{j})$, $G_{i,j}=g(x_{i},y_{j})$, $E_{i,j}=e^{u(x_{i},y_{j})}$, for $x_{i},y_{j}$ grid points, $i,j=1,\dots,\sqrt{n_{h}}$.

The MARC algorithm is then used on the nonlinear minimization problem

[TABLE]

which is equivalent to the system $Au+e^{u}=g$. The coarse approximations to the objective function arise from a coarser discretization of the problem. Each coarse two-dimensional grid has a dimension that is four times lower than the dimension of the grid on the corresponding upper level.

The prolongation operators $P_{\ell}$ from level $\ell-1$ to $\ell$ are based on the nine-point interpolation scheme defined by the stencil $\begin{pmatrix}\frac{1}{4}\;\frac{1}{2}\;\frac{1}{4}\\ \frac{1}{2}\;1\;\frac{1}{2}\\ \frac{1}{4}\;\frac{1}{2}\;\frac{1}{4}\end{pmatrix}$ and the full weighting operators defined by $R_{\ell}=\frac{1}{4}P_{\ell}^{T}$ are used as restriction operators [6].

We compare the one-level ARC with MARC. Parameters common to both methods are set as: $\epsilon^{l_{\max}}=10^{-7}$, $\gamma_{1}=0.85$, $\gamma_{2}=0.5$, $\gamma_{3}=2$, $\lambda_{0}=0.05$ $\eta_{1}=0.1$, $\eta_{2}=0.75$. For MARC we set $\kappa=0.1$ and $\epsilon^{\ell}=\epsilon^{\ell_{\max}}$ for all $\ell$.

At each iteration we find an approximate minimizer of the cubic models as described in [8, §6.2]. This requires a sequence of Cholesky factorizations, which represents the dominant cost per nonlinear iteration. We measure the performance of the methods in terms of total number of floating point iterations required for these factorizations.

We study the effect of the multilevel strategy on the convergence of the method for problems of fixed dimension $n_{h}$. We then consider the solution of problem (33) using two different discretizations with $n_{h}=4096$ (Table 1) and $n_{h}=16384$ (Table 2), respectively. We allow $4$ levels in MARC. We report the results of the average of ten simulations with different random initial guesses of the form $u_{0}=a\,\mathrm{rand}(n_{h},1)$, for different values of $a$. In each simulation the random starting guess is the same for the two considered methods. All the quantities reported in Tables 1 and 2 are the average of the values obtained over the ten simulations. $it_{T}$ denotes the number of total iterations, $it_{f}$ denotes the number of iterations in which the Taylor model has been used, ${\tt RMSE}$ is the root-mean square error with respect to the true solution and save is the ratio between the total number of floating point operations required for the Cholesky factorizations by ARC and MARC, respectively. For the save quantity, we report three values: the minimum, the average and the maximum value obtained over the ten simulations.

The results reported in Tables 1 and 2 confirm the relevance of MARC as compared to ARC. The numerical experiments highlight the different convergence properties of both algorithms. The use of MARC is especially convenient when the initial guess is not so close to the true solution. Indeed, the performance of ARC deteriorates as the distance of the initial guess from the true solution increases, while MARC seems to be much less sensible to this choice. For the problem of smaller dimension, ARC still manages to find a solution for further initial guesses, even if this requires an higher number of iterations, while for the problem of larger dimension the method fails to find a solution in feasible time. The new multilevel approach is found to lead to considerable computational savings in terms of floating point operations compared to the classical one-level strategy.

6 Conclusions

We have introduced a family of multilevel methods of order $q\geq 1$ for unconstrained minimization. These methods represent an extension of the higher-order methods presented in [2] and of the multilevel trust-region method proposed in [13]. We have proposed a unifying framework to analyse these methods, which is useful to prove their convergence properties and evaluate their worst-case complexity to reach first-order stationary points. As expected, we show that the local rate of convergence and the complexity bound depend on $q$ and high values of $q$ allow both fast local convergence and lower complexity bounds.

We believe this represents a contribution in the optimization field, as the use of multilevel ideas allows to reduce the major cost per iteration of the high-order methods. This gives a first answer to the question posed in [2] about whether the approach presented there can have practical implications, in applications for which computing $q$ derivatives is feasible.

We have implemented the multilevel method corresponding to $q=2$ and presented numerical results that show the considerable benefits of the multilevel strategy in terms of savings in floating point operations. Additional numerical results can be found in [7], where the authors apply the multilevel method in the family corresponding to $q=1$ to problems arising in the training of artificial neural networks for the approximate solution of partial differential equations. This case is particularly interesting as it allows to show the efficiency of multilevel methods even for problems without an underlying geometrical structure.

Appendix A Proof of Lemma 4.5

In this appendix we report the proof of Lemma 2. We restate it here for convenience of the reader.

Lemma** (Lemma 4.5).**

Let Assumptions 1 and 2 hold. Let $\{x_{k}^{h}\}$ be the sequence generated by Algorithm 2 and $x_{k}^{*}$ be a projection point of $x_{k}^{h}$ onto $\mathcal{X}$. Assume that it exists a strictly positive constant $\underline{\rho}$ such that $\{x_{k}^{h}\}\in\mathcal{B}(x_{k}^{*},\underline{\rho})$ and that $\nabla_{x}^{2}f$ is Lipschitz continuous in $\mathcal{B}(x_{k}^{*},\underline{\rho})$ with Lipschitz constant $L_{2}$. Then, it holds:

[TABLE]

with

[TABLE]

with $\lambda_{\max}$ and $\kappa_{R}$ defined respectively in (23) and Assumption 1.

Proof.

The proof is divided into two parts. We first consider the case in which $s_{k}^{h}$ has been obtained from the approximate minimization of the Taylor model, and then the case in which it has been obtained as prolongation of the step obtained from the approximate minimization of the coarse model.

Le us then assume that the Taylor model has been employed. Reminding that $\nabla_{x}f^{h}(x^{*})=0$ for each $x^{*}\in\mathcal{X}$, and definition (2) we obtain:

[TABLE]

with

[TABLE]

Some algebraic manipulations (adding $\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k+1}^{h}-x_{k}^{*})$ to both sides of (34) and expressing $(x_{k+1}^{h}-x_{k}^{*})=s_{k}^{h}+(x_{k}^{h}-x_{k}^{*})$) lead to:

[TABLE]

Using the fact that $\nabla_{x}^{2}f^{h}(x_{k}^{*})\succeq 0$, the stopping criterion (18), and the triangle inequality, it follows

[TABLE]

Using the Lipschitz continuity of $\nabla_{x}^{2}f$ in $\mathcal{B}(x_{k}^{*},\underline{\rho})$, the relation (16) with $q=2$ and the triangle inequality $\|x_{k+1}^{h}-x_{k}^{*}\|\geq\|x_{k+1}^{h}-x_{k}^{h}\|-\|x_{k}^{h}-x_{k}^{*}\|=\|s_{k}^{h}\|-\|x_{k}^{h}-x_{k}^{*}\|$, we obtain:

[TABLE]

Notice that

[TABLE]

We can then study when the inequality holds

[TABLE]

The right hand side of the inequality is expressed as a polynomial of $\|s_{k}^{h}\|$ of order $1$ with positive value in 0, so the inequality will be true if $\|s_{k}^{h}\|$ is small enough. We can then assume $\|s_{k}^{h}\|<1$, so that $\|s_{k}^{h}\|^{q}\leq\|s_{k}^{h}\|^{2}$ if $q\geq 2$. Then, we have that

[TABLE]

We can then solve

[TABLE]

The solution leads to

[TABLE]

Let us now consider the case in which the lower level model is used. The idea is similar as in the previous case. Reminding (25) and that $R\nabla_{x}f^{h}(x_{k}^{*})=0$, we have:

[TABLE]

Algebraic manipulations (adding $R\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k+1}^{h}-x_{k}^{*})$ to both sides and expressing $(x_{k+1}^{h}-x_{k}^{*})=s_{k}^{h}+(x_{k}^{h}-x_{k}^{*})$) lead to:

[TABLE]

where

[TABLE]

Further, we can write $R\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k+1}^{h}-x_{k}^{*})=R\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k+1}^{h}-x_{k}^{h})+R\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k}^{h}-x_{k}^{*})=R\nabla_{x}^{2}f^{h}(x_{k}^{*})Ps_{k}^{H}+R\nabla_{x}^{2}f^{h}(x_{k}^{*})(x_{k}^{h}-x_{k}^{*})$:

[TABLE]

We can again use relation (16) (applied to $f^{H},T_{q,k}^{H}$ with constant $L_{H}$ and to $f,T_{2,k}^{h}$ with constant $L_{2}$), (18), the fact that $R\nabla_{x}^{2}f^{h}(x_{k}^{*})P$ is still positive definite, and Assumption 1 together with relation $s_{k}^{h}=Ps_{k}^{H}$, to deduce that:

[TABLE]

We remark that

[TABLE]

As previously, we can solve the following inequality:

[TABLE]

and conclude that:

[TABLE]

We can then use the fact that $\lambda_{k}\leq\lambda_{\max}$ for all $k$ and that $\|s_{k}^{h}\|\leq\kappa_{R}\|s_{k}^{H}\|$ to conclude that in all cases it exists a constant $C$ such that $\|s_{k}^{h}\|\leq C\|x_{k}^{h}-x_{k}^{*}\|$. ∎