Model Construction for Convex-Constrained Derivative-Free Optimization

Lindon Roberts School of Mathematics and Statistics, University of Sydney, Camperdown NSW 2006, Australia ([email protected]). This work was supported by the Australian Research Council Discovery Early Career Award DE240100006.

Abstract

We develop a new approximation theory for linear and quadratic interpolation models, suitable for use in convex-constrained derivative-free optimization (DFO). Most existing model-based DFO methods for constrained problems assume the ability to construct sufficiently accurate approximations via interpolation, but the standard notions of accuracy (designed for unconstrained problems) may not be achievable by only sampling feasible points, and so may not give practical algorithms. This work extends the theory of convex-constrained linear interpolation developed in [Hough & Roberts, SIAM J. Optim, 32:4 (2022), pp. 2552–2579] to the case of linear regression models and underdetermined quadratic interpolation models.

Keywords: derivative-free optimization, convex constraints, interpolation

Mathematics Subject Classification: 65K05, 90C30, 90C56

1 Introduction

This work is concerned with model-based methods for derivative-free optimization (DFO). Such methods are widely studied and have proven to be popular in practice [16, 1, 25]. Most commonly, model-based methods construct a local interpolant for the objective from selected evaluations which is then used in place of Taylor series in a standard framework, most commonly trust-region methods [13]. The convergence and worst-case complexity of model-based DFO in the unconstrained setting is well-established. Here, we are interested in nonconvex minimization over a convex set,

\displaystyle\min_{\bm{x}\in C}f(\bm{x}),

(1.1)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is smooth (but nonconvex), and $C\subseteq\mathbb{R}^{n}$ is a closed, convex set with nonempty interior. Such problems appear in the derivative-free setting, for example bound and linear constraints appearing in problems from climate science [31, 30], deep learning [33] and systems design engineering [23]. Although a model-based DFO framework for solving (1.1) was recently developed in [22], it only considered the case of constructing linear interpolants from function evaluations (at feasible points). Outside of structured problems such as nonlinear least-squares [8], quadratic interpolants are the most common models used in model-based DFO [28, 29, 16]. In this work, we develop a comprehensive approximation theory for linear regression and underdetermined quadratic interpolation models based on sampling only from feasible points. This theory fits the requirements of the method from [22], and so allows this framework to be used with more practical model choices.

1.1 Existing work

Model-based DFO methods for unconstrained problems have well-established convergence [13, 16] and worst-case complexity theory [17, 25]. This theory relies on building (typically linear or quadratic) interpolants for the objective function based on samples of the objective at specific points, which must be “fully linear”, i.e. satisfying specific error bounds. Unlike much approximation theory in numerical analysis (e.g. [32]), it is typically assumed that we have limited ability to select the sample points and instead should use as many pre-existing evaluations as possible (motivated by problems with computationally expensive objectives, an important DFO use case). In this setting, the model accuracy can be bounded by quantities only depending on the choice of sampling points (and information about the smoothness of the objective). Originally in [13], the concept of “ $\Lambda$ -poised” interpolation sets was introduced as a sufficient condition for constructing fully linear interpolation models, as well as developing procedures for improving the choice of sample points (while retaining as many existing points as possible) to construct a $\Lambda$ -poised set. This was extended in [14] for polynomial regression and underdetermined polynomial interpolation models. A specific form of underdetermined quadratic interpolation, namely minimum Frobenius interpolation, was introduced and extensively studied in [26, 27] and is the foundation of several state-of-the-art codes [28, 29, 6]. An analysis of minimum Frobenius interpolation in terms of $\Lambda$ -poisedness/fully linear error bounds was given in [16]. Alternate model constructions have been proposed based on weighted linear regression [5], underdetermined quadratic interpolation based on Sobolev norms [37], radial basis functions (RBF) [35] and combinations of polynomial and RBFs [2], for example.

Many model-based DFO methods have been proposed for constrained optimization; see [25, Section 7] for a thorough survey. We note in particular the following strictly feasible methods: [10] studies convex-constrained problems, [29, 19] and [34, Section 6.3] consider bound constraints, and [20] considers linear inequality constraints. Strictly feasible methods for nonconvex-constrainted problems have also been considered in [11, 9, 36]. Of these, where convergence was established it was based on the assumption that fully linear models could be constructed using the same $\Lambda$ -poised interpolation sets as in the unconstrained case. However, this is potentially impractical, as in some settings “it may be impossible to obtain a fully linear model using only feasible points” [25, p. 362].

More recently, [22] introduces a model-based DFO method for solving (1.1) by adapting the unconstrained approach from [15] based on gradient-based trust-region methods for convex-constrained optimization [12, Chapter 12]. The method has convergence and worst-case complexity results which align with unconstrained model-based DFO [17] and convex-constrained gradient-based trust-region methods [7]. To make this approach practical, it introduced a weaker notion of fully linear models which are sufficient for algorithm convergence, but which can be attained using linear interpolation over a $\Lambda$ -poised set of only feasible points. It also demonstrated how to modify an interpolation set to make it $\Lambda$ -poised (without changing too many points). However, linear models are not typically used in practice, except for special cases (such as nonlinear least-squares minimization [8]), with quadratic models being preferred [28, 29, 6].

We additionally note that the BOBYQA code [29] for bound-constrained optimization uses maximization of Lagrange polynomials inside the feasible region—similar to how we will define $\Lambda$ -poisedness in this setting—but this was a heuristic with no theoretical analysis.¹¹1In fact, our results here give some theoretical justification to this approach.

1.2 Contributions

The contributions of this work are:

•

A generalization of the linear interpolation theory (based on sampling feasible points) from [22] to the case of linear regression models. This includes fully linear error bounds and procedures to generate $\Lambda$ -poised interpolation sets (which are simpler than those developed in [22]).
•

A full approximation theory for underdetermined quadratic interpolation (in the style of [26, 27]) based on sampling only feasible points. This includes derivation of fully linear error bounds, defining the relevant notion of $\Lambda$ -poisedness in this setting, and procedures for efficiently generating $\Lambda$ -poised sets. By construction, in the case of a full interpolation set of $(n+1)(n+2)/2$ points for quadratic interpolation, this approach produces the same result as standard quadratic interpolation, and so our results also apply to this setting.

While the linear regression theory is of independent interest for DFO applied to noisy objectives [24], it is primarily included as a prerequisite for deriving the quadratic results. The results proven in both cases are sufficient to enable these model constructions to be used in the algorithmic framework from [22] (with the associated convergence and worst-case complexity guarantees).

While the formulations and results are very similar to those for the unconstrained case from [14], the constrained setting necessitates new approaches for proving the results (which when applied to the $C=\mathbb{R}^{n}$ case provide simpler proofs of existing results).

Organization

A summary of the modified fully linear model accuracy requirement and associated algorithm from [22] is given in Section 2. That $\Lambda$ -poisedness (only over the feasible region) guarantees fully linear interpolation models is first shown for linear regression models in Section 3 and then for (underdetermined) quadratic interpolation models in Section 4. Lastly, Section 5 provides concrete procedures for identifying and constructing $\Lambda$ -poised quadratic interpolation models using only feasible points.

Notation

We let $\|\cdot\|$ be the Euclidean norm of vectors and the operator 2-norm of matrices, and use $B(\bm{x},\Delta)$ for $\bm{x}\in\mathbb{R}^{n}$ and $\Delta>0$ for the closed ball $\{\bm{y}\in\mathbb{R}^{n}:\|\bm{y}-\bm{x}\|\leq\Delta\}$ .

2 Background

This section summarizes the key background and results from [22], which provide a model-based DFO algorithm for solving (1.1).

We now outline how problem (1.1) can be solved using model-based DFO methods. We begin with outlining our problem assumptions.

Assumption 2.1.

The objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ is bounded below by $f_{\textnormal{low}}$ and continuously differentiable, and $\nabla f$ is $L_{\nabla f}$ -Lipschitz continuous.

Although $\nabla f$ exists, we assume that it is not available to the algorithm as an oracle.

Assumption 2.2.

The feasible set $C\subseteq\mathbb{R}^{n}$ is closed and convex with nonempty interior.

Assumption 2.2 implies that the Euclidean projection operator for $C$ ,

\displaystyle\operatorname{proj}_{C}(\bm{y}):=\operatorname*{arg\,min}_{\bm{x}\in C}\|\bm{y}-\bm{x}\|^{2},

(2.1)

is a well-defined function $\operatorname{proj}_{C}:\mathbb{R}^{n}\to C$ [3, Theorem 6.25]. Our framework is based on a setting in which $\operatorname{proj}_{C}$ is inexpensive to evaluate (at least compared to evaluations of the objective $f$ ). Several examples of feasible sets $C$ which have simple analytic expressions for $\operatorname{proj}_{C}$ are given in [3, Table 6.1].

In this work, as in [22], we aim to find a first-order optimal solution to (1.1). To measure this, we use the following first-order criticality measure from [12, Section 12.1.4], defined for all $\bm{x}\in C$ :

\displaystyle\pi^{f}(\bm{x}):=\left|\min_{\begin{subarray}{c}\bm{x}+\bm{d}\in C\\ \|\bm{d}\|\leq 1\end{subarray}}\nabla f(\bm{x})^{T}\bm{d}\right|.

(2.2)

From [12, Theorem 12.1.6], we have that $\pi^{f}$ is a non-negative, continuous function of $\bm{x}$ and $\pi^{f}(\bm{x}^{*})=0$ if and only if $\bm{x}^{*}$ is a first-order critical point for (1.1). In the unconstrained case $C=\mathbb{R}^{n}$ , then (2.2) simplifies to $\pi^{f}(\bm{x})=\|\nabla f(\bm{x})\|$ . For notational convenience, we define $\pi^{f}_{k}:=\pi^{f}(\bm{x}_{k})$ , where $\bm{x}_{k}\in\mathbb{R}^{n}$ is the $k$ -th iterate of our algorithm.

2.1 Model Construction & Accuracy

We will consider a model-based DFO approach for solving (1.1). Specifically, at each iteration $k$ we construct a local quadratic approximation to $f$ which we hope to be accurate near the current iterate $\bm{x}_{k}\in\mathbb{R}^{n}$ :

\displaystyle f(\bm{y})\approx m_{k}(\bm{y}):=c_{k}+\bm{g}_{k}^{T}(\bm{y}-\bm{x}_{k})+\frac{1}{2}(\bm{y}-\bm{x}_{k})^{T}H_{k}(\bm{y}-\bm{x}_{k}),

(2.3)

for some $c_{k}\in\mathbb{R}$ , $\bm{g}_{k}\in\mathbb{R}^{n}$ and $H_{k}\in\mathbb{R}^{n\times n}$ with $H_{k}=H_{k}^{T}$ .

Assumption 2.3.

There exists $\kappa_{H}\geq 1$ such that $\|H_{k}\|\leq\kappa_{H}-1$ for all $k$ .

Convergence of our algorithm will require that the local models $m_{k}$ (2.3) are sufficiently accurate approximations of the objective $f$ (at least on some iterations). The required notion of ‘sufficiently accurate’ is the following:

Definition 2.4.

The model $m_{k}$ (2.3) is $C$ -fully linear in $B(\bm{x}_{k},\Delta_{k})$ if there exist $\kappa_{\textnormal{ef}},\kappa_{\textnormal{eg}}>0$ , independent of $k$ , such that


$\displaystyle\max_{\begin{subarray}{c}\bm{x}_{k}+\bm{d}\in C\\ \\|\bm{d}\\|\leq\Delta_{k}\end{subarray}}\|f(\bm{x}_{k}+\bm{d})-m_{k}(\bm{x}_{k}+\bm{d})\|$	$\displaystyle\leq\kappa_{\textnormal{ef}}\Delta_{k}^{2},$	(2.4a)
$\displaystyle\max_{\begin{subarray}{c}\bm{x}_{k}+\bm{d}\in C\\ \\|\bm{d}\\|\leq 1\end{subarray}}\left\|(\nabla f(\bm{x}_{k})-\bm{g}_{k})^{T}\bm{d}\right\|$	$\displaystyle\leq\kappa_{\textnormal{eg}}\Delta_{k}.$	(2.4b)

When we wish to refer to this notion in the abstract sense, we use the term “ $C$ -full linearity”.

We note in particular that (2.4a) uses the constraint $\|\bm{d}\|\leq\Delta_{k}$ and (2.4b) uses $\|\bm{d}\|\leq 1$ . This essentially corresponds to the standard definition of fully linear models [15, Definition 3.1] in the unconstrained case $C=\mathbb{R}^{n}$ .

Our algorithm will require the existence of two procedures related to $C$ -full linearity: given a model $m_{k}$ , iterate $\bm{x}_{k}$ and trust-region radius $\Delta_{k}$ , we assume that we can

•

Verify whether or not $m_{k}$ is $C$ -fully linear in $B(\bm{x}_{k},\Delta_{k})$ ; and
•

If $m_{k}$ is not $C$ -fully linear, create a new model (typically a small modification of $m_{k}$ ) which is $C$ -fully linear in $B(\bm{x}_{k},\Delta_{k})$ .

Lastly, our model $m_{k}$ (2.3) induces an approximate first-order criticality measure, namely

\displaystyle\pi^{m}(\bm{x}):=\left|\min_{\begin{subarray}{c}\bm{x}+\bm{d}\in C\\ \|\bm{d}\|\leq 1\end{subarray}}\bm{g}_{k}^{T}\bm{d}\right|.

(2.5)

As above, for convenience we will define $\pi^{m}_{k}:=\pi^{m}(\bm{x}_{k})$ .

2.2 Algorithm Specification

Our DFO method is based on a trust-region framework. That is, given the local model $m_{k}$ (2.3), we calculate a potential new iterate as $\bm{x}_{k}+\bm{s}_{k}$ , where $\bm{s}_{k}$ is an approximate minimizer of the trust-region subproblem

\displaystyle\min_{\begin{subarray}{c}\bm{x}+\bm{s}\in C\\ \|\bm{s}\|\leq\Delta_{k}\end{subarray}}m_{k}(\bm{x}_{k}+\bm{s}),

(2.6)

where $\Delta_{k}>0$ is a parameter updated dynamically inside the algorithm. Formally, we require the following, which can be achieved using a Goldstein-type linesearch method [12, Algorithm 12.2.2 & Theorem 12.2.2].

Assumption 2.5.

There exists a constant $c_{1}\in(0,1)$ such that the computed step $\bm{s}_{k}$ satisfies $\bm{x}_{k}+\bm{s}_{k}\in C$ , $\|\bm{s}_{k}\|\leq\Delta_{k}$ and the generalized Cauchy decrease condition:

\displaystyle m_{k}(\bm{x}_{k})-m_{k}(\bm{x}_{k}+\bm{s}_{k})\geq c_{1}\pi^{m}_{k}\min\left(\frac{\pi^{m}_{k}}{1+\|H_{k}\|},\Delta_{k},1\right).

(2.7)

The full derivative-free algorithm for solving (1.1) from [22], called CDFO-TR, is given in Algorithm 1. The overall structure is similar to other model-based DFO methods, such as [15, Algorithm 4.1] for unconstrained minimization.

1:Starting point

\bm{x}_{0}\in\mathbb{R}^{n}

and trust-region radius

\Delta_{0}>0

2:Parameters: maximum trust-region radius

\Delta_{\max}\geq\Delta_{0}

, scaling factors

0<\gamma_{\textnormal{dec}}<1<\gamma_{\textnormal{inc}}

, criticality constants

\epsilon_{C},\mu>0

, and acceptance threshold

\eta\in(0,1)

3:Build an initial model

m_{0}

(2.3).

4:for

k=0,1,2,\ldots

5: if

\pi^{m}_{k}<\epsilon_{C}

and (

\pi^{m}_{k}<\mu^{-1}\Delta_{k}

m_{k}

is not

C

-fully linear in

B(\bm{x}_{k},\Delta_{k})

) then

6: Criticality step: Set

\bm{x}_{k+1}=\bm{x}_{k}

. If

m_{k}

C

-fully linear in

B(\bm{x}_{k},\Delta_{k})

, set

\Delta_{k+1}=\gamma_{\textnormal{dec}}\Delta_{k}

, otherwise set

\Delta_{k+1}=\Delta_{k}

. Construct

m_{k+1}

to be

C

-fully linear in

B(\bm{x}_{k+1},\Delta_{k+1})

7: else

\leftarrow

\pi^{m}_{k}\geq\epsilon_{C}

or (

\pi^{m}_{k}\geq\mu^{-1}\Delta_{k}

and

m_{k}

C

-fully linear in

B(\bm{x}_{k},\Delta_{k})

)

8: Approximately solve (2.6) to get a step

\bm{s}_{k}

9: Evaluate

f(\bm{x}_{k}+\bm{s}_{k})

and calculate ratio

\displaystyle\rho_{k}:=\frac{f(\bm{x}_{k})-f(\bm{x}_{k}+\bm{s}_{k})}{m_{k}(\bm{x}_{k})-m_{k}(\bm{x}_{k}+\bm{s}_{k})}.

(2.8)

10: if

\rho_{k}\geq\eta

then

11: Successful step: Set

\bm{x}_{k+1}=\bm{x}_{k}+\bm{s}_{k}

and

\Delta_{k+1}=\min(\gamma_{\textnormal{inc}}\Delta_{k},\Delta_{\max})

. Form

m_{k+1}

in any manner.

12: else if

m_{k}

is not

C

-fully linear in

B(\bm{x}_{k},\Delta_{k}

) then

13: Model-improving step: Set

\bm{x}_{k+1}=\bm{x}_{k}

and

\Delta_{k+1}=\Delta_{k}

, and construct

m_{k+1}

to be

C

-fully linear in

B(\bm{x}_{k+1},\Delta_{k+1})

14: else

\leftarrow

\rho_{k}<\eta

and

m_{k}

C

-fully linear in

B(\bm{x}_{k},\Delta_{k})

15: Unsuccessful step: Set

\bm{x}_{k+1}=\bm{x}_{k}

and

\Delta_{k+1}=\gamma_{\textnormal{dec}}\Delta_{k}

. Form

m_{k+1}

in any manner.

16: end if

17: end if

18:end for

Algorithm 1 CDFO-TR: model-based DFO method for (1.1), from [22].

We have the following global convergence and worst-case complexity results for Algorithm 1 from [22].

Theorem 2.6 (Theorem 3.10 & Corollary 3.15, [22]).

If Assumptions 2.1, 2.3 and 2.5 hold, then $\lim_{k\to\infty}\pi^{f}_{k}=0$ . Moreover, if $\epsilon\in(0,1]$ and $\epsilon_{C}\geq c_{2}\epsilon$ for some constant $c_{2}>0$ , then the number of iterations $k$ before Algorithm 1 produces an iterate with $\pi^{f}_{k}<\epsilon$ is at most $\mathcal{O}(\kappa_{H}\kappa_{d}^{2}\epsilon^{-2})$ , where $\kappa_{d}:=\max(\kappa_{\textnormal{ef}},\kappa_{\textnormal{eg}})$ .

In [22], details are given on how to construct fully linear models based on linear interpolation to feasible points. Although linear models can be practical for some structured problems, such as nonlinear least-squares objectives (e.g. [8]), quadratic models are generally preferred. The remainder of this paper is devoted to the construction of fully linear quadratic models by only sampling the objective at feasible points.

3 Linear Regression Models

We first consider how to extend the linear interpolation approximation theory from [22] to the case of linear regression models (with the slight extra generalization that the base point $\bm{x}$ does not need to be an interpolation point). The purpose of this is twofold: regression models can be useful for modelling noisy objectives (e.g. [5]), and this theory will be necessary to develop the corresponding quadratic interpolation theory in Section 4. In the unconstrained case $C=\mathbb{R}^{n}$ , these results were originally proven in [14].

Suppose we have sampled the function $f$ at $p$ points $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset\mathbb{R}^{n}$ (where $p\geq n+1$ ), and given this information we wish to find a linear model

\displaystyle f(\bm{y})\approx m(\bm{y}):=c+\bm{g}^{T}(\bm{y}-\bm{x}),

(3.1)

by solving

\displaystyle\min_{c,\bm{g}\in\mathbb{R}\times\mathbb{R}^{n}}\>\sum_{t=1}^{p}(f(\bm{y}_{t})-m(\bm{y}_{t}))^{2}.

(3.2)

Equivalently, we can find the $c$ and $\bm{g}$ for our model by finding the least-squares solution to the $p\times(n+1)$ system

\displaystyle M\begin{bmatrix}c\\ \bm{g}\end{bmatrix}:=\begin{bmatrix}1&(\bm{y}_{1}-\bm{x})^{T}\\ \vdots&\vdots\\ 1&(\bm{y}_{p}-\bm{x})^{T}\end{bmatrix}\begin{bmatrix}c\\ \bm{g}\end{bmatrix}=\begin{bmatrix}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\end{bmatrix},

(3.3)

which may be written as

\displaystyle\begin{bmatrix}c\\ \bm{g}\end{bmatrix}=M^{\dagger}\begin{bmatrix}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\end{bmatrix},

(3.4)

where $M^{\dagger}\in\mathbb{R}^{(n+1)\times p}$ is the Moore-Penrose pseudoinverse of $M$ .

Later we will require the following standard properties of $M^{\dagger}$ (see, e.g. [18, Section 5.5.4]).

Lemma 3.1.

For $p\geq n+1$ , the Moore-Penrose pesudoinverse $M^{\dagger}$ of $M$ (3.3) satisfies $(M^{T})^{\dagger}=(M^{\dagger})^{T}$ and $M^{\dagger}MM^{\dagger}=M^{\dagger}$ . The minimal-norm solution to the underdetermined system $M^{T}\bm{u}=\bm{v}$ is $\bm{u}=(M^{T})^{\dagger}\bm{v}$ .

The quality of the choice of interpolation points will be assessed by considering the associated set of Lagrange polynomials. In this case, the Lagrange polynomials associated with our interpolation set are the linear functions

\displaystyle\ell_{t}(\bm{y}):=c_{t}+\bm{g}_{t}^{T}(\bm{y}-\bm{x}),\qquad\forall t=1,\ldots,p,

(3.5)

each defined by the least-squares regression problem

\displaystyle\min_{c_{t},\bm{g}_{t}\in\mathbb{R}\times\mathbb{R}^{n}}\>\sum_{s=1}^{p+1}(\delta_{s,t}-\ell_{t}(\bm{y}_{s}))^{2},\qquad\forall t=1,\ldots,p,

(3.6)

or equivalently

\displaystyle\begin{bmatrix}c_{t}\\ \bm{g}_{t}\end{bmatrix}=M^{\dagger}\bm{e}_{t},\qquad\forall t=1,\ldots,p.

(3.7)

Our notion of the sampled points $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ being a ‘good’ choice is given by the Lagrange polynomials having small magnitude in the region of interest. This is formalized in the following notion, which generalizes [14, Definition 2.7] to the convex feasible region $C$ .

Definition 3.2.

Given $\Lambda\geq 1$ , the set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset C$ is $\Lambda$ -poised for linear regression in $B(\bm{x},\Delta)\cap C$ if $\{\bm{y}_{2}-\bm{y}_{1},\ldots,\bm{y}_{p}-\bm{y}_{1}\}$ spans $\mathbb{R}^{n}$ and

\displaystyle\max_{t=1,\ldots,p}|\ell_{t}(\bm{y})|\leq\Lambda,\qquad\forall\bm{y}\in C\cap B(\bm{x},\min(\Delta,1)).

(3.8)

We note that if $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ is $\Lambda$ -poised for linear regression, then $M$ has full column rank (since $\{\bm{y}_{2}-\bm{y}_{1},\ldots,\bm{y}_{p}-\bm{y}_{1}\}$ is assumed to span $\mathbb{R}^{n}$ ). This requirement also necessitates that $p\geq n+1$ .

Assuming the $\Lambda$ -poisedness of $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ will be sufficient for the associated linear regression model to be $C$ -fully linear. The proof of this follows a similar structure to [22, Lemma 4.3 & Theorem 4.4], but with an increased complexity to the proofs coming from using regression instead of interpolation.

Lemma 3.3.

Suppose $f$ satisfies Assumption 2.1 and $C$ satisfies Assumption 2.2. Then if $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ is $\Lambda$ -poised for linear regression in $B(\bm{x},\Delta)\cap C$ and $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ for all $t=1,\ldots,p$ and some $\beta>0$ , we have

\displaystyle|m(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})|\leq\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2},

(3.9)

for all $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ . If we also have $\bm{x}\in C$ , this in turn implies error bounds on $c$ and $\bm{g}$ individually, namely

\displaystyle|c-f(\bm{x})|\leq\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2},

(3.10)

and

\displaystyle|(\bm{y}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))|\leq p\Lambda L_{\nabla f}\beta^{2}\min(\Delta,1)^{2},

(3.11)

for all $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ .

Proof.

We begin by defining the residual of the least-squares problem (3.3) as

\displaystyle\bm{r}:=M\begin{bmatrix}c\\ \bm{g}\end{bmatrix}-\begin{bmatrix}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\end{bmatrix}=(MM^{\dagger}-I)\begin{bmatrix}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\end{bmatrix}\in\mathbb{R}^{p}.

(3.12)

Then for all $t=1,\ldots,p$ have

	$\displaystyle m(\bm{y}_{t})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x})$	$\displaystyle=\begin{bmatrix}c\\ \bm{g}\end{bmatrix}^{T}\underbrace{\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix}}_{=M^{T}\bm{e}_{t}}-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x}),$		(3.13)
		$\displaystyle=\bm{r}^{T}\bm{e}_{t}+f(\bm{y}_{t})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x}).$		(3.14)

Now, fix $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ . Since $M$ has full column rank, its rows span $\mathbb{R}^{n+1}$ and so there exist constants $\alpha_{t}(\bm{y})$ such that

\displaystyle\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}=\sum_{t=1}^{p}\alpha_{t}(\bm{y})\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix}=M^{T}\bm{\alpha}(\bm{y}).

(3.15)

Of these, we take the $\bm{\alpha}(\bm{y})$ with minimal norm (c.f. Lemma 3.1),

\displaystyle\bm{\alpha}(\bm{y})=(M^{T})^{\dagger}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}.

(3.16)

There are two key properties of $\bm{\alpha}(\bm{y})$ : firstly,

\displaystyle\ell_{t}(\bm{y})=\begin{bmatrix}c_{t}\\ \bm{g}_{t}\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}=\bm{e}_{t}^{T}(M^{\dagger})^{T}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}=\alpha_{t}(\bm{y}),

(3.17)

where the last equality follows $(M^{T})^{\dagger}=(M^{\dagger})^{T}$ (Lemma 3.1). Hence $|\alpha_{t}(\bm{y})|=|\ell_{t}(\bm{y})|\leq\Lambda$ from the $\Lambda$ -poisedness condition. Secondly, we have

\displaystyle\bm{\alpha}(\bm{y})^{T}\bm{r}=\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}^{T}M^{\dagger}(MM^{\dagger}-I)\begin{bmatrix}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\end{bmatrix}=0,

(3.18)

from $M^{\dagger}MM^{\dagger}=M^{\dagger}$ (Lemma 3.1). All together, we have

$\displaystyle\|m(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|$	$\displaystyle=\left\|\begin{bmatrix}c\\ \bm{g}\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}-\begin{bmatrix}f(\bm{x})\\ \nabla f(\bm{x})\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}\right\|,$	(3.19)
	$\displaystyle=\left\|\sum_{t=1}^{p}\alpha_{t}(\bm{y})\left(\begin{bmatrix}c\\ \bm{g}\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix}-\begin{bmatrix}f(\bm{x})\\ \nabla f(\bm{x})\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix}\right)\right\|,$	(3.20)
	$\displaystyle=\left\|\sum_{t=1}^{p}\alpha_{t}(\bm{y})\left\{\bm{r}^{T}\bm{e}_{t}+f(\bm{y}_{t})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x})\right\}\right\|,$	(3.21)
	$\displaystyle\leq\left\|\bm{\alpha}(\bm{y})^{T}\bm{r}\right\|+\frac{L_{\nabla f}}{2}\sum_{t=1}^{p}\|\alpha_{t}(\bm{y})\|\cdot\\|\bm{y}_{t}-\bm{x}\\|^{2},$	(3.22)
	$\displaystyle\leq\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2},$	(3.23)

and we recover (3.9), where we used (3.18), $|\alpha_{t}(\bm{y})|\leq\Lambda$ and $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ to get the last inequality.

If $\bm{x}\in C$ , we can take $\bm{y}=\bm{x}$ in (3.9) to get (3.10). Combining (3.10) with (3.9) we get

	$\displaystyle\|(\bm{y}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle\leq\|c+\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|+\|c-f(\bm{x})\|,$		(3.24)
		$\displaystyle\leq p\Lambda L_{\nabla f}\beta^{2}\min(\Delta,1)^{2},$		(3.25)

and we get (3.11). ∎

Theorem 3.4.

Suppose the assumptions of Lemma 3.3 hold and $\bm{x}\in C$ . Then the regression model (3.1) is $C$ -fully linear in $B(\bm{x},\Delta)$ with constants

\displaystyle\kappa_{\textnormal{ef}}=p\Lambda L_{\nabla f}\beta^{2}+\frac{L_{\nabla f}}{2},\qquad\text{and}\qquad\kappa_{\textnormal{eg}}=p\Lambda L_{\nabla f}\beta^{2},

(3.26)

in (2.4).

Proof.

Fix $\bm{y}\in B(\bm{x},\Delta)\cap C$ . We first derive $\kappa_{\textnormal{ef}}$ by considering the cases $\Delta>1$ and $\Delta\leq 1$ separately.

If $\Delta>1$ , then $\hat{\bm{y}}=\bm{x}+\frac{1}{\Delta}(\bm{y}-\bm{x})\in B(\bm{x},1)=B(\bm{x},\min(\Delta,1))$ . Since $C$ is convex and $\bm{x},\bm{y}\in C$ we also have $\hat{\bm{y}}\in C$ . Hence (3.9) gives

\displaystyle|c+\bm{g}^{T}(\hat{\bm{y}}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\hat{\bm{y}}-\bm{x})|

\displaystyle\leq\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2},

(3.27)

and so

\displaystyle\left|c+\frac{1}{\Delta}\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\frac{1}{\Delta}\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right|

\displaystyle\leq\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2}.

(3.28)

This gives us

$\displaystyle\|f(\bm{y})-m(\bm{y})\|$	$\displaystyle\leq\|f(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|+\|m(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|,$	(3.29)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\\|\bm{y}-\bm{x}\\|^{2}+\Delta\left\|\frac{1}{\Delta}c+\frac{1}{\Delta}\bm{g}^{T}(\bm{y}-\bm{x})-\frac{1}{\Delta}f(\bm{x})-\frac{1}{\Delta}\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right\|,$	(3.30)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\\|\bm{y}-\bm{x}\\|^{2}+\Delta\left[\left\|c+\frac{1}{\Delta}\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\frac{1}{\Delta}\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right\|\right.$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\left.+\left\|\left(\frac{1}{\Delta}-1\right)(c-f(\bm{x}))\right\|\right],$	(3.31)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\Delta^{2}+\Delta\left[\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2}\right.$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\left.+\left(1-\frac{1}{\Delta}\right)\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2}\right],$	(3.32)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\Delta^{2}+p\Lambda L_{\nabla f}\beta^{2}\Delta\min(\Delta,1)^{2},$	(3.33)
	$\displaystyle=\frac{L_{\nabla f}}{2}\Delta^{2}+p\Lambda L_{\nabla f}\beta^{2}\Delta^{2},$	(3.34)

where we use (3.10) to get the fourth inequality, and the last line follows from $\Delta>1$ . Instead if $\Delta\leq 1$ , then $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ already, and so (3.9) immediately gives

$\displaystyle\|f(\bm{y})-m(\bm{y})\|$	$\displaystyle=\|f(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|+\left\|m(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right\|,$	(3.35)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\Delta^{2}+\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\min(\Delta,1)^{2},$	(3.36)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\Delta^{2}+\frac{p\Lambda L_{\nabla f}\beta^{2}}{2}\Delta^{2}.$	(3.37)

Either way, we get the desired value of $\kappa_{\textnormal{ef}}$ .

To get $\kappa_{\textnormal{eg}}$ we now fix an arbitrary $\widetilde{\bm{y}}\in B(\bm{x},1)\cap C$ and again consider the cases $\Delta\geq 1$ and $\Delta<1$ separately. First, if $\Delta\geq 1$ , then $\widetilde{\bm{y}}\in B(\bm{x},\min(\Delta,1))\cap C$ . From (3.11) we get

\displaystyle|(\widetilde{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))|\leq p\Lambda L_{\nabla f}\beta^{2}\min(\Delta,1)^{2}\leq p\Lambda L_{\nabla f}\beta^{2}\Delta,

(3.38)

where the second inequality follows from $\Delta\geq 1$ . Alternatively, if $\Delta<1$ then the convexity of $C$ implies that $\hat{\bm{y}}:=\bm{x}+\Delta(\widetilde{\bm{y}}-\bm{x})\in B(\bm{x},\Delta)\cap C=B(\bm{x},\min(\Delta,1))\cap C$ . Again from (3.11) we get

$\displaystyle\|(\widetilde{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle=\frac{1}{\Delta}\|(\hat{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|,$	(3.39)
	$\displaystyle\leq p\Lambda L_{\nabla f}\beta^{2}\Delta^{-1}\min(\Delta,1)^{2},$	(3.40)
	$\displaystyle=p\Lambda L_{\nabla f}\beta^{2}\Delta,$	(3.41)

where the last line follows from $\Delta<1$ . Again, either way we get the desired value of $\kappa_{\textnormal{eg}}$ . ∎

4 Underdetermined Quadratic Interpolation Models

We now consider the case of forming $C$ -fully linear quadratic interpolation models. Our approach follows that of [26] for the unconstrained case, although the fully linear error bounds for this approach were shown later in [16]. We note that this approach is different to the underdetermined quadratic interpolation used in [14].

Here, we aim to construct a quadratic model

\displaystyle f(\bm{y})\approx m(\bm{y}):=c+\bm{g}^{T}(\bm{y}-\bm{x})+\frac{1}{2}(\bm{y}-\bm{x})^{T}H(\bm{y}-\bm{x}),

(4.1)

where $c\in\mathbb{R}$ , $\bm{g}\in\mathbb{R}^{n}$ and $H\in\mathbb{R}^{n\times n}$ with $H=H^{T}$ . We assume that our interpolation set is $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ with $p\in\{n+2,\ldots,(n+1)(n+2)/2\}$ . The case $p=(n+1)(n+2)/2$ corresponds to full quadratic interpolation [13, Section 4.2]. We exclude the case $p=n+1$ which, from (4.2) below, corresponds to linear interpolation and was analyzed (in the convex-constrained case) in [22, Section 4].

If $n+2\leq p<(n+1)(n+2)/2$ then there are infinitely many models satisfying the interpolation conditions $f(\bm{y}_{t})=m(\bm{y}_{t})$ for all $t=1,\ldots,p$ . So, following [26] we choose the model with minimum Frobenius norm Hessian by solving


$\displaystyle\min_{c,\bm{g},H\in\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times n}}$	$\displaystyle\>\>\frac{1}{4}\\|H\\|_{F}^{2},$	(4.2a)
s.t.	$\displaystyle f(\bm{y}_{t})=m(\bm{y}_{t}),\qquad\forall t=1,\ldots,p.$	(4.2b)

This is a convex quadratic program, and (as shown in [26]), reduces to solving the $(p+n+1)\times(p+n+1)$ linear system

\displaystyle F\left[\begin{array}[]{c}\lambda_{1}\\ \vdots\\ \lambda_{p}\\ \hline\cr c\\ \bm{g}\end{array}\right]:=\left[\begin{array}[]{c|c}Q&M\\ \hline\cr M^{T}&0\end{array}\right]\left[\begin{array}[]{c}\lambda_{1}\\ \vdots\\ \lambda_{p}\\ \hline\cr c\\ \bm{g}\end{array}\right]=\left[\begin{array}[]{c}f(\bm{y}_{1})\\ \vdots\\ f(\bm{y}_{p})\\ \hline\cr 0\\ \bm{0}\end{array}\right]\in\mathbb{R}^{p+n+1},

(4.20)

where $M\in\mathbb{R}^{p\times(n+1)}$ is from the linear regression problem (3.3) and $Q\in\mathbb{R}^{p\times p}$ has entries $Q_{i,j}=\frac{1}{2}[(\bm{y}_{i}-\bm{x})^{T}(\bm{y}_{j}-\bm{x})]^{2}$ for $i,j=1,\ldots,p$ . The solution to (4.20) immediately gives us $c$ and $\bm{g}$ ; the (symmetric) model Hessian is given by $H=\sum_{t=1}^{p}\lambda_{t}(\bm{y}_{t}-\bm{x})(\bm{y}_{t}-\bm{x})^{T}$ . We note that $Q$ is symmetric positive semidefinite [26, eq. 2.10], and $\lambda_{1},\ldots,\lambda_{p}$ are the Lagrange multipliers associated with constraints (4.2b).

The Lagrange polynomials associated with our interpolation set are

\displaystyle\ell_{t}(\bm{y}):=c_{t}+\bm{g}_{t}^{T}(\bm{y}-\bm{x})+\frac{1}{2}(\bm{y}-\bm{x})^{T}H_{t}(\bm{y}-\bm{x}),\qquad\forall t=1,\ldots,p,

(4.21)

where $c_{t}$ , $\bm{g}_{t}$ and $H_{t}$ come from solving (4.2) with (4.2b) replaced by $m(\bm{y}_{s})=\delta_{s,t}$ for $s=1,\ldots,p$ . Equivalently, we have

\displaystyle F\begin{bmatrix}\bm{\lambda}_{t}\\ c_{t}\\ \bm{g}_{t}\end{bmatrix}=\bm{e}_{t},\qquad\forall t=1,\ldots,p,

(4.22)

and $H_{t}=\sum_{s=1}^{p}[\bm{\lambda}_{t}]_{s}(\bm{y}_{s}-\bm{x})(\bm{y}_{s}-\bm{x})^{T}$ . This gives

$\displaystyle\ell_{t}(\bm{y})$	$\displaystyle=c_{t}+\bm{g}_{t}(\bm{y}-\bm{x})+\frac{1}{2}\sum_{s=1}^{p}[\bm{\lambda}_{t}]_{s}\left[(\bm{y}-\bm{x})^{T}(\bm{y}_{s}-\bm{x})\right]^{2},$	(4.23)
	$\displaystyle=\begin{bmatrix}\bm{\lambda}_{t}\\ c_{t}\\ \bm{g}_{t}\end{bmatrix}^{T}\underbrace{\begin{bmatrix}\{\frac{1}{2}\left[(\bm{y}-\bm{x})^{T}(\bm{y}_{s}-\bm{x})\right]^{2}\}_{s=1,\ldots,p}\\ 1\\ \bm{y}-\bm{x}\end{bmatrix}}_{=:\bm{\phi}(\bm{y})},$	(4.24)
	$\displaystyle=\bm{e}_{t}^{T}F^{-1}\bm{\phi}(\bm{y}).$	(4.25)

Given (4.20) and (4.22), we conclude that

\displaystyle\begin{bmatrix}\lambda_{1}\\ \vdots\\ \lambda_{p}\end{bmatrix}=\sum_{t=1}^{p}f(\bm{y}_{t})\bm{\lambda}_{t},\qquad c=\sum_{t=1}^{p}f(\bm{y}_{t})c_{t},\quad\text{and}\quad\bm{g}=\sum_{t=1}^{p}f(\bm{y}_{t})\bm{g}_{t},

(4.26)

which gives

\displaystyle m(\bm{y})=\sum_{t=1}^{p}f(\bm{y}_{t})\ell_{t}(\bm{y}).

(4.27)

We now have our new definition of $\Lambda$ -poisedness:

Definition 4.1.

Given $\Lambda\geq 1$ , the interpolation set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset C$ is $\Lambda$ -poised for minimum Frobenius norm interpolation in $B(\bm{x},\Delta)\cap C$ if $F$ is invertible, and

\displaystyle\max_{t=1,\ldots,p}|\ell_{t}(\bm{y})|\leq\Lambda,\qquad\forall\bm{y}\in C\cap B(\bm{x},\min(\Delta,1)).

(4.28)

Lemma 4.2.

If the set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset C$ is $\Lambda$ -poised for minimum Frobenius norm interpolation, then it is $\sqrt{p}\,\Lambda$ -poised for linear regression.

Proof.

Since $F$ is invertible, the sub-matrix $M$ has full column rank by [4, Theorem 3.3], and so $\{\bm{y}_{2}-\bm{y}_{1},\ldots,\bm{y}_{p}-\bm{y}_{1}\}$ spans $\mathbb{R}^{n}$ . The remainder of this proof is based on the argument in [16, p. 83]. We note that since $H=0$ is a global minimizer of the objective function (4.2), that if $f$ is linear then we have exact interpolation, $m=f$ . Applying this to the functions $f(\bm{y})=1$ and $f(\bm{y})=(\bm{y}-\bm{x})^{T}\bm{e}_{i}$ for $i=1,\ldots,n$ , we get from (4.27) that

\displaystyle 1=\sum_{t=1}^{p}\ell_{t}(\bm{y}),\qquad(\bm{y}-\bm{x})^{T}\bm{e}_{i}=\sum_{t=1}^{p}(\bm{y}_{t}-\bm{x})^{T}\bm{e}_{i}\cdot\ell_{t}(\bm{y}),\quad\forall i=1,\ldots,n.

(4.29)

Denoting $\bm{\ell}(\bm{y})\in\mathbb{R}^{p}$ as the vector of all $\ell_{1}(\bm{y}),\ldots,\ell_{p}(\bm{y})$ , these are equivalent to

\displaystyle M^{T}\bm{\ell}(\bm{y})=\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}.

(4.30)

Hence $\bm{\ell}(\bm{y})$ is another solution to the (underdetermined) system (3.15). Since the minimal norm solution to (3.15) was $\bm{\alpha}(\bm{y})$ , we must have $\|\bm{\alpha}(\bm{y})\|\leq\|\bm{\ell}(\bm{y})\|$ . However from (3.17) we know $\bm{\alpha}(\bm{y})=\bm{\ell}^{\text{reg}}(\bm{y})$ , where $\bm{\ell}^{\text{reg}}(\bm{y})$ are the Lagrange polynomials associated with linear regression for $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ . Thus $\|\bm{\ell}^{\text{reg}}(\bm{y})\|\leq\|\bm{\ell}(\bm{y})\|$ . Now, fixing any $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ we get

\displaystyle\|\bm{\ell}^{\text{reg}}(\bm{y})\|_{\infty}\leq\|\bm{\ell}^{\text{reg}}(\bm{y})\|\leq\|\bm{\ell}(\bm{y})\|\leq\sqrt{p}\,\|\bm{\ell}(\bm{y})\|_{\infty}\leq\sqrt{p}\,\Lambda,

(4.31)

and we are done. ∎

We are now in a position to construct our fully linear error bounds. We will begin by proving a bound on the size of the model Hessian, which requires the following technical result.

Lemma 4.3.

Fix $t>1$ and $c_{1},c_{2}\geq 0$ . If $a,b\in\mathbb{R}$ satisfy $|a+b|\leq c_{1}$ and $|ta+b|\leq c_{2}$ then $|a|\leq(c_{1}+c_{2})/(t-1)$ and $|b|\leq(tc_{1}+c_{2})/(t-1)$ .

Proof.

We first prove the bound on $|a|$ . To find a contradiction, first suppose that $a>(c_{1}+c_{2})/(t-1)$ and so $a>0$ . Since $a+b\geq-c_{1}$ we have $b\geq-c_{1}-a$ , which means $ta+b\geq(t-1)a-c_{1}>c_{2}$ , contradicting $|ta+b|\leq c_{2}$ . Instead, suppose that $a<-(c_{1}+c_{2})/(t-1)$ and so $a<0$ . Then $a+b\leq c_{1}$ means $b\leq c_{1}-a$ and so $ta+b\leq(t-1)a+c_{1}<-c_{2}$ , again contradicting $|ta+b|\leq c_{2}$ .

The bound on $|b|$ follows from similar reasoning. First suppose that $b>(tc_{1}+c_{2})/(t-1)$ and so $b>0$ . Since $a+b\leq c_{1}$ we have $a\leq c_{1}-b$ and so $ta+b\leq tc_{1}-(t-1)b<-c_{2}$ , contradicting $|ta+b|\leq c_{2}$ . Instead suppose that $b<-(tc_{1}+c_{2})/(t-1)$ and so $b<0$ . Then since $a+b\geq-c_{1}$ we have $a\geq-c_{1}-b$ and so $ta+b\geq-tc_{1}-(t-1)b>c_{2}$ , again contradicting $|ta+b|\leq c_{2}$ . ∎

We can now give our bound on the model Hessian. In the existing (unconstrained) theory, this is presented as a bound on $\|H\|$ [16, Theorem 5.7], but here we only need to consider specific Rayleigh-type quotients.

Lemma 4.4.

Suppose $f$ and $C$ satisfy Assumptions 2.1 and 2.2 respectively. Then if $\bm{x}\in C$ , $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ is $\Lambda$ -poised for underdetermined quadratic interpolation in $B(\bm{x},\Delta)\cap C$ and $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ for all $t=1,\ldots,p$ and some $\beta>0$ , then the model $m$ generated by (4.2) has Hessian $H$ satisfying

\displaystyle\max_{s,t=1,\ldots,p}\frac{|(\bm{y}_{s}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})|}{\beta^{2}\min(\Delta,1)^{2}}\leq\kappa_{H}:=L_{\nabla f}p\left[8\Lambda\beta^{2}+36\Lambda\beta+58\Lambda+6\right].

(4.32)

Proof.

First, fix $u\in\{1,\ldots,p\}$ and consider the associated Lagrange polynomial

\displaystyle\ell_{u}(\bm{y})=c_{u}+\bm{g}_{u}^{T}(\bm{y}-\bm{x})+\frac{1}{2}(\bm{y}-\bm{x})^{T}H_{u}(\bm{y}-\bm{x}).

(4.33)

Additionally, fix $s,t\in\{1,\ldots,p\}$ , and we will first provide a bound on

\displaystyle|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})|.

(4.34)

To do this, we consider the value of $\ell_{u}$ at five different points: $\bm{x}$ , $\bm{y}_{s}$ , $\bm{y}_{t}$ , plus

\displaystyle\hat{\bm{y}}_{s}:=\bm{x}+\frac{1}{\hat{\beta}}(\bm{y}_{s}-\bm{x}),\qquad\text{and}\qquad\hat{\bm{y}}_{t}:=\bm{x}+\frac{1}{\hat{\beta}}(\bm{y}_{t}-\bm{x}),

(4.35)

where $\hat{\beta}:=\max(\beta,2)$ . Since $\|\bm{y}_{s}-\bm{x}\|\leq\beta\min(\Delta,1)$ we have $\|\hat{\bm{y}}_{s}-\bm{x}\|\leq\frac{\beta}{\hat{\beta}}\min(\Delta,1)\leq\min(\Delta,1)$ , and $\hat{\bm{y}}_{s}\in C$ by convexity, since $\bm{x},\bm{y}_{s}\in C$ . Similarly, we also have $\hat{\bm{y}}_{t}\in B(\bm{x},\min(\Delta,1))\cap C$ . For these points, from the definition of $\Lambda$ -poisedness we know $|\ell_{u}(\bm{x})|,|\ell_{u}(\hat{\bm{y}}_{s})|,|\ell_{u}(\hat{\bm{y}}_{t})|\leq\Lambda$ and by definition of Lagrange polynomials we have $\ell_{u}(\bm{y}_{s}),\ell_{u}(\bm{y}_{t})\in\{0,1\}$ and so $|\ell_{u}(\bm{y}_{s})|,|\ell_{u}(\bm{y}_{t})|\leq 1$ .

From $|\ell_{u}(\bm{x})|\leq\Lambda$ we have $|c_{u}|\leq\Lambda$ , and so $|\ell_{u}(\hat{\bm{y}}_{s})|\leq\Lambda$ implies

$\displaystyle\Lambda$	$\displaystyle\geq\left\|c_{u}+\bm{g}_{u}^{T}(\hat{\bm{y}}_{s}-\bm{x})+\frac{1}{2}(\hat{\bm{y}}_{s}-\bm{x})^{T}H_{u}(\hat{\bm{y}}_{s}-\bm{x})\right\|,$	(4.36)
	$\displaystyle=\left\|c_{u}+\frac{1}{\hat{\beta}}\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})+\frac{1}{2\hat{\beta}^{2}}(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\right\|,$	(4.37)
	$\displaystyle\geq\frac{1}{\hat{\beta}^{2}}\left\|\hat{\beta}\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})+\frac{1}{2}(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\right\|-\|c_{u}\|,$	(4.38)

where the last line follows from the reverse triangle inequality. Together with $|\ell_{u}(\bm{y}_{s})|\leq 1$ , we get

	$\displaystyle\left\|\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})+\frac{1}{2}(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\right\|$	$\displaystyle\leq 1+\|c_{u}\|\leq\Lambda+1,$		(4.39)
	$\displaystyle\left\|\hat{\beta}\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})+\frac{1}{2}(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\right\|$	$\displaystyle\leq 2\Lambda\hat{\beta}^{2}.$		(4.40)

Since $\hat{\beta}>1$ by definition, applying Lemma 4.3 we conclude

	$\displaystyle\left\|\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})\right\|$	$\displaystyle\leq\frac{\Lambda+1+2\Lambda\hat{\beta}^{2}}{\hat{\beta}-1},$		(4.41)
	$\displaystyle\frac{1}{2}\left\|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\right\|$	$\displaystyle\leq\frac{\hat{\beta}(\Lambda+1)+2\Lambda\hat{\beta}}{\hat{\beta}-1}.$		(4.42)

Since $s$ was arbitrary, the same inequalities hold with $\bm{y}_{s}$ replaced by $\bm{y}_{t}$ .

Now, consider the point

\displaystyle\hat{\bm{y}}_{s,t}:=\bm{x}+\frac{1}{2}(\hat{\bm{y}}_{s}-\bm{x})+\frac{1}{2}(\hat{\bm{y}}_{t}-\bm{x}).

(4.43)

Since $\bm{x},\hat{\bm{y}}_{s},\hat{\bm{y}}_{t}\in C$ , we have $\hat{\bm{y}}_{s,t}\in C$ , and also

\displaystyle\|\hat{\bm{y}}_{s,t}-\bm{x}\|\leq\frac{1}{2}\|\hat{\bm{y}}_{s}-\bm{x}\|+\frac{1}{2}\|\hat{\bm{y}}_{t}-\bm{x}\|\leq\min(\Delta,1),

(4.44)

and so $|\ell_{u}(\hat{\bm{y}}_{s,t})|\leq\Lambda$ . Written in full, this is

	$\displaystyle\left\|c_{u}+\frac{1}{2\hat{\beta}}\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})+\frac{1}{2\hat{\beta}}\bm{g}_{u}^{T}(\bm{y}_{t}-\bm{x})\right.$
	$\displaystyle\qquad\qquad\left.+\frac{1}{2}\left(\frac{1}{2\hat{\beta}}(\bm{y}_{s}-\bm{x})+\frac{1}{2\hat{\beta}}(\bm{y}_{t}-\bm{x})\right)^{T}H_{u}\left(\frac{1}{2\hat{\beta}}(\bm{y}_{s}-\bm{x})+\frac{1}{2\hat{\beta}}(\bm{y}_{t}-\bm{x})\right)\right\|\leq\Lambda.$		(4.45)

That is,

	$\displaystyle\frac{1}{4\hat{\beta}^{2}}\|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})\|$	$\displaystyle\leq\Lambda+\|c_{u}\|+\frac{1}{2\hat{\beta}}\left(\|\bm{g}_{u}^{T}(\bm{y}_{s}-\bm{x})\|+\|\bm{g}_{u}^{T}(\bm{y}_{t}-\bm{x})\|\right)$
		$\displaystyle\qquad\qquad+\frac{1}{8\hat{\beta}^{2}}\left(\|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{s}-\bm{x})\|+\|(\bm{y}_{t}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})\|\right).$		(4.46)

Applying $|c_{u}|\leq\Lambda$ , (4.41) and (4.42), we conclude

\displaystyle\frac{1}{4\hat{\beta}^{2}}|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})|

\displaystyle\leq 2\Lambda+\frac{\Lambda+1+2\Lambda\hat{\beta}^{2}}{\hat{\beta}(\hat{\beta}-1)}+\frac{\hat{\beta}(\Lambda+1)+2\Lambda\hat{\beta}}{2\hat{\beta}^{2}(\hat{\beta}-1)},

(4.47)

\displaystyle|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})|

\displaystyle\leq\tilde{\kappa},

(4.48)

for all $s,t,u=1,\ldots,p$ , where

\displaystyle\tilde{\kappa}:=8\Lambda\hat{\beta}^{2}+\frac{4\Lambda\hat{\beta}+4\hat{\beta}+8\Lambda\hat{\beta}^{3}}{\hat{\beta}-1}+\frac{2\hat{\beta}(\Lambda+1)+4\Lambda\hat{\beta}}{\hat{\beta}-1}=8\Lambda\hat{\beta}^{2}+\frac{8\Lambda\hat{\beta}^{3}+10\Lambda\hat{\beta}+6\hat{\beta}}{\hat{\beta}-1}.

(4.49)

More simply, we have

$\displaystyle\tilde{\kappa}$	$\displaystyle=8\Lambda\hat{\beta}^{2}+\frac{8\Lambda(\hat{\beta}-1)^{3}+24\Lambda(\hat{\beta}-1)^{2}+(34\Lambda+6)(\hat{\beta}-1)+18\Lambda+6}{\hat{\beta}-1},$	(4.50)
	$\displaystyle\leq 8\Lambda(\beta+2)^{2}+8\Lambda(\beta+1)^{2}+24\Lambda(\beta+1)+34\Lambda+6+18\Lambda+6,$	(4.51)
	$\displaystyle=16\Lambda\beta^{2}+72\Lambda\beta+116\Lambda+12,$	(4.52)

using $2\leq\hat{\beta}=\max(\beta,2)\leq\beta+2$ to get the inequality. Now let us turn our attention to the model $m$ (4.2). Note that we may add/subtract a linear function to $f(\bm{y})$ without changing $H$ (since this just changes $c$ and $\bm{g}$ in (4.2)). So, we consider the model $\widetilde{m}$ generated by interpolation to $\widetilde{f}(\bm{y}):=f(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})$ . From (4.27) we may write

\displaystyle H=\widetilde{H}=\sum_{u=1}^{p}\widetilde{f}(\bm{y}_{u})H_{u},

(4.53)

where $H_{u}$ are the Hessians of the Lagrange polynomials. From Assumption 2.1, we also have $|\widetilde{f}(\bm{y})|\leq\frac{L_{\nabla f}}{2}\|\bm{y}-\bm{x}\|^{2}$ . Then for any $s,t=1,\ldots,p$ we conclude

$\displaystyle\|(\bm{y}_{s}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})\|$	$\displaystyle\leq\sum_{u=1}^{p}\|\widetilde{f}(\bm{y}_{u})\|\cdot\|(\bm{y}_{s}-\bm{x})^{T}H_{u}(\bm{y}_{t}-\bm{x})\|,$	(4.54)
	$\displaystyle\leq\sum_{u=1}^{p}\frac{L_{\nabla f}}{2}\tilde{\kappa}\\|\bm{y}_{u}-\bm{x}\\|^{2},$	(4.55)
	$\displaystyle\leq\frac{L_{\nabla f}}{2}\tilde{\kappa}p\beta^{2}\min(\Delta,1)^{2}.$	(4.56)

and we are done after applying (4.52). ∎

Remark 4.5.

In the unconstrained case, [16, Theorem 5.7] gives the bound

\displaystyle\|H\|\leq 4(p+1)\sqrt{\frac{(n+1)(n+2)}{2}}\>L_{\nabla f}\Lambda\max(1,\Delta_{\max}^{2}),

(4.57)

where $\Delta_{\max}$ is an upper bound for $\Delta$ . This result assumes $\|\bm{y}_{t}-\bm{x}\|\leq\Delta$ , and so to match our assumptions it requires $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ where $\beta=\Delta_{\max}$ . Thus we see that Lemma 4.4 improves on [16, Theorem 5.7] by a factor of $\mathcal{O}(n)$ .

Next, we can prove an analogous result to Lemma 3.3.

Lemma 4.6.

Suppose the assumptions of Lemma 4.4 are satisfied. Then the model $m$ generated by (4.2) satisfies

\displaystyle\left|c+\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right|\leq\frac{1}{2}p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2},

(4.58)

for all $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ . Furthermore if $\bm{x}\in C$ , we have

\displaystyle|c-f(\bm{x})|\leq\frac{1}{2}p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2},

(4.59)

and

\displaystyle|(\bm{g}-\nabla f(\bm{x}))^{T}(\bm{y}-\bm{x})|\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2},

(4.60)

for all $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ .

Proof.

Fix $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ . From Lemma 4.2, our interpolation set is $\sqrt{p}\,\Lambda$ -poised for linear regression. Therefore there exist constants $\{\alpha_{t}(\bm{y})\}_{t=1}^{p}$ such that

\displaystyle\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}=\sum_{t=1}^{p}\alpha_{t}(\bm{y})\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix},

(4.61)

where $|\alpha_{t}(\bm{y})|=|\ell_{t}^{\text{reg}}(\bm{y})|\leq\sqrt{p}\,\Lambda$ (provided the minimal norm solution is taken as in (3.15)). Thus we have

	$\displaystyle\left\|c+\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\right\|$
	$\displaystyle\qquad\qquad=\left\|\begin{bmatrix}c-f(\bm{x})\\ \bm{g}-\nabla f(\bm{x})\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}-\bm{x}\end{bmatrix}\right\|,$		(4.62)
	$\displaystyle\qquad\qquad\leq\sum_{t=1}^{p}\sqrt{p}\>\Lambda\left\|\begin{bmatrix}c-f(\bm{x})\\ \bm{g}-\nabla f(\bm{x})\end{bmatrix}^{T}\begin{bmatrix}1\\ \bm{y}_{t}-\bm{x}\end{bmatrix}\right\|,$		(4.63)
	$\displaystyle\qquad\qquad=\sqrt{p}\>\Lambda\sum_{t=1}^{p}\left\|m(\bm{y}_{t})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x})-\frac{1}{2}(\bm{y}_{t}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})\right\|,$		(4.64)
	$\displaystyle\qquad\qquad\leq\sqrt{p}\>\Lambda\sum_{t=1}^{p}\left[\left\|f(\bm{y}_{t})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}_{t}-\bm{x})\right\|+\frac{1}{2}\left\|(\bm{y}_{t}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})\right\|\right],$		(4.65)
	$\displaystyle\qquad\qquad\leq p^{3/2}\Lambda\left(\frac{L_{\nabla f}}{2}\beta^{2}\min(\Delta,1)^{2}+\frac{\kappa_{H}}{2}\beta^{2}\min(\Delta,1)^{2}\right),$		(4.66)

using Assumption 2.1 and Lemma 4.4 in the last line, and we have (4.58). To get (4.59), we just take $\bm{y}=\bm{x}$ in (4.58). Finally, to get (4.60), we use

\displaystyle|(\bm{g}-\nabla f(\bm{x}))^{T}(\bm{y}-\bm{x})|\leq|c+\bm{g}^{T}(\bm{y}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})|+|c-f(\bm{x})|,

(4.67)

together with (4.58) and (4.59). ∎

Finally, we can give our fully linear error bounds for underdetermined quadratic interpolation models.

Theorem 4.7.

\displaystyle\kappa_{\textnormal{ef}}=\frac{1}{2}L_{\nabla f}+\frac{3}{2}p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}+\frac{1}{2}p\Lambda^{2}\kappa_{H}\beta^{2},\quad\text{and}\quad\kappa_{\textnormal{eg}}=p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2},

(4.68)

in (2.4), where $\kappa_{H}$ is defined in Lemma 4.4.

Proof.

First fix $\bm{y}\in B(\bm{x},\Delta)\cap C$ and we will derive $\kappa_{\textnormal{ef}}$ by considering the cases $\Delta\geq 1$ and $\Delta<1$ separately. If $\Delta\geq 1$ then $\hat{\bm{y}}:=\bm{x}+\Delta^{-1}(\bm{y}-\bm{x})$ is in $B(\bm{x},1)=B(\bm{x},\min(\Delta,1))$ and $\hat{\bm{y}}\in C$ since $C$ is convex. We then apply Lemma 4.6 to get

$\displaystyle\|(\bm{y}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle=\Delta\|(\hat{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|,$	(4.69)
	$\displaystyle\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta\min(\Delta,1)^{2},$	(4.70)
	$\displaystyle\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta^{2},$	(4.71)

where the last inequality follows from $\min(\Delta,1)^{2}=1\leq\Delta$ . Also, from (4.61) we have $\hat{\bm{y}}-\bm{x}=\sum_{t=1}^{p}\alpha_{t}(\hat{\bm{y}})(\bm{y}_{t}-\bm{x})$ with $|\alpha_{t}(\bm{y})|\leq\sqrt{p}\,\Lambda$ , and so from Lemma 4.4,

$\displaystyle\|(\bm{y}-\bm{x})^{T}H(\bm{y}-\bm{x})\|$	$\displaystyle=\Delta^{2}\|(\hat{\bm{y}}-\bm{x})^{T}H(\hat{\bm{y}}-\bm{x})\|,$	(4.72)
	$\displaystyle\leq\Delta^{2}\sum_{s,t=1}^{p}\left\|\alpha_{s}(\hat{\bm{y}})\alpha_{t}(\hat{\bm{y}})(\bm{y}_{s}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})\right\|,$	(4.73)
	$\displaystyle\leq p\Lambda^{2}\kappa_{H}\beta^{2}\Delta^{2}\min(\Delta,1)^{2},$	(4.74)
	$\displaystyle=p\Lambda^{2}\kappa_{H}\beta^{2}\Delta^{2},$	(4.75)

again using $\Delta\geq 1$ in the last line.

Instead, if $\Delta<1$ then $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ and we apply Lemma 4.6 directly to get

	$\displaystyle\|(\bm{y}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle\leq\frac{1}{2}p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2},$		(4.76)
		$\displaystyle=p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta^{2}.$		(4.77)

Also, using Lemma 4.4 we have

$\displaystyle\|(\bm{y}-\bm{x})^{T}H(\bm{y}-\bm{x})\|$	$\displaystyle\leq\sum_{t=1}^{p}\left\|\alpha_{s}(\bm{y})\alpha_{t}(\bm{y})(\bm{y}_{s}-\bm{x})^{T}H(\bm{y}_{t}-\bm{x})\right\|,$	(4.78)
	$\displaystyle\leq p\Lambda^{2}\kappa_{H}\beta^{2}\min(\Delta,1)^{2},$	(4.79)
	$\displaystyle=p\Lambda^{2}\kappa_{H}\beta^{2}\Delta^{2}.$	(4.80)

In either case, we use (4.71) and (4.75), or (4.77) and (4.80), with Assumption 2.1 and Lemma 4.6 to get

$\displaystyle\|f(\bm{y})-m(\bm{y})\|$	$\displaystyle\leq\|f(\bm{y})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|$
	$\displaystyle\qquad+\|c+\bm{g}^{T}(\bm{y}-\bm{x})+\frac{1}{2}(\bm{y}-\bm{x})^{T}H(\bm{y}-\bm{x})-f(\bm{x})-\nabla f(\bm{x})^{T}(\bm{y}-\bm{x})\|,$	(4.81)
	$\displaystyle\leq\frac{1}{2}L_{\nabla f}\Delta^{2}+\|c-f(\bm{x})\|+\|(\bm{y}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|+\frac{1}{2}\|(\bm{y}-\bm{x})^{T}H(\bm{y}-\bm{x})\|,$	(4.82)
	$\displaystyle\leq\frac{1}{2}L_{\nabla f}\Delta^{2}+\frac{1}{2}p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2}$
	$\displaystyle\qquad\qquad\qquad\qquad+p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta^{2}+\frac{1}{2}p\Lambda^{2}\kappa_{H}\beta^{2}\Delta^{2},$	(4.83)

and we get the value of $\kappa_{\textnormal{ef}}$ after $\min(\Delta,1)\leq\Delta$ .

	$\displaystyle\|(\widetilde{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\min(\Delta,1)^{2},$		(4.84)
		$\displaystyle\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta,$		(4.85)

since $\min(\Delta,1)^{2}=1\leq\Delta$ . Alternatively, if $\Delta<1$ then the convexity of $C$ implies that $\hat{\bm{y}}:=\bm{x}+\Delta(\widetilde{\bm{y}}-\bm{x})\in B(\bm{x},\Delta)\cap C=B(\bm{x},\min(\Delta,1))\cap C$ . Again we apply Lemma 4.6 and get

$\displaystyle\|(\widetilde{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|$	$\displaystyle=\Delta^{-1}\|(\hat{\bm{y}}-\bm{x})^{T}(\bm{g}-\nabla f(\bm{x}))\|,$	(4.86)
	$\displaystyle\leq p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta^{-1}\min(\Delta,1)^{2},$	(4.87)
	$\displaystyle=p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H})\beta^{2}\Delta.$	(4.88)

The value for $\kappa_{\textnormal{eg}}$ then follows from (4.85) and (4.88). ∎

Remark 4.8.

The fully linear error bounds for the unconstrained case [14, 16] give $\kappa_{\textnormal{ef}},\kappa_{\textnormal{eg}}=\mathcal{O}(p^{3/2}\Lambda(L_{\nabla f}+\kappa_{H}))$ , so Theorem 4.7 matches the bound for $\kappa_{\textnormal{eg}}$ up to a factor of $\beta^{2}$ (but with a value of $\kappa_{H}$ which is $\mathcal{O}(n)$ smaller than the unconstrained case; see Remark 4.5). Our value of $\kappa_{\textnormal{ef}}$ has an extra term of size $\mathcal{O}(p\Lambda^{2}\kappa_{H}\beta^{2})$ and so may be larger depending on the relative sizes of $p$ and $\Lambda$ .

5 Constructing $\bm{\Lambda}$ -Poised Quadratic Models

To meet all the requirements of Algorithm 1, we require procedures which can verify whether or not a given model is $C$ -fully linear, and if not, modify it to be $C$ -fully linear. From Theorem 4.7 it is clear that it suffices to verify/ensure that $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ is $\Lambda$ -poised for underdetermined quadratic interpolation in $B(\bm{x},\Delta)\cap C$ and $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ for all $t=0,\ldots,p$ .

In the first case—verifying whether not the interpolation set is $\Lambda$ -poised—we may follow the approach in [22, Section 4.3.1] and simply maximize/minimize each $\ell_{t}$ in $B(\bm{x},\min(\Delta,1))\cap C$ . Checking $\|\bm{y}_{t}-\bm{x}\|\leq\beta\min(\Delta,1)$ is straightforward. Although maximizing/minimizing $\ell_{t}$ , a possibly nonconvex quadratic function, in $B(\bm{x},\min(\Delta,1))\cap C$ may not be straightforward depending on the structure of $C$ , we only need to know if the solution is above/below $\Lambda$ and not the exact solution. Existing solvers, including potentially global optimization solvers, are sufficient for this.

We now consider the situation where $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ is not $\Lambda$ -poised, and we wish to construct a new interpolation set which is $\Lambda$ -poised. There are two possible scenarios:

•

The associated matrix $F$ in (4.20) is not invertible; or
•

The matrix $F$ is invertible but $|\ell_{t}(\bm{y})|>\Lambda$ for some $t=1,\ldots,p$ and $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ .

We will describe how to handle both of these situations in Section 5.2, but we will first need a technical result about how the matrix $F$ in (4.20) changes as interpolation points are changed.

5.1 Updating the matrix $\bm{F}$

In this section we analyze changes in the determinant of the matrix $F$ as interpolation points are updated. These results are based on ideas from [26, Section 4], which considers changes in $F^{-1}$ from updating interpolation points, rather than $\det(F)$ .

Lemma 5.1.

If $A$ is a symmetric, invertible matrix and we form $\widetilde{A}$ by changing the $t$ -th row and column of $A$ to $\widetilde{\bm{v}}$ , then

\displaystyle\det(\widetilde{A})=\left[(\bm{e}_{t}^{T}A^{-1}\widetilde{\bm{v}})^{2}+(\bm{e}_{t}^{T}A^{-1}\bm{e}_{t})\widetilde{\bm{v}}^{T}(\bm{e}_{t}-A^{-1}\widetilde{\bm{v}})\right]\det(A).

(5.1)

Proof.

Denote the $t$ -th row and column of $A$ by $\bm{v}:=A\bm{e}_{t}$ . We first note that $\widetilde{A}$ can be written as a symmetric rank-3 update of $A$ :

\displaystyle\widetilde{A}=A+(\widetilde{\bm{v}}-\bm{v})\bm{e}_{t}^{T}+\bm{e}_{t}(\widetilde{\bm{v}}-\bm{v})^{T}-(\widetilde{v}_{t}-v_{t})\bm{e}_{t}\bm{e}_{t}^{T}.

(5.2)

The first two update terms add $\widetilde{\bm{v}}-\bm{v}$ to the $t$ -th column and row respectively, and the last term corrects for the double update of $A_{t,t}$ from the first two updates. We can then apply the matrix determinant lemma²²2That is, $\det(A+UV^{T})=\det(I+V^{T}A^{-1}U)\det(A)$ if $A$ is invertible, e.g. [21, Corollary 18.1.4]. to get

	$\displaystyle\frac{\det(\widetilde{A})}{\det(A)}$	$\displaystyle=\det\left(\begin{bmatrix}1+\bm{e}_{t}^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})&\bm{e}_{t}^{T}A^{-1}\bm{e}_{t}&\bm{e}_{t}^{T}A^{-1}\bm{e}_{t}\\ (\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})&1+(\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}\bm{e}_{t}&(\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}\bm{e}_{t}\\ -(\widetilde{v}_{t}-v_{t})\bm{e}_{t}^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})&-(\widetilde{v}_{t}-v_{t})\bm{e}_{t}^{T}A^{-1}\bm{e}_{t}&1-(\widetilde{v}_{t}-v_{t})\bm{e}_{t}^{T}A^{-1}\bm{e}_{t}\end{bmatrix}\right),$		(5.3)
		$\displaystyle=\det\left(\begin{bmatrix}1+\tau&\alpha&\alpha\\ \beta&1+\tau&\tau\\ -\gamma\tau&-\gamma\alpha&1-\gamma\alpha\end{bmatrix}\right),$		(5.4)

where $\alpha:=\bm{e}_{t}^{T}A^{-1}\bm{e}_{t}$ , $\beta:=(\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})$ , $\gamma:=\widetilde{v}_{t}-v_{t}$ and $\tau:=\bm{e}_{t}^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})=(\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}\bm{e}_{t}$ since $A^{-1}$ is symmetric. Directly computing, we get

$\displaystyle\frac{\det(\widetilde{A})}{\det(A)}$	$\displaystyle=(1+\tau)\left[(1+\tau)(1-\gamma\alpha)+\tau\gamma\alpha\right]-\beta\left[\alpha(1-\gamma\alpha)+\gamma\alpha^{2}\right]-\gamma\tau\left[\alpha\tau-\alpha(1+\tau)\right],$	(5.5)
	$\displaystyle=(1+\tau)\left[1+\tau-\gamma\alpha\right]-\beta\alpha+\gamma\tau\alpha,$	(5.6)
	$\displaystyle=(1+\tau)^{2}-\alpha(\gamma+\beta).$	(5.7)

Finally, since $\bm{v}=A\bm{e}_{t}$ , we get

\displaystyle\tau=\bm{e}_{t}^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v})=\bm{e}_{t}^{T}A^{-1}\widetilde{\bm{v}}-\bm{e}_{t}^{T}\underbrace{A^{-1}\bm{v}}_{=\bm{e}_{t}}=\bm{e}_{t}^{T}A^{-1}\widetilde{\bm{v}}-1,

(5.8)

and

$\displaystyle\gamma+\beta$	$\displaystyle=(\widetilde{\bm{v}}-\bm{v})^{T}\bm{e}_{t}+(\widetilde{\bm{v}}-\bm{v})^{T}A^{-1}(\widetilde{\bm{v}}-\bm{v}),$	(5.9)
	$\displaystyle=\widetilde{\bm{v}}^{T}\bm{e}_{t}-\bm{v}^{T}\bm{e}_{t}+\widetilde{\bm{v}}^{T}A^{-1}\widetilde{\bm{v}}-2\widetilde{\bm{v}}^{T}A^{-1}\bm{v}+\bm{v}^{T}A^{-1}\bm{v},$	(5.10)
	$\displaystyle=\widetilde{\bm{v}}^{T}\bm{e}_{t}-\bm{v}^{T}\bm{e}_{t}+\widetilde{\bm{v}}^{T}A^{-1}\widetilde{\bm{v}}-2\widetilde{\bm{v}}^{T}\bm{e}_{t}+\bm{v}^{T}\bm{e}_{t},$	(5.11)
	$\displaystyle=\widetilde{\bm{v}}^{T}(A^{-1}\widetilde{\bm{v}}-\bm{e}_{t}).$	(5.12)

Our result then follows from (5.7) combined with the definition of $\alpha$ , (5.8) and (5.12). ∎

Theorem 5.2.

Suppose the interpolation set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\in\mathbb{R}^{n}$ is such that the matrix $F$ in (4.20) is invertible. If the point $\bm{y}_{t}$ for some $t\in\{1,\ldots,p\}$ is changed to a new point $\bm{y}\in\mathbb{R}^{n}$ , then the new interpolation set yields a new matrix $\widetilde{F}$ in (4.20) which satisfies

\displaystyle|\det(\widetilde{F})|\geq\ell_{t}(\bm{y})^{2}|\det(F)|.

(5.13)

Proof.

If $\ell_{t}(\bm{y})=0$ , then the result holds trivially, so we assume without loss of generality that $\ell_{t}(\bm{y})\neq 0$ .

Changing the interpolation point $\bm{y}_{t}$ to $\bm{y}$ requires replacing the $t$ -th row and column of $F$ by $\widetilde{\bm{v}}$ , where


$\displaystyle\widetilde{\bm{v}}_{s}$	$\displaystyle=\frac{1}{2}\left[(\bm{y}-\bm{x})^{T}(\bm{y}_{s}-\bm{x})\right]^{2},\qquad\text{if $s\in\{1,\ldots,p\}$ except for $s=t$},$	(5.14a)
$\displaystyle\widetilde{\bm{v}}_{t}$	$\displaystyle=\frac{1}{2}\left[(\bm{y}-\bm{x})^{T}(\bm{y}-\bm{x})\right]^{2},$	(5.14b)
$\displaystyle\widetilde{\bm{v}}_{p+2}$	$\displaystyle=1,$	(5.14c)
$\displaystyle\widetilde{\bm{v}}_{p+2+i}$	$\displaystyle=[\bm{y}-\bm{x}]_{i},\qquad\forall i=1,\ldots,n.$	(5.14d)

That is, $\widetilde{\bm{v}}$ is the same as $\bm{\phi}(\bm{y})$ from (4.25) except for the $t$ -th entry. Specifically, we have

\displaystyle\widetilde{\bm{v}}

\displaystyle=\bm{\phi}(\bm{y})+\eta_{t}\bm{e}_{t},

(5.15)

where

\displaystyle\eta_{t}

\displaystyle=\frac{1}{2}\left[(\bm{y}-\bm{x})^{T}(\bm{y}-\bm{x})\right]^{2}-\frac{1}{2}\left[(\bm{y}-\bm{x})^{T}(\bm{y}_{t}-\bm{x})\right]^{2}=\frac{1}{2}\|\bm{y}-\bm{x}\|^{4}-\bm{\phi}(\bm{y})^{T}\bm{e}_{t}.

(5.16)

Given this, Lemma 5.1 gives us (denoting $\alpha_{t}:=\bm{e}_{t}^{T}F^{-1}\bm{e}_{t}$ for notational convenience)

$\displaystyle\frac{\det(\widetilde{F})}{\det(F)}$	$\displaystyle=\left(\bm{e}_{t}^{T}F^{-1}\bm{\phi}(\bm{y})+\eta_{t}\bm{e}_{t}^{T}F^{-1}\bm{e}_{t}\right)^{2}+(\bm{e}_{t}^{T}F^{-1}\bm{e}_{t})(\bm{\phi}(\bm{y})+\eta_{t}\bm{e}_{t})^{T}(\bm{e}_{t}-F^{-1}(\bm{\phi}(\bm{y})+\eta_{t}\bm{e}_{t})),$	(5.17)
	$\displaystyle=\left(\ell_{t}(\bm{y})+\eta_{t}\alpha_{t}\right)^{2}+\alpha_{t}\left(\bm{\phi}(\bm{y})^{T}\bm{e}_{t}-\bm{\phi}(\bm{y})^{T}F^{-1}\bm{\phi}(\bm{y})-\eta_{t}\ell_{t}(\bm{y})+\eta_{t}-\eta_{t}\ell_{t}(\bm{y})-\eta_{t}^{2}\alpha_{t}\right),$	(5.18)
	$\displaystyle=\ell_{t}(\bm{y})^{2}+\alpha_{t}\left(\bm{\phi}(\bm{y})^{T}\bm{e}_{t}-\bm{\phi}(\bm{y})^{T}F^{-1}\bm{\phi}(\bm{y})+\eta_{t}\right).$	(5.19)

So from the definition of $\eta_{t}$ , we get

\displaystyle\det(\widetilde{F})

\displaystyle=\left(\ell_{t}(\bm{y})^{2}+\alpha_{t}\beta_{t}\right)\>\det(F),

(5.20)

where $\beta_{t}:=\frac{1}{2}\|\bm{y}-\bm{x}\|^{4}-\bm{\phi}(\bm{y})^{T}F^{-1}\bm{\phi}(\bm{y})$ . Finally, [26, Lemma 2] gives us $\alpha_{t},\beta_{t}\geq 0$ provided $\ell_{t}(\bm{y})\neq 0$ , which gives us the result.³³3That $\ell_{t}(\bm{y})\neq 0$ is not explicitly assumed in the statement of [26, Lemma 2] but is required as the proof relies on [26, Theorem, p. 196]. ∎

Remark 5.3.

Theorem 5.2 does not require any of the points $\bm{y}_{1},\ldots,\bm{y}_{p}$ or $\bm{y}$ to be in $C$ . We will use this when constructing interpolation sets for which $F$ is invertible (Algorithm 2).

Remark 5.4.

A similar analysis to the above holds for linear interpolation ( $p=n$ ), where Theorem 5.2 simplifies to $|\det(\widetilde{F})|=|\ell_{t}(\bm{y})|\cdot|\det(F)|$ as shown in [8, Section 4.2]. This allows us to use an approach like Algorithm 3 below to construct $\Lambda$ -poised sets for linear interpolation, a simpler method than the one outlined in [22, Section 4.3].

5.2 Constructing $\bm{\Lambda}$ -Poised Sets

We begin by describing how to construct an initial interpolation set, contained in $B(\bm{x},\min(\Delta,1))\cap C$ for which $F$ is invertible. The full process for this is given in Algorithm 2. We first generate an initial set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ in $B(\bm{x},\min(\Delta,1))$ , but not necessarily in $C$ , for which $F$ is invertible. We can do this using the approach outlined in [29, Section 2], which chooses points of the form $\bm{x}\pm\min(\Delta,1)\bm{e}_{t}$ or $\bm{x}\pm\min(\Delta,1)\bm{e}_{s}\pm\min(\Delta,1)\bm{e}_{t}$ for some $s,t\in\{1,\ldots,n\}$ in a specific way. We then replace any of these points which are not in $C$ by new points in $B(\bm{x},\min(\Delta,1))\cap C$ while maintaining an invertibile $F$ .

1:Trust region center

\bm{x}\in C

and radius

\Delta>0

, number of interpolation points

p\in\{n+2,\ldots,(n+1)(n+2)/2\}

2:Generate an initial collection of

p

points

\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset B(\bm{x},\min(\Delta,1))

using the method from [29, Section 2].

3:while there exists

t\in\{1,\ldots,p\}

with

\bm{y}_{t}\notin C

4: Calculate Lagrange polynomial

\ell_{t}

and find

\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C

with

\ell_{t}(\bm{y})\neq 0

5: Replace the interpolation point

\bm{y}_{t}

\bm{y}

6:end while

7:return interpolation set

\{\bm{y}_{1},\ldots,\bm{y}_{p}\}

Algorithm 2 Construct an underdetermined quadratic interpolation set with invertible

F

(4.20).

Theorem 5.5.

For any $\bm{x}\in C$ , $\Delta>0$ and $p\in\{n+2,\ldots.(n+1)(n+2)/2\}$ , the output of Algorithm 2 is an interpolation set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset B(\bm{x},\min(\Delta,1))\cap C$ for which $F$ (4.20) is invertible.

Proof.

Firstly, Algorithm 2 terminates after at most $p$ iterations, since an index $t$ can only ever be selected once (after which the new $\bm{y}_{t}\in C$ ).

Next, all points in the output set were either originally generated in line 2, in which case they must have originally been in both $B(\bm{x},\min(\Delta,1))$ (from the construction in line 2) and $C$ (in order to never be replaced in the main loop), or were a $\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C$ chosen in line 4. Hence $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset B(\bm{x},\min(\Delta,1))\cap C$ .

Finally, the initial collection $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ from [29, Section 2] generates an invertible $F$ , with an explicit formula for $F^{-1}$ given in [29, eqns. 2.11–2.15]. Hence the initial $F$ satisfies $|\det(F)|>0$ . Then at each iteration replace $\bm{y}_{t}$ with a point $\bm{y}$ for which $\ell_{t}(\bm{y})\neq 0$ , and so the new $F$ also satisfies $|\det(F)|>0$ by Theorem 5.2. Thus the output interpolation set has an invertible $F$ . ∎

1:Trust region center

\bm{x}\in C

and radius

\Delta>0

, number of interpolation points

p\in\{n+2,\ldots,(n+1)(n+2)/2\}

, optional initial interpolation set

\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset C

, poisedness constant

\Lambda>1

2:if initial interpolation set not provided or yields singular

F

(4.20) or

\max_{t=1,\ldots,p}\|\bm{y}_{t}-\bm{x}\|>\min(\Delta,1)

then

3: Set

\{\bm{y}_{1},\ldots,\bm{y}_{p}\}

to be the output of Algorithm 2.

4:end if

5:while

\max_{\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C}\max_{t=1,\ldots,p}|\ell_{t}(\bm{y})|>\Lambda

6: Find

t\in\{1,\ldots,p\}

and

\bm{y}\in B(\bm{x},\min(\Delta,1))\cap C

with

|\ell_{t}(\bm{y})|>\Lambda

7: Replace the interpolation point

\bm{y}_{t}

\bm{y}

and recalculate Lagrange polynomials.

8:end while

9:return interpolation set

\{\bm{y}_{1},\ldots,\bm{y}_{p}\}

Algorithm 3 Construct a

\Lambda

-poised underdetermined quadratic interpolation set.

We can now present our algorithm to generate a $\Lambda$ -poised set for minimum Frobenius norm interpolation, given in Algorithm 3. This algorithm is very similar to existing methods for unconstrained linear/quadratic (e.g. [16, Algorithm 6.3]) or constrained linear poisedness ([22, Algorithm 4.2]). Essentially we first ensure we have an interpolation set with invertible $F$ (so we can construct the associated Lagrange polynomials), and we then iteratively find indices $t$ and points $\bm{y}$ with $|\ell_{t}(\bm{y})|>\Lambda$ , replacing $\bm{y}_{t}$ with $\bm{y}$ .

Theorem 5.6.

For any $\bm{x}\in C$ , $\Delta>0$ , $p\in\{n+2,\ldots.(n+1)(n+2)/2\}$ and $\Lambda>1$ , Algorithm 3 produces an interpolation set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ which is $\Lambda$ -poised for minimum Frobenius norm interpolation in $B(\bm{x},\Delta)\cap C$ and for which $\|\bm{y}_{t}-\bm{x}\|\leq\min(\Delta,1)$ for all $t=1,\ldots,p$ .

Proof.

If Algorithm 3 terminates, then by the definition of the loop check it must satisfy the requirements of the theorem, so it suffices to show that the ‘while’ loop terminates.

At the start of the ‘while’ loop in line 5, we have that $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}\subset B(\bm{x},\min(\Delta,1))\cap C$ and $F$ is invertible (either because the ‘if’ statement in line 2 evaluates to false or by Theorem 5.5). Hence at the start of the ‘while’ loop we have $|\det(F)|>0$ . Let us call this value $d_{0}>0$ .

In each iteration of the ‘while’ loop, we increase $|\det(F)|$ by a factor of at least $\Lambda^{2}$ by Theorem 5.2. Hence after $k$ iterations of the loop, we have $|\det(F)|\geq\Lambda^{2k}d_{0}$ . However, since by construction our interpolation set is contained in the compact set $B(\bm{x},\min(\Delta,1))\cap C$ , there is a global upper bound $d_{\max}$ on $|\det(F)|$ (which is a continuous function of the interpolation points). Hence the loop must terminate after at most $\log(d_{\max}/d_{0})/(2\log\Lambda)$ iterations. ∎

Remark 5.7.

In the unconstrained case $C=\mathbb{R}^{n}$ , Theorem 5.6 and in particular its reliance on Theorem 5.2 clarifies the argument behind [16, Theorem 6.6].

Theorems 4.7 (with $\beta=1$ ) and 5.6 then give us the following.

Corollary 5.8.

Suppose $f$ and $C$ satisfy Assumptions 2.1 and 2.2 respectively. Then for any $\bm{x}\in C$ , $\Delta>0$ , $p\in\{n+2,\ldots.(n+1)(n+2)/2\}$ and $\Lambda>1$ , Algorithm 3 produces an interpolation set $\{\bm{y}_{1},\ldots,\bm{y}_{p}\}$ for which the model $m$ generated by (4.2) is fully linear in $B(\bm{x},\Delta)$ .

6 Conclusions and Future Work

In [22], a weaker notion of fully linear models was used to construct convergent algorithms for (1.1) (namely Algorithm 1), and it was shown that linear interpolation models constructed only with feasible points can yield such models. Here, we show that Algorithm 1 can also be implemented using underdetermined quadratic interpolation models in the minimum Frobenius sense of [26], which is a much more useful construction in practice. As an extra benefit, it also establishes the same results for linear regression models, as a strict generalization of the analysis in [22] that is independently useful in the case where only noisy objective evaluations are available.

The natural extension of this work is to the construction and analysis of second-order accurate models (“fully quadratric” in the language of [16]) built only using feasible points. A concrete implementation of Algorithm 1 with quadratic interpolation models would also be a useful way to extend the implementation from [22] (which only considers nonlinear least-squares objectives, for which linear interpolation models can be practically useful).

References

[1] C. Audet and W. Hare, Derivative-Free and Blackbox Optimization, Springer Series in Operations Research and Financial Engineering, Springer, Cham, Switzerland, 2017.
[2] F. Augustin and Y. M. Marzouk, A trust-region method for derivative-free nonlinear constrained stochastic optimization, arXiv preprint arXiv:1703.04156, (2017).
[3] A. Beck, First-Order Methods in Optimization, MOS-SIAM Series on Optimization, MOS/SIAM, Philadelphia, 2017.
[4] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems, Acta Numerica, 14 (2005), pp. 1–137.
[5] S. C. Billups, J. Larson, and P. Graf, Derivative-free optimization of expensive functions with computational error using weighted regression, SIAM Journal on Optimization, 23 (2013), pp. 27–53.
[6] C. Cartis, J. Fiala, B. Marteau, and L. Roberts, Improving the flexibility and robustness of model-based derivative-free optimization solvers, ACM Transactions on Mathematical Software, 45 (2019), pp. 32:1–41.
[7] C. Cartis, N. I. M. Gould, and P. L. Toint, An adaptive cubic regularization algorithm for nonconvex optimization with convex constraints and its function-evaluation complexity, IMA Journal of Numerical Analysis, 32 (2012), pp. 1662–1695.
[8] C. Cartis and L. Roberts, A derivative-free Gauss-Newton method, Mathematical Programming Computation, 11 (2019), pp. 631–674.
[9] P. Conejo, E. Karas, and L. Pedroso, A trust-region derivative-free algorithm for constrained optimization, Optimization Methods and Software, 30 (2015), pp. 1126–1145.
[10] P. Conejo, E. Karas, L. Pedroso, A. Ribeiro, and M. Sachine, Global convergence of trust-region algorithms for convex constrained minimization without derivatives, Applied Mathematics and Computation, 220 (2013), pp. 324–330.
[11] A. Conn, K. Scheinberg, and Ph. Toint, A derivative free optimization algorithm in practice, in 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, USA, 1998, American Institute of Aeronautics and Astronautics.
[12] A. R. Conn, N. I. M. Gould, and P. L. Toint, Trust-Region Methods, vol. 1 of MPS-SIAM Series on Optimization, MPS/SIAM, Philadelphia, 2000.
[13] A. R. Conn, K. Scheinberg, and L. N. Vicente, Geometry of interpolation sets in derivative free optimization, Mathematical Programming, 111 (2007), pp. 141–172.
[14] , Geometry of sample sets in derivative-free optimization: Polynomial regression and underdetermined interpolation, IMA Journal of Numerical Analysis, 28 (2008), pp. 721–748.
[15] , Global convergence of general derivative-free trust-region algorithms to first- and second-order critical points, SIAM Journal on Optimization, 20 (2009), pp. 387–415.
[16] , Introduction to Derivative-Free Optimization, vol. 8 of MPS-SIAM Series on Optimization, MPS/SIAM, Philadelphia, 2009.
[17] R. Garmanjani, D. Júdice, and L. N. Vicente, Trust-region methods without using derivatives: Worst case complexity and the nonsmooth case, SIAM Journal on Optimization, 26 (2016), pp. 1987–2011.
[18] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 3rd ed., 1996.
[19] S. Gratton, P. L. Toint, and A. Tröltzsch, An active-set trust-region method for derivative-free nonlinear bound-constrained optimization, Optimization Methods and Software, 26 (2011), pp. 873–894.
[20] E. A. E. Gumma, M. H. A. Hashim, and M. M. Ali, A derivative-free algorithm for linearly constrained optimization problems, Computational Optimization and Applications, 57 (2014), pp. 599–621.
[21] D. A. Harville, Matrix Algebra From a Statistician’s Perspective, Springer, New York, 2008.
[22] M. Hough and L. Roberts, Model-based derivative-free methods for convex-constrained optimization, SIAM Journal on Optimization, 32 (2022), pp. 2552–2579.
[23] K. Kulshreshtha, B. Jurgelucks, F. Bause, J. Rautenberg, and C. Unverzagt, Increasing the sensitivity of electrical impedance to piezoelectric material parameters with non-uniform electrical excitation, Journal of Sensors and Sensor Systems, 4 (2015), pp. 217–227.
[24] J. Larson and S. C. Billups, Stochastic derivative-free optimization using a trust region framework, Computational Optimization and Applications, 64 (2016), pp. 619–645.
[25] J. Larson, M. Menickelly, and S. M. Wild, Derivative-free optimization methods, Acta Numerica, 28 (2019), pp. 287–404.
[26] M. J. D. Powell, Least Frobenius norm updating of quadratic models that satisfy interpolation conditions, Mathematical Programming, 100 (2004), pp. 183–215.
[27] , On the use of quadratic models in unconstrained minimization without derivatives, Optimization Methods and Software, 19 (2004), pp. 399–411.
[28] , The NEWUOA software for unconstrained optimization without derivatives, in Large-Scale Nonlinear Optimization, P. Pardalos, G. Di Pillo, and M. Roma, eds., vol. 83, Springer US, Boston, MA, 2006, pp. 255–297.
[29] , The BOBYQA algorithm for bound constrained optimization without derivatives, Tech. Rep. DAMTP 2009/NA06, University of Cambridge, 2009.
[30] M. Scheuerer and T. M. Hamill, Statistical postprocessing of ensemble precipitation forecasts by fitting censored, shifted gamma distributions, Monthly Weather Review, 143 (2015), pp. 4578 – 4596.
[31] S. F. B. Tett, J. M. Gregory, N. Freychet, C. Cartis, M. J. Mineter, and L. Roberts, Does model calibration reduce uncertainty in climate projections?, Journal of Climate, 35 (2022), pp. 2585–2602.
[32] L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, Philadelphia, 2020.
[33] G. Ughi, V. Abrol, and J. Tanner, An empirical study of derivative-free-optimization algorithms for targeted black-box attacks in deep neural networks, Optimization and Engineering, 23 (2022), pp. 1319–1346.
[34] S. M. Wild, Derivative-Free Optimization Algorithms for Computationally Expensive Functions, PhD thesis, Cornell University, 2009.
[35] S. M. Wild, R. G. Regis, and C. A. Shoemaker, ORBIT: Optimization by radial basis function interpolation in trust-regions, SIAM Journal on Scientific Computing, 30 (2008), pp. 3197–3219.
[36] M. Q. Xuan and J. Nocedal, A feasible method for constrained derivative-free optimization, arXiv preprint arXiv:2402.11920, (2024).
[37] Z. Zhang, Sobolev seminorm of quadratic functions with applications to derivative-free optimization, Mathematical Programming, 146 (2014), pp. 77–96.

Model Construction for Convex-Constrained Derivative-Free Optimization

Abstract

1 Introduction

1.1 Existing work

1.2 Contributions

Organization

Notation

2 Background

Assumption 2.1.

Assumption 2.2.

2.1 Model Construction & Accuracy

Assumption 2.3.

Definition 2.4.

2.2 Algorithm Specification

Assumption 2.5.

Theorem 2.6 (Theorem 3.10 & Corollary 3.15, [22]).

3 Linear Regression Models

Lemma 3.1.

Definition 3.2.

Lemma 3.3.

Proof.

Theorem 3.4.

Proof.

4 Underdetermined Quadratic Interpolation Models

Definition 4.1.

Lemma 4.2.

Proof.

Lemma 4.3.

Proof.

Lemma 4.4.

Proof.

Remark 4.5.

Lemma 4.6.

Proof.

Theorem 4.7.

Proof.

Remark 4.8.

5 Constructing 𝚲\bm{\Lambda}-Poised Quadratic Models

5.1 Updating the matrix 𝑭\bm{F}

Lemma 5.1.

Proof.

Theorem 5.2.

Proof.

Remark 5.3.

Remark 5.4.

5.2 Constructing 𝚲\bm{\Lambda}-Poised Sets

Theorem 5.5.

Proof.

Theorem 5.6.

Proof.

Remark 5.7.

Corollary 5.8.

6 Conclusions and Future Work

References

5 Constructing $\bm{\Lambda}$ -Poised Quadratic Models

5.1 Updating the matrix $\bm{F}$

5.2 Constructing $\bm{\Lambda}$ -Poised Sets