# Scalable Subspace Methods for Derivative-Free Nonlinear Least-Squares
Optimization

Coralia Cartis Mathematical Institute, University of Oxford, Radcliffe
Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, United Kingdom
(cartis@maths.ox.ac.uk). This work was supported by the EPSRC Centre for
Doctoral Training in Industrially Focused Mathematical Modelling
(EP/L015803/1) in collaboration with the Numerical Algorithms Group Ltd.
Lindon Roberts Mathematical Sciences Institute, Building 145, Science Road,
Australian National University, Canberra ACT 2601, Australia
(lindon.roberts@anu.edu.au). We note that the work in Sections 4 and 5
originally appeared in this author’s thesis [65, Chapter 7].

(27th August 2024)

###### Abstract

We introduce a general framework for large-scale model-based derivative-free
optimization based on iterative minimization within random subspaces. We
present a probabilistic worst-case complexity analysis for our method, where
in particular we prove high-probability bounds on the number of iterations
before a given optimality is achieved. This framework is specialized to
nonlinear least-squares problems, with a model-based framework based on the
Gauss-Newton method. This method achieves scalability by constructing local
linear interpolation models to approximate the Jacobian, and computes new
steps at each iteration in a subspace with user-determined dimension. We then
describe a practical implementation of this framework, which we call DFBGN. We
outline efficient techniques for selecting the interpolation points and search
subspace, yielding an implementation that has a low per-iteration linear
algebra cost (linear in the problem dimension) while also achieving fast
objective decrease as measured by evaluations. Extensive numerical results
demonstrate that DFBGN has improved scalability, yielding strong performance
on large-scale nonlinear least-squares problems.

Keywords: derivative-free optimization, large-scale optimization, nonlinear
least-squares, worst case complexity.

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

##  1 Introduction

An important class of nonlinear optimization methods is so-called derivative-
free optimization (DFO). In DFO, we consider problems where derivatives of the
objective (and/or constraints) are not available to be evaluated, and we only
have access to function values. This topic has received growing attention in
recent years, and is primarily used for objectives which are black-box (so
analytic derivatives or algorithmic differentiation are not available), and
expensive to evaluate or noisy (so finite differencing is impractical or
inaccurate). There are many types of DFO methods, such as model-based, direct
and pattern search, implicit filtering and others (see [50] for a recent
survey), and these techniques have been used in a variety of applications [1].

Here, we consider model-based DFO methods for unconstrained optimization,
which are based on iteratively constructing and minimizing interpolation
models for the objective. We also specialize these methods for nonlinear
least-squares problems, by constructing interpolation models for each residual
term rather than for the full objective [79, 73, 18].

This paper aims to provide a method that attempts to answer a key question
regarding model-based DFO: how to improve the scalability of this class.
Existing model-based DFO techniques are primarily designed for small- to
medium-scale problems, as the linear algebra cost of each iteration—largely
due to the cost of constructing interpolation models—means that their runtime
increases rapidly for large problems. There are several settings where
scalable DFO algorithms may be useful, such as data assimilation [9, 3],
machine learning [67, 35], generating adversarial examples for deep neural
networks [2, 68], image analysis [30], and as a possible proxy for global
optimization methods [19].

To address this, we introduce RSDFO, a scalable algorithmic framework for
model-based DFO. At each iteration of RSDFO we select a random low-dimensional
subspace, build and minimize a model to compute a step in this space, then
change the subspace at the next iteration. We provide a probabilistic worst-
case complexity analysis of RSDFO. To our knowledge, this is the first
_subspace model-based DFO method with global complexity and convergence
guarantees_. We then describe how this general framework can be specialized to
the case of nonlinear least-squares minimization through a model construction
technique inspired by the Gauss-Newton method, yielding a new algorithm RSDFO-
GN with associated worst-case complexity bounds. We then present an efficient
implementation of RSDFO-GN, which we call DFBGN. DFBGN is available on
Github111https://github.com/numericalalgorithmsgroup/dfbgn and includes
several algorithmic features that yield strong performance on large-scale
problems and a low per-iteration linear algebra cost that is typically linear
in the problem dimension.

###  1.1 Existing Literature

The contributions in this paper are connected to several areas of research. We
briefly review these topics below.

##### Block Coordinate Descent

There is a large body of work on (derivative-based) block coordinate descent
(BCD) methods, typically motivated by machine learning applications. BCD
extends coordinate search methods [75] by updating a subset of the variables
at each iteration, typically using a coordinate-wise variant of a first-order
method. For nonconvex problems, the first convergence result for a randomized
coordinate descent method based on proximal gradient descent was given in
[58]. Here, the sampling of coordinates was uniform and required step sizes
based on Lipschitz constants associated with the objective. This was extended
in [52] to general randomized block selection with a nonomonotone linesearch-
type method (to allow for unknown Lipschitz constants), and to a (possibly
deterministic) ‘essentially cyclic’ block selection and extrapolation (but
requiring Lipschitz constants) in [77]. Several extensions of this approach
have been developed, including the use of stochastic gradients [76], parallel
block updating [32] and inexact step calculations [32, 78].

BCD methods have been extended to nonlinear least-squares problems, leading to
so-called Subspace Gauss-Newton methods. These are derivative-based methods
where a Gauss-Newton step is computed for a subset of variables. This approach
was initially proposed in [70] for parameter estimation in climate
models—where derivative estimates were computed using implicit filtering
[49]—and analyzed in quadratic regularization and trust-region settings for
general unconstrained objectives in [15, 16].

##### Sketching

Sketching is an alternative dimensionality reduction technique for least-
squares problems, reducing the number of residuals rather than the number of
variables. Sketching ideas have been applied to linear [46, 53, 74] and
nonlinear [31] least-squares problems, as well as model-based DFO for
nonlinear least-squares [13], as well as subsampling algorithms for finite
sum-of-functions minimization such as Newton’s method [66, 6].

There are also alternative approaches to sketching which lead to subspace-type
methods, where local gradients and Hessians are estimated only within a
subspace (possibly used in conjunction with random subsampling). Sketching in
this context has been applied to, for example, Newton’s method [59, 40, 6],
BFGS [39], and SAGA [41], as well as to trust-region and quadratic
regularization methods [15, 16].

##### Random Embeddings for Global Optimization

Some global optimization methods have been proposed which randomly project a
high-dimensional problem into a low-dimensional subspace and solve this
smaller problem using existing (global or local) methods. Though applicable to
general global optimization problems (as a more sophisticated variant of
random search), this technique has been explored particularly for defeating
the curse of dimensionality when optimising functions which have low effective
dimensionality [64, 72, 17]. For the latter class, often only one random
subspace projection is needed, though the addition of constraints leads to
multiple embeddings being required [17]. Our approach here differs from these
works in both theoretical and numerical aspects, as it is focused on a
specific random subspace technique for local optimization.

##### Probabilistic Model-Based DFO

For model-based DFO, several algorithms have been developed and analyzed where
the local model at each iteration is only sufficiently accurate with a certain
probability [5, 21, 11]. Similar analysis also exists for derivative-based
algorithms [20, 43]. Our approach is based on deterministic model-based DFO
within subspaces, and we instead require a very weak probabilistic condition
on the (randomly chosen) subspaces (Assumption 2.9).

##### Randomized Direct Search DFO

In randomized direct search methods, iterates are perturbed in a random subset
of directions (rather than a positive spanning set) when searching for local
improvement. In this framework, effectively only a random subspace is searched
in each iteration. Worst-case complexity bounds for this technique are given
under predetermined step length regimes in [10, 36], and with adaptive step
sizes in [42, 44], where [44] extends [42] to linearly constrained problems.

##### Large-Scale DFO

There have been several alternative approaches considered for improving the
scalability of DFO. These often consider problems with specific structure
which enable efficient model construction, such as partial separability [24,
60], sparse Hessians [4], and minimization over the convex hull of finitely
many points [27]. On the other hand, there is a growing body of literature on
‘gradient sampling’ techniques for machine learning problems. These methods
typically consider stochastic first-order methods but with a gradient
approximation based on finite differencing in random directions [56]. This
framework has lead to variants of methods such as stochastic gradient descent
[34], SVRG [51] and Adam [22], for example. We note that linear interpolation
to orthogonal directions—more similar to traditional model-based DFO—has been
shown to outperform gradient sampling as a gradient estimation technique [8,
7].

##### Subspace DFO Methods

A model-based DFO method with similarities to our subspace approach is the
moving ridge function method from [45]. Here, existing objective evaluations
are used to determine an ‘active subspace’ which captures the largest
variability in the objective and build an interpolation model within this
subspace. We also note the VXQR method from [57], which performs line searches
along a direction chosen from a subspace determined by previous iterates. Both
of these methods do not include convergence theory. By comparison, aside from
our focus on nonlinear least-squares problems, both our general theoretical
framework and our implemented method select their working subspaces randomly,
and we provide (probabilistic) convergence guarantees.

###  1.2 Contributions

We introduce RSDFO (Randomized Subspace Derivative-Free Optimization), a
generic model-based DFO framework that relies on constructing a model in a
subspace at each iteration. Our novel approach enables model-based DFO methods
to be applied in a large-scale regime by giving the user explicit control over
the subspace dimension, and hence control over the per-iteration linear
algebra cost of the method. This framework is then specialized to the case of
nonlinear least-squares problems, yielding a new algorithm RSDFO-GN
(Randomized Subspace DFO with Gauss-Newton). The subspace model construction
framework of RSDFO-GN is based on DFO Gauss-Newton methods [18, 14], and
retains the same theoretical guarantees as RSDFO. We then describe a practical
implementation of RSDFO-GN, which we call DFBGN (Derivative-Free Block Gauss-
Newton).222Technically, DFBGN is not a block method as its subspaces are not
coordinate-aligned, but has already been released with this name. Compared to
existing methods, DFBGN reduces the linear algebra cost of model construction
and the initial objective evaluation cost by allowing fewer interpolation
points at every iteration. In order for DFBGN to have both scalability and a
similar evaluation efficiency to existing methods (i.e. objective reduction
achieved for a given number of objective evaluations), several modifications
to the theoretical framework, regarding the selection of interpolation points
and the search subspace, are necessary.

##### Theoretical Results

We consider a generic theoretical framework RSDFO, where the subspace
dimension is a user-chosen algorithm hyperparameter, and no specific model
construction approach is specified. Our framework is not specific to a least-
squares problem structure, and holds for any objective with Lipschitz
continuous gradient, and allows for a general class of random subspace
constructions (not relying on a specific class of embeddings or projections).
The theoretical results here extend the approach and techniques in [15, 16] to
model-based DFO methods. In particular, we use the notion of a well-aligned
subspace (Definition 2.8) from [15, 16], one in which sufficient decrease is
achievable, and assume that our search subspace is well-aligned with some
probability (Assumption 2.9). This is achieved provided we select a
sufficiently large subspace dimension (depending on the desired failure
probability and subspace alignment quality).

We derive a high probability worst-case complexity bound for RSDFO.
Specifically, our main bounds are of the form $\mathbb{P}\left[\min_{j\leq
k}\|\nabla f(\bm{x}_{j})\|\leq Ck^{-1/2}\right]\geq 1-e^{-ck}$ and
$\mathbb{P}\left[K_{\epsilon}\leq C\epsilon^{-2}\right]\leq
1-e^{-c\epsilon^{-2}}$, where $K_{\epsilon}$ is the first iteration to achieve
first-order optimality $\epsilon$ (see Theorem 2.18 and Corollary 2.19). This
result then implies a variety of alternative convergence results, such as
expectation bounds and almost-sure convergence. Based on [15, 16], we give
several constructions for determining our random subspace, and show that for
constructions based on Johnson-Lindenstrauss transformations, that we can
achieve convergence with a subspace dimension that is _independent of the
ambient dimension_.

Our analysis matches the standard deterministic $\mathcal{O}(\epsilon^{-2})$
complexity bounds for model-based DFO methods (e.g. [33, 18]). Compared to the
analysis of derivative-based methods (e.g. BCD [77] and probabilistically
accurate models [20]) we need to incorporate the possibility that the
interpolation model is not accurate (not fully linear, see Definition 2.1).
However, unlike [5, 21, 11] we do not assume that full linearity is a
stochastic property; instead, our stochasticity comes from the subspace
selection and we explicitly handle non-fully linear models similar to [26,
18]. This gives us a framework which is similar to standard model-based DFO
and with weak probabilistic conditions. Compared to the analysis of
derivative-based random subspace methods in [15, 16], our analysis is
complicated substantially by the possibility of inaccurate models and the
intricacies of model-based DFO algorithms.

We then consider RSDFO-GN, which explicitly describes how interpolation models
can be constructed for nonlinear least-squares problems, thus providing a
concrete implementation of RSDFO in this context. We prove that RSDFO-GN
retains the same complexity bounds as RSDFO.

##### Implementation

Although RSDFO-GN gives a reduced linear algebra cost of model construction
compared to existing methods, it is important that we have an implementation
that achieves this reduced cost without sacrificing overall performance (in
the sense of objective decrease achieved within a given number of objective
evaluations).

We introduce a practical, implementable variant of RSDFO-GN called DFBGN,
which is based on the solver DFO-LS [14]. DFBGN includes an efficient,
geoemtry-aware approach for selecting interpolation points, and hence
directions in our subspace, for removal, an adaptive randomized approach for
selecting new interpolation points/subspace directions. We study the per-
iteration linear algebra cost of DFBGN, and show that it is _linear in the
problem dimension_ , a substantial improvement over existing methods, which
are cubic in the problem dimension. Our per-iteration linear algebra costs are
also linear in the number of residuals, the same as existing methods, but with
a substantially smaller constant (quadratic in the subspace dimension, which
is user-determined, rather than quadratic in the problem dimension).

##### Numerical Results

We compare DFBGN with DFO-LS (which itself is shown to have state-of-the-art
performance in [14]) on collections of both medium-scale (approx. 100
dimensions) and large-scale test problems (approx. 1000 dimensions). We show
that DFBGN with a full-sized subspace has similar performance to DFO-LS in
terms of objective evaluations, but shows improved performance on runtime. As
the dimension of the subspace reduces (i.e. the size of the interpolation set
reduces), we demonstrate a tradeoff between reduced linear algebra costs and
increased evaluation counts required to achieve a given objective reduction.
The flexibility of DFBGN allows this tradeoff to be explicitly managed. When
tested on large-scale problems, DFO-LS frequently reaches a reasonable runtime
limit without making substantial progress, whereas DFBGN with small subspace
size can perform many more iterations and hence make better progress than DFO-
LS. In the case of expensive objectives with small evaluation budgets, we show
that DFBGN can make progress with few objective evaluations in a similar way
to DFO-LS (which has a mechanism to make progress from as few as 2 objective
evaluations independent of problem dimension), but with substantially lower
linear algebra costs.

##### Structure of paper

In Section 2 we describe RSDFO and provide our probabilistic worst-case
complexity analysis. We specialize RSDFO to RSDFO-GN in Section 3. Then we
describe the practical implementation DFBGN and its features in Section 4. Our
numerical results are given in Section 5.

##### Implementation

A Python implementation of DFBGN is available on
Github.333https://github.com/numericalalgorithmsgroup/dfbgn

##### Notation

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

##  2 Random Subspace Model-Based DFO

In this section we outline our general model-based DFO algorithmic framework
based on minimization in random subspaces. We consider the nonconvex problem

$\displaystyle\min_{\bm{x}\in\mathbb{R}^{n}}f(\bm{x}),$ (2.1)

where we assume that $f:\mathbb{R}^{n}\to\mathbb{R}$ is continuously
differentiable, but that access to its gradient is not possible (e.g. for the
reasons described in Section 1). In a standard model-based DFO framework (e.g.
[26, 50]), at each iteration $k$ we construct a quadratic model
$m_{k}:\mathbb{R}^{n}\to\mathbb{R}$ which approximates $f$ near our iterate
$\bm{x}_{k}$:

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

for some $\bm{g}_{k}\in\mathbb{R}^{n}$ and $H_{k}\in\mathbb{R}^{n\times n}$
symmetric. Based on this model, we build a globally convergent algorithm using
a trust-region framework [25]. This algorithmic framework is suitable
providing that—when necessary—we can guarantee $m_{k}$ is a sufficiently
accurate model for $f$ near $\bm{x}_{k}$. Details about how to construct
sufficiently accurate models based on interpolation are given in [26].

Our core idea here is to construct interpolation models which only approximate
the objective in a subspace, rather than in the full space $\mathbb{R}^{n}$.
This allows us to use interpolation sets with fewer points, since we do not
have to capture the objective’s behaviour outside our subspace, which improves
the scalability of the method.

In this section, we outline our general algorithmic framework and provide a
worst-case complexity analysis showing convergence to first-order stationary
points with high probability. We then describe how this framework may be
specialized to the case of nonlinear least-squares minimization.

###  2.1 RSDFO Algorithm

In our general framework, which we call RSDFO (Randomized Subspace DFO), we
modify the above approach by randomly choosing a $p$-dimensional subspace
(where $p<n$ is user-chosen) and constructing an interpolation model defined
only in that subspace.444Formally, we define our model in an affine space, but
we call it a subspace throughout as this fits with an intuitive view of what
RSDFO aims to achieve. Specifically, in each iteration $k$ we randomly choose
a $p$-dimensional affine space $\mathcal{Y}_{k}\subset\mathbb{R}^{n}$ given by
the range of $Q_{k}\in\mathbb{R}^{n\times p}$, i.e.

$\mathcal{Y}_{k}=\\{\bm{x}_{k}+Q_{k}\hat{\bm{s}}:\hat{\bm{s}}\in\mathbb{R}^{p}\\}.$
(2.3)

We then construct a model which interpolate $f$ at points in $\mathcal{Y}_{k}$
and ultimately construct a local quadratic model for $f$ only on
$\mathcal{Y}_{k}$. That is, given $Q_{k}$, we assume that we have
$\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ given by

$f(\bm{x}_{k}+Q_{k}\hat{\bm{s}})\approx\hat{m}_{k}(\hat{\bm{s}}):=f(\bm{x}_{k})+\hat{\bm{g}}_{k}^{T}\hat{\bm{s}}+\frac{1}{2}\hat{\bm{s}}^{T}\hat{H}_{k}\hat{\bm{s}},$
(2.4)

where $\hat{\bm{g}}_{k}\in\mathbb{R}^{p}$ and
$\hat{H}_{k}\in\mathbb{R}^{p\times p}$ are the low-dimensional model gradient
and Hessian respectively, adopting the convention of using hats on variables
to denote low-dimensional quantities. In Section 3 we specialize this to a
model construction process for nonlinear least-squares problems.

For our trust-region algorithm, we (approximately) minimize $\hat{m}_{k}$
inside the trust region to get a tentative step

$\hat{\bm{s}}_{k}\approx\operatorname*{arg\,min}_{\hat{\bm{s}}\in\mathbb{R}^{p}}\hat{m}_{k}(\hat{\bm{s}}),\quad\text{s.t.}\quad\|\hat{\bm{s}}\|\leq\Delta_{k},$
(2.5)

for the current trust-region radius $\Delta_{k}>0$, yielding a tentative step
$\bm{s}_{k}=Q_{k}\hat{\bm{s}}_{k}\in\mathbb{R}^{n}$. We thus also get the
computational advantage coming from solving a $p$-dimensional trust-region
subproblem.

In our setting we are only interested in the approximation properties of
$\hat{m}_{k}$ in the space $\mathcal{Y}_{k}$, and so we introduce the
following notion of a “sufficiently accurate” model:

######  Definition 2.1.

Given $Q\in\mathbb{R}^{n\times p}$, a model
$\hat{m}:\mathbb{R}^{p}\to\mathbb{R}$ is $Q$-fully linear in
$B(\bm{x},\Delta)\subset\mathbb{R}^{n}$ if

$\displaystyle|f(\bm{x}+Q\hat{\bm{s}})-\hat{m}(\hat{\bm{s}})|$
$\displaystyle\leq\kappa_{\rm ef}\Delta^{2},$ (2.6a)
$\displaystyle\|Q^{T}\nabla
f(\bm{x}+Q\hat{\bm{s}})-\nabla\hat{m}(\hat{\bm{s}})\|$
$\displaystyle\leq\kappa_{\rm eg}\Delta,$ (2.6b)

for all $\bm{s}\in\mathbb{R}^{p}$ with $\|\hat{\bm{s}}\|\leq\Delta$. The
constants $\kappa_{\rm ef}$ and $\kappa_{\rm eg}$ must be independent of $Q$,
$\hat{m}$, $\bm{x}$ and $\Delta$.

The gradient condition (2.6b) comes from noting that if
$\hat{f}(\hat{\bm{s}}):=f(\bm{x}+Q\hat{\bm{s}})$ then
$\nabla\hat{f}(\hat{\bm{s}})=Q^{T}\nabla f(\bm{x}+Q\hat{\bm{s}})$. We note
that if we have full-dimensional subspaces $p=n$ and take $Q=I$, then we
recover the standard notion of fully linear models [26, Definition 6.1].

1:Starting point $\bm{x}_{0}\in\mathbb{R}^{n}$, initial trust region radius
$\Delta_{0}>0$, and subspace dimension $p\in\\{1,\ldots,n\\}$.

2:Parameters: maximum trust-region radius $\Delta_{\rm max}\geq\Delta_{0}$,
trust-region radius scalings $0<\gamma_{\rm dec}<1<\gamma_{\rm
inc}\leq\overline{\gamma}_{\rm inc}$, criticality constants
$\epsilon_{C},\mu>0$ and trust-region scaling $0<\gamma_{C}<1$, safety step
threshold $\beta_{F}>0$ and trust-region scaling $0<\gamma_{F}<1$, and
acceptance thresholds $0<\eta_{1}\leq\eta_{2}<1$.

3:Set flag CHECK_MODEL=FALSE.

4:for $k=0,1,2,\ldots$ do

5: if CHECK_MODEL=TRUE then

6: Set $Q_{k}=Q_{k-1}$.

7: Construct a reduced model $\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ (2.4)
which is $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$.

8: else

9: Define a subspace by randomly sampling $Q_{k}\in\mathbb{R}^{n\times p}$.

10: Construct a reduced model $\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ (2.4)
which need not be $Q_{k}$-fully linear.

11: end if

12: if $\|\hat{\bm{g}}_{k}\|<\epsilon_{C}$ and
($\|\hat{\bm{g}}_{k}\|<\mu^{-1}\Delta_{k}$ or $\hat{m}_{k}$ is not
$Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$) then

13: (Criticality step) Set $\bm{x}_{k+1}=\bm{x}_{k}$,
$\Delta_{k+1}=\gamma_{C}\Delta_{k}$ and CHECK_MODEL=TRUE.

14: else $\leftarrow$ $\|\hat{\bm{g}}_{k}\|\geq\epsilon_{C}$ or
($\|\hat{\bm{g}}_{k}\|\geq\mu^{-1}\Delta_{k}$ and $\hat{m}_{k}$ is
$Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$)

15: Approximately solve the subspace trust-region subproblem in
$\mathbb{R}^{p}$ (2.5) and calculate the step
$\bm{s}_{k}=Q_{k}\hat{\bm{s}}_{k}\in\mathbb{R}^{n}$.

16: if $\|\hat{\bm{s}}_{k}\|<\beta_{F}\Delta_{k}$ then

17: (Safety step) Set $\bm{x}_{k+1}=\bm{x}_{k}$ and
$\Delta_{k+1}=\gamma_{F}\Delta_{k}$.

18: If $\hat{m}_{k}$ is $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$,
then set CHECK_MODEL=FALSE, otherwise CHECK_MODEL=TRUE.

19: else

20: Evaluate $f(\bm{x}_{k}+\bm{s}_{k})$ and calculate ratio
$\rho_{k}:=\frac{f(\bm{x}_{k})-f(\bm{x}_{k}+\bm{s}_{k})}{\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})}.$
(2.7)

21: Accept/reject step and update trust region radius: set
$\bm{x}_{k+1}=\begin{cases}\bm{x}_{k}+\bm{s}_{k},&\rho_{k}\geq\eta_{1},\\\
\bm{x}_{k},&\rho_{k}<\eta_{1},\end{cases}\quad\text{and}\quad\Delta_{k+1}=\begin{cases}\min(\max(\gamma_{\rm
inc}\Delta_{k},\overline{\gamma}_{\rm inc}\|\hat{\bm{s}}_{k}\|),\Delta_{\rm
max}),&\rho_{k}\geq\eta_{2},\\\ \max(\gamma_{\rm
dec}\Delta_{k},\|\hat{\bm{s}}_{k}\|),&\eta_{1}\leq\rho_{k}<\eta_{2},\\\
\min(\gamma_{\rm
dec}\Delta_{k},\|\hat{\bm{s}}_{k}\|),&\rho_{k}<\eta_{1}.\end{cases}$ (2.8)

22: If $\rho_{k}\geq\eta_{2}$ or $\hat{m}_{k}$ is $Q_{k}$-fully linear in
$B(\bm{x}_{k},\Delta_{k})$, then set CHECK_MODEL=FALSE, otherwise set
CHECK_MODEL=TRUE.

23: end if

24: end if

25:end for

Algorithm 1 RSDFO (Randomized Subspace Derivative-Free Optimization) for
solving (2.1).

##### Complete RSDFO Algorithm

The complete RSDFO algorithm is stated in Algorithm 1. The overall structure
is common to model-based DFO methods [26]. In particular, we assume that we
have procedures to verify whether or not a model is $Q_{k}$-fully linear in
$B(\bm{x}_{k},\Delta_{k})$ and (if not) to generate a $Q_{k}$-fully linear
model. When we specialize RSDFO to nonlinear least-squares problems in Section
3, we will describe how we can obtain such procedures.

After defining our subspace and model, we first perform a criticality step,
which guarantees that—whenever we suspect we are close to first-order
stationarity, as measured by $\|\hat{\bm{g}}_{k}\|$—we have an accurate model
and an appropriately sized trust-region radius. Often this criticality step is
formulated as its own subroutine with an extra inner loop [26], but following
[33] we only perform one criticality step per iteration and avoid nested
loops.

The remainder of Algorithm 1 broadly follows standard trust-region methods. If
the computed step is much shorter than the trust-region radius, we enter a
safety step—originally due to Powell [61]—which is similar to an unsuccessful
iteration (i.e. $\rho_{k}<\eta_{1}$, so we reject the step and decrease
$\Delta_{k}$) but without evaluating $f(\bm{x}_{k}+\bm{s}_{k})$ and hence
saving one objective evaluation. Our updating mechanism for $\Delta_{k}$ (2.8)
also takes into account the computed step size, and is the same as in [63,
18].

An important feature of RSDFO is that in some iterations, we reuse the
previous subspace, $\mathcal{Y}_{k}=\mathcal{Y}_{k-1}$, corresponding to the
flag CHECK_MODEL=TRUE. In this case, we had an inaccurate model in iteration
$k-1$ and require that our new model $\hat{m}_{k}$ is accurate ($Q_{k}$-fully
linear). This mechanism essentially ensures that $\Delta_{k}$ is not decreased
too quickly as a result of inaccurate models, and is mostly decreased to
achieve sufficient objective reduction.

We now give our convergence and worst-case complexity analysis of Algorithm 1.

###  2.2 Assumptions and Preliminary Results

We begin our analysis with some basic assumptions and preliminary results.

######  Assumption 2.2 (Smoothness).

The objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ is bounded below by
$f_{\rm low}$ and continuously differentiable, and $\nabla f$ is $L_{\nabla
f}$-Lipschitz continuous in the extended level set
$\\{\bm{y}\in\mathbb{R}^{n}:\|\bm{y}-\bm{x}\|\leq\Delta_{\max}\>\text{for
some}\>f(\bm{x})\leq f(\bm{x}_{0})\\}$, for some constant $L_{\nabla f}>0$.

We also need two standard assumptions for trust-region methods: uniformly
bounded above model Hessians and sufficiently accurate solutions to the trust-
region subproblem (2.5).

######  Assumption 2.3 (Bounded model Hessians).

We assume that $\|\hat{H}_{k}\|\leq\kappa_{H}$ for all $k$, for some
$\kappa_{H}\geq 1$.

######  Assumption 2.4 (Cauchy decrease).

Our method for solving the trust-region subproblem (2.5) gives a step
$\hat{\bm{s}}_{k}$ satisfying the sufficient decrease condition

$\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})\geq
c_{1}\|\hat{\bm{g}}_{k}\|\min\left(\Delta_{k},\frac{\|\hat{\bm{g}}_{k}\|}{\max(\|\hat{H}_{k}\|,1)}\right),$
(2.9)

for some $c_{1}\in[1/2,1]$ independent of $k$.

A useful consequence, needed for the analysis of our trust-region radius
updating scheme, is the following.

######  Lemma 2.5 (Lemma 3.6, [18]).

Suppose Assumption 2.4 holds. Then

$\|\hat{\bm{s}}_{k}\|\geq
c_{2}\min\left(\Delta_{k},\frac{\|\hat{\bm{g}}_{k}\|}{\max(\|\hat{H}_{k}\|,1)}\right),$
(2.10)

where $c_{2}:=2c_{1}/(1+\sqrt{1+2c_{1}})$.

######  Lemma 2.6.

Suppose Assumptions 2.3 and 2.4 hold, and we run RSDFO with $\beta_{F}\leq
c_{2}$. If $\hat{m}_{k}$ is $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$
and

$\Delta_{k}\leq c_{0}\|\hat{\bm{g}}_{k}\|,\qquad\text{where}\qquad
c_{0}:=\min\left(\mu,\frac{1}{\kappa_{H}},\frac{c_{1}(1-\eta_{2})}{2\kappa_{\rm
ef}}\right),$ (2.11)

then the criticality and safety steps are not called, and
$\rho_{k}\geq\eta_{2}$.

###### Proof.

Since $\hat{m}_{k}$ is $Q_{k}$-fully linear and
$\Delta_{k}\leq\mu\|\hat{\bm{g}}_{k}\|$, the criticality step is not called.
From Lemma 2.5 and $\Delta_{k}\leq\|\hat{\bm{g}}_{k}\|/\kappa_{H}$, we have
$\|\hat{\bm{s}}_{k}\|\geq c_{2}\Delta_{k}\geq\beta_{F}\Delta_{k}$ and so the
safety step is not called.

From Assumptions 2.3 and 2.4, we have

$\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})\geq
c_{1}\|\hat{\bm{g}}_{k}\|\min\left(\Delta_{k},\frac{\|\hat{\bm{g}}_{k}\|}{\kappa_{H}}\right)=c_{1}\|\hat{\bm{g}}_{k}\|\Delta_{k},$
(2.12)

since $\Delta_{k}\leq\|\hat{\bm{g}}_{k}\|/\kappa_{H}$ by assumption. Next,
since $\hat{m}_{k}$ is $Q_{k}$-fully linear, from (2.6a) we have

$\displaystyle|f(\bm{x}_{k})-\hat{m}_{k}(\bm{0})|$
$\displaystyle=|f(\bm{x}_{k}+Q_{k}\bm{0})-\hat{m}_{k}(\bm{0})|\leq\kappa_{\rm
ef}\Delta_{k}^{2},$ (2.13)
$\displaystyle|f(\bm{x}_{k}+\bm{s}_{k})-\hat{m}_{k}(\hat{\bm{s}}_{k})|$
$\displaystyle=|f(\bm{x}_{k}+Q_{k}\hat{\bm{s}}_{k})-\hat{m}_{k}(\hat{\bm{s}}_{k})|\leq\kappa_{\rm
ef}\Delta_{k}^{2}.$ (2.14)

Hence we have

$\displaystyle|\rho_{k}-1|$
$\displaystyle\leq\frac{|f(\bm{x}_{k})-\hat{m}_{k}(\bm{0})|}{|\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})|}+\frac{|f(\bm{x}_{k}+\bm{s}_{k})-\hat{m}_{k}(\hat{\bm{s}}_{k})|}{|\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})|}\leq\frac{2\kappa_{\rm
ef}\Delta_{k}^{2}}{c_{1}\|\hat{\bm{g}}_{k}\|\Delta_{k}}\leq 1-\eta_{2},$
(2.15)

since $\Delta_{k}\leq c_{0}\|\hat{\bm{g}}_{k}\|\leq
c_{1}(1-\eta_{2})\|\hat{\bm{g}}_{k}\|/(2\kappa_{\rm ef})$. Thus
$\rho_{k}\geq\eta_{2}$, which is the claim of the lemma. ∎

######  Remark 2.7.

The requirement $\beta_{F}\leq c_{2}$ in Lemma 2.6 is not restrictive. Since
have $c_{1}\geq 1/2$ in Assumption 2.4, it suffices to choose
$\beta_{F}\leq\sqrt{2}-1$, for example.

Our key new assumption is on the quality of our subspace selection, as
suggested in [15, 16]:

######  Definition 2.8.

The matrix $Q_{k}$ is well-aligned if

$\|Q_{k}^{T}\nabla f(\bm{x}_{k})\|\geq\alpha_{Q}\|\nabla f(\bm{x}_{k})\|,$
(2.16)

for some $\alpha_{Q}\in(0,1)$ independent of $k$.

######  Assumption 2.9 (Subspace quality).

Our subspace selection (determined by $Q_{k}$) satisfies the following two
properties:

  1. (a)

At each iteration $k$ of RSDFO in which CHECK_MODEL = FALSE, our subspace
selection $Q_{k}$ is well-aligned for some fixed $\alpha_{Q}\in(0,1)$ with
probability at least $1-\delta_{S}$, for some $\delta_{S}\in(0,1)$,
independently of $\\{Q_{0},\ldots,Q_{k-1}\\}$.

  2. (b)

$\|Q_{k}\|\leq Q_{\max}$ for all $k$ and some $Q_{\max}>0$.

Of these two properties, (a) is needed for our complexity analysis, while (b)
is only needed in order to construct $Q_{k}$-fully linear models (in Section
3). We will discuss how to achieve Assumption 2.9 in more detail in Section
2.6.

######  Lemma 2.10.

In all iterations $k$ of RSDFO where the criticality step is not called, we
have $\|\hat{\bm{g}}_{k}\|\geq\min(\epsilon_{C},\mu^{-1}\Delta_{k})$. If the
criticality step is not called in iteration $k$, $Q_{k}$ is well-aligned and
$\|\nabla f(\bm{x}_{k})\|\geq\epsilon$, then

$\|\hat{\bm{g}}_{k}\|\geq\epsilon_{g}(\epsilon):=\min\left(\epsilon_{C},\frac{\alpha_{Q}\epsilon}{\kappa_{\rm
eg}\mu+1}\right)>0.$ (2.17)

###### Proof.

The first part follows immediately from the entry condition of the criticality
step. To prove (2.17), suppose the criticality step is not called in iteration
$k$ and $\|\hat{\bm{g}}_{k}\|<\epsilon_{C}$. Then we have
$\Delta_{k}\leq\mu\|\hat{\bm{g}}_{k}\|$ and $\hat{m}_{k}$ is $Q_{k}$-fully
linear, and so from (2.6b) we have

$\|Q_{k}^{T}\nabla f(\bm{x}_{k})\|\leq\|Q_{k}^{T}\nabla
f(\bm{x}_{k})-\hat{\bm{g}}_{k}\|+\|\hat{\bm{g}}_{k}\|\leq\kappa_{\rm
eg}\Delta_{k}+\|\hat{\bm{g}}_{k}\|\leq(\kappa_{\rm
eg}\mu+1)\|\hat{\bm{g}}_{k}\|.$ (2.18)

Since $Q_{k}$ is well-aligned, we conclude from (2.16) and (2.18) that

$\alpha_{Q}\|\nabla f(\bm{x}_{k})\|\leq\|Q_{k}^{T}\nabla
f(\bm{x}_{k})\|\leq(\kappa_{\rm eg}\mu+1)\|\hat{\bm{g}}_{k}\|,$ (2.19)

and we are done, since $\|\nabla f(\bm{x}_{k})\|\geq\epsilon$. ∎

###  2.3 Counting iterations

We now provide a series of results counting the number of iterations of RSDFO
of different types, following the style of analysis from [20, 16]. First we
introduce some notation to enumerate our iterations. Suppose we run RSDFO
until the end of iteration $K$. We then define the following subsets of
$\\{0,\ldots,K\\}$:

  * •

$\mathcal{C}$ is the set of iterations in $\\{0,\ldots,K\\}$ where the
criticality step is called.

  * •

$\mathcal{F}$ is the set of iterations in $\\{0,\ldots,K\\}$, where the safety
step is called (i.e. $\|\hat{\bm{s}}_{k}\|<\beta_{F}\Delta_{k}$).

  * •

$\mathcal{VS}$ is the set of very successful iterations in $\\{0,\ldots,K\\}$,
where $\rho_{k}\geq\eta_{2}$.

  * •

$\mathcal{S}$ is the set of successful iterations in $\\{0,\ldots,K\\}$, where
$\rho_{k}\geq\eta_{1}$. Note that $\mathcal{VS}\subset\mathcal{S}$.

  * •

$\mathcal{U}$ is the set of unsuccessful iterations in $\\{0,\ldots,K\\}$,
where $\rho_{k}<\eta_{1}$.

  * •

$\mathcal{A}$ is the set of well-aligned iterations in $\\{0,\ldots,K\\}$,
where (2.16) holds.

  * •

$\mathcal{A}^{C}$ is the set of poorly aligned iterations in
$\\{0,\ldots,K\\}$, where $\eqref{eq_sketched_gradient}$ does not hold.

  * •

$\mathcal{D}(\Delta)$ is the set of iterations in $\\{0,\ldots,K\\}$ where
$\Delta_{k}\geq\Delta$ for some $\Delta>0$.

  * •

$\mathcal{D}^{C}(\Delta)$ is the set of iterations in $\\{0,\ldots,K\\}$ where
$\Delta_{k}<\Delta$.

  * •

$\mathcal{L}$ is the set of iterations in $\\{0,\ldots,K\\}$ where
$\hat{m}_{k}$ is $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$.

  * •

$\mathcal{L}^{C}$ is the set of iterations in $\\{0,\ldots,K\\}$ where
$\hat{m}_{k}$ is not $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$.

In particular, we have the partitions, for any $\Delta>0$,

$\\{0,\ldots,K\\}=\mathcal{C}\cup\mathcal{F}\cup\mathcal{S}\cup\mathcal{U}=\mathcal{A}\cup\mathcal{A}^{C}=\mathcal{D}(\Delta)\cup\mathcal{D}^{C}(\Delta)=\mathcal{L}\cup\mathcal{L}^{C}.$
(2.20)

First, we bound the number of successful iterations with large $\Delta_{k}$
using standard arguments from trust-region methods. Throughout, we use
$\\#(\cdot)$ to refer to the cardinality of a set of iterations.

######  Lemma 2.11.

Suppose Assumptions 2.2, 2.3 and 2.4 hold. If $\|\nabla
f(\bm{x}_{k})\|\geq\epsilon$ for all $k=0,\ldots,K$, then

$\\#(\mathcal{A}\cap\mathcal{D}(\Delta)\cap\mathcal{S})\leq\phi(\Delta,\epsilon):=\frac{f(\bm{x}_{0})-f_{\rm
low}}{\eta_{1}c_{1}\epsilon_{g}(\epsilon)\min(\epsilon_{g}(\epsilon)/\kappa_{H},\Delta)},$
(2.21)

for all $\Delta>0$.

###### Proof.

Since $\|\nabla f(\bm{x}_{k})\|\geq\epsilon$, from Lemma 2.10 we have
$\|\hat{\bm{g}}_{k}\|\geq\epsilon_{g}(\epsilon)$ for all
$k\in\mathcal{A}\cap\mathcal{S}$ (noting that $k\in\mathcal{S}$ implies
$k\in\mathcal{C}^{C}$). Then since $\rho_{k}\geq\eta_{1}$, from Assumptions
2.3 and 2.4 we get

$\displaystyle f(\bm{x}_{k})-f(\bm{x}_{k+1})$
$\displaystyle\geq\eta_{1}[\hat{m}_{k}(\bm{0})-\hat{m}_{k}(\hat{\bm{s}}_{k})],$
(2.22)
$\displaystyle\geq\eta_{1}c_{1}\|\hat{\bm{g}}_{k}\|\min\left(\Delta_{k},\frac{\|\hat{\bm{g}}_{k}\|}{\kappa_{H}}\right),$
(2.23)
$\displaystyle\geq\eta_{1}c_{1}\epsilon_{g}(\epsilon)\min\left(\Delta,\frac{\epsilon_{g}(\epsilon)}{\kappa_{H}}\right),$
(2.24)

where the last line follows from
$\|\hat{\bm{g}}_{k}\|\geq\epsilon_{g}(\epsilon)$ and $\Delta_{k}\geq\Delta$
(from $k\in\mathcal{D}(\Delta)$). Since our step acceptance guarantees our
algorithm is monotone (i.e. $f(\bm{x}_{k+1})\leq f(\bm{x}_{k})$ for all $k$),
we get

$f(\bm{x}_{0})-f_{\rm
low}\geq\sum_{k\in\mathcal{A}\cap\mathcal{D}(\Delta)\cap\mathcal{S}}[f(\bm{x}_{k})-f(\bm{x}_{k+1})]\geq\left[\eta_{1}c_{1}\epsilon_{g}(\epsilon)\min\left(\Delta,\frac{\epsilon_{g}(\epsilon)}{\kappa_{H}}\right)\right]\cdot\\#(\mathcal{A}\cap\mathcal{D}(\Delta)\cap\mathcal{S}),$
(2.25)

from which the result follows. ∎

######  Lemma 2.12.

Suppose Assumptions 2.2, 2.3 and 2.4 hold, and $\beta_{F}\leq c_{2}$. If
$\|\nabla f(\bm{x}_{k})\|\geq\epsilon$ for all $k=0,\ldots,K$, then

$\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}\setminus\mathcal{VS})=0,$
(2.26)
$\Delta\leq\Delta^{*}(\epsilon):=\min\left(c_{0}\epsilon_{g}(\epsilon),\frac{\alpha_{Q}\epsilon}{\kappa_{\rm
eg}+\mu^{-1}}\right).$ (2.27)

###### Proof.

To find a contradiction, first suppose
$k\in\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}\cap\mathcal{C}^{C}\setminus\mathcal{VS}$.
Then since $k\in\mathcal{A}\cap\mathcal{C}^{C}$ and $\|\nabla
f(\bm{x}_{k})\|\geq\epsilon$ by assumption, we have
$\|\hat{\bm{g}}_{k}\|\geq\epsilon_{g}(\epsilon)$ from Lemma 2.10. Since
$k\in\mathcal{D}^{C}(\Delta)$, we have
$\Delta_{k}<\Delta\leq\Delta^{*}(\epsilon)\leq c_{0}\epsilon_{g}(\epsilon)\leq
c_{0}\|\hat{\bm{g}}_{k}\|$ by definition of $\Delta^{*}(\epsilon)$. From this
and $k\in\mathcal{L}$, the assumptions of Lemma 2.6 are met, so
$k\notin\mathcal{F}$ and $\rho_{k}\geq\eta_{2}$; that is, $k\in\mathcal{VS}$,
a contradiction. Hence we have
$\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}\cap\mathcal{C}^{C}\setminus\mathcal{VS})=0$.

Next, we suppose
$k\in\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}\cap\mathcal{C}$ and
again look for a contradiction. In this case, we have
$\Delta_{k}<\Delta\leq\Delta^{*}(\epsilon)\leq\alpha_{Q}\epsilon/(\kappa_{\rm
eg}+\mu^{-1})$, and so from $k\in\mathcal{A}\cap\mathcal{L}$ and $\|\nabla
f(\bm{x}_{k})\|\geq\epsilon$ we have

$\|\hat{\bm{g}}_{k}\|\geq\|Q_{k}^{T}\nabla f(\bm{x}_{k})\|-\|Q_{k}^{T}\nabla
f(\bm{x}_{k})-\hat{\bm{g}}_{k}\|\geq\alpha_{Q}\epsilon-\kappa_{\rm
eg}\Delta_{k}>\mu^{-1}\Delta_{k}.$ (2.28)

This means we have $\|\hat{\bm{g}}_{k}\|>\mu^{-1}\Delta_{k}$ and
$k\in\mathcal{L}$, so the criticality step is not entered; i.e.
$k\in\mathcal{C}^{C}$, a contradiction. Hence we have
$\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}\cap\mathcal{C})=0$
and we are done. ∎

######  Lemma 2.13.

Suppose Assumptions 2.2, 2.3 and 2.4 hold. Then we have

$\displaystyle\\#(\mathcal{D}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\Delta)\setminus\mathcal{S})$ $\displaystyle\leq
C_{1}\\#(\mathcal{D}(\overline{\gamma}_{\rm
inc}^{-1}\Delta)\cap\mathcal{S})+C_{2},$ (2.29)

for all $\Delta\leq\Delta_{0}$, where

$C_{1}:=\frac{\log(\overline{\gamma}_{\rm
inc})}{\log(1/\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec}))}\qquad\text{and}\qquad
C_{2}:=\frac{\log(\Delta_{0}/\Delta)}{\log(1/\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec}))}.$ (2.30)

###### Proof.

If $k\in\mathcal{S}$, we always have $\Delta_{k+1}\leq\overline{\gamma}_{\rm
inc}\Delta_{k}$. On the other hand, if $k\in\mathcal{U}$, we have

$\Delta_{k+1}=\min(\gamma_{\rm
dec}\Delta_{k},\|\hat{\bm{s}}_{k}\|)\leq\gamma_{\rm dec}\Delta_{k},$ (2.31)

Hence,

$\Delta_{k+1}\leq\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})\Delta_{k},\qquad\text{for all
$k\in\mathcal{C}\cup\mathcal{F}\cup\mathcal{U}$}.$ (2.32)

We now consider the value of $\log(\Delta_{k})$ for $k=0,\ldots,K$, so at each
iteration we have an additive change:

  * •

Since $\Delta\leq\Delta_{0}$, the threshold value $\log(\Delta)$ is
$\log(\Delta_{0}/\Delta)$ below the starting value $\log(\Delta_{0})$.

  * •

If $k\in\mathcal{S}$, then $\log(\Delta_{k})$ increases by at most
$\log(\overline{\gamma}_{\rm inc})$. In particular, $\Delta_{k+1}\geq\Delta$
is only possible if $\Delta_{k}\geq\overline{\gamma}_{\rm inc}^{-1}\Delta$.

  * •

If $k\notin\mathcal{S}=\mathcal{C}\cup\mathcal{F}\cup\mathcal{U}$, then
$\log(\Delta_{k})$ decreases by at least
$|\log(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec}))|=\log(1/\max(\gamma_{C},\gamma_{F},\gamma_{\rm dec}))$.

Now, any decrease in $\Delta_{k}$ coming from
$k\in\mathcal{D}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\Delta)\setminus\mathcal{S}$ yields $\Delta_{k+1}\geq\Delta$. Hence
the total decrease in $\log(\Delta_{k})$ must be fully matched by the initial
gap $\log(\Delta_{0}/\Delta)$ plus the maximum possible amount that
$\log(\Delta_{k})$ can be increased above $\log(\Delta)$. That is, we must
have

$\displaystyle\log(1/\max(\gamma_{C},\gamma_{F},\gamma_{\rm dec}))$
$\displaystyle\cdot\\#(\mathcal{D}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\Delta)\setminus\mathcal{S})$
$\displaystyle\qquad\leq\log(\Delta_{0}/\Delta)+\log(\overline{\gamma}_{\rm
inc})\cdot\\#(\mathcal{D}(\overline{\gamma}_{\rm
inc}^{-1}\Delta)\cap\mathcal{S}),$ (2.33)

which gives us (2.29). ∎

######  Lemma 2.14.

Suppose Assumptions 2.2, 2.3 and 2.4 hold. Then

$\\#(\mathcal{D}^{C}(\gamma_{\rm inc}^{-1}\Delta)\cap\mathcal{VS})\leq
C_{3}\cdot\\#(\mathcal{D}^{C}(\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-1}\Delta)\setminus\mathcal{VS}),$ (2.34)

for all $\Delta\leq\min(\Delta_{0},\gamma_{\rm inc}^{-1}\Delta_{\max})$, where

$C_{3}:=\frac{\log(1/\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F}))}{\log(\gamma_{\rm inc})}.$ (2.35)

###### Proof.

We follow a similar reasoning to the proof of Lemma 2.13. For every iteration
$k\in\mathcal{VS}\cap\mathcal{D}^{C}(\Delta)$, we increase $\Delta_{k}$ by a
factor of at least $\gamma_{\rm inc}$, since $\Delta_{k}<\Delta\leq\gamma_{\rm
inc}^{-1}\Delta_{\max}$. Equivalently, we increase $\log(\Delta_{k})$ by at
least $\log(\gamma_{\rm inc})$. In particular, if $\Delta_{k}<\gamma_{\rm
inc}^{-1}\Delta$, then $\Delta_{k+1}<\Delta$.

Alternatively, if $k\in\mathcal{S}\setminus\mathcal{VS}$, we set

$\Delta_{k+1}=\max(\gamma_{\rm
dec}\Delta_{k},\|\hat{\bm{s}}_{k}\|)\geq\gamma_{\rm dec}\Delta_{k}.$ (2.36)

If $k\in\mathcal{U}$ we set

$\Delta_{k+1}=\min(\gamma_{\rm
dec}\Delta_{k},\|\hat{\bm{s}}_{k}\|)\geq\min(\gamma_{\rm
dec},\beta_{F})\Delta_{k},$ (2.37)

since $\|\hat{\bm{s}}_{k}\|\geq\beta_{F}\Delta_{k}$ from $k\notin\mathcal{F}$.
Hence, for every iteration $k\notin\mathcal{VS}$, we decrease $\Delta_{k}$ by
a factor of at most $\min(\gamma_{C},\gamma_{F},\gamma_{\rm dec},\beta_{F})$,
or equivalently we decrease $\log(\Delta_{k})$ by at most the amount
$|\log(\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F}))|=\log(1/\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F}))$. Then, to have $\Delta_{k+1}<\Delta$ we require
$\Delta_{k}<\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-1}\Delta$.

Therefore, since $\Delta_{0}\geq\Delta$, the total increase in
$\log(\Delta_{k})$ from $k\in\mathcal{VS}\cap\mathcal{D}^{C}(\gamma_{\rm
inc}^{-1}\Delta)$ must be fully matched by the total decrease in
$\log(\Delta_{k})$ from
$k\in\mathcal{D}^{C}(\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-1}\Delta)\setminus\mathcal{VS}$. That is,

$\displaystyle\log(\gamma_{\rm inc})$
$\displaystyle\\#(\mathcal{VS}\cap\mathcal{D}^{C}(\gamma_{\rm
inc}^{-1}\Delta))$
$\displaystyle\qquad\qquad\leq\log(1/\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F}))\\#(\mathcal{D}^{C}(\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-1}\Delta)\setminus\mathcal{VS}),$ (2.38)

and we are done. ∎

######  Lemma 2.15.

Suppose Assumptions 2.2, 2.3 and 2.4 hold. Then

$\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}^{C}\setminus\mathcal{VS})\leq\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L})+1,$
(2.39)

for all $\Delta>0$.

###### Proof.

After every iteration $k$ where $\hat{m}_{k}$ is not $Q_{k}$-fully linear and
either the criticality step is called or $\rho_{k}<\eta_{2}$, we always set
$\Delta_{k+1}\leq\Delta_{k}$ and CHECK_MODEL=TRUE. This means that
$Q_{k+1}=Q_{k}$, so $Q_{k+1}$ is well-aligned if and only if $Q_{k}$ is well-
aligned. Hence if
$k\in\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}^{C}\setminus\mathcal{VS}$
then either $k=K$ or
$k+1\in\mathcal{A}\cap\mathcal{D}^{C}(\Delta)\cap\mathcal{L}$, and we are
done. ∎

We are now in a position to bound the total number of well-aligned iterations.

######  Lemma 2.16.

Suppose Assumptions 2.2, 2.3 and 2.4 hold, and both $\beta_{F}\leq c_{2}$ and
$\gamma_{\rm inc}>\min(\gamma_{C},\gamma_{F},\gamma_{\rm dec},\beta_{F})^{-2}$
hold. Then if $\|\nabla f(\bm{x}_{k})\|\geq\epsilon$ for all $k=0,\ldots,K$,
we have

$\displaystyle\\#(\mathcal{A})$
$\displaystyle\leq\psi(\epsilon)+\frac{C_{4}}{1+C_{4}}(K+1),$ (2.40)

where

$\displaystyle\psi(\epsilon)$
$\displaystyle:=\frac{1}{1+C_{4}}\left[(C_{1}+2)\phi(\Delta_{\min}(\epsilon),\epsilon)+\frac{4\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min}(\epsilon),\epsilon)}{1-2C_{3}}\right.$ (2.41)
$\displaystyle\qquad\qquad\qquad\qquad\left.+C_{2}+\frac{2}{1-2C_{3}}+1\right],$
(2.42) $\displaystyle\Delta_{\min}(\epsilon)$
$\displaystyle:=\min\left(\overline{\gamma}_{\rm
inc}^{-1}\Delta_{0},\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})\overline{\gamma}_{\rm inc}^{-1}\Delta^{*}(\epsilon)\right),$ (2.43)
$\displaystyle C_{4}$
$\displaystyle:=\max\left(C_{1},\frac{4C_{3}}{1-2C_{3}}\right)>0.$ (2.44)

In these expressions, the values $C_{1}$ and $C_{2}$ are defined in Lemma
2.13, $C_{3}$ is defined in Lemma 2.14, $\phi(\cdot,\epsilon)$ is defined in
Lemma 2.11, and $\epsilon_{g}(\epsilon)$ and $\Delta^{*}(\epsilon)$ are
defined in Lemmas 2.10 and 2.12 respectively.

###### Proof.

For ease of notation, we will write $\Delta_{\min}$ in place of
$\Delta_{\min}(\epsilon)$. We begin by noting that $\gamma_{\rm
inc}>\min(\gamma_{C},\gamma_{F},\gamma_{\rm dec},\beta_{F})^{-2}$ implies that
$C_{3}\in(0,1/2)$, which we will use later.

Next, we have

$\displaystyle\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min}))$
$\displaystyle=\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\setminus\mathcal{S}),$
(2.45)
$\displaystyle\leq\phi(\Delta_{\min},\epsilon)+\\#(\mathcal{A}\cap\mathcal{D}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm inc}\Delta_{\min})\setminus\mathcal{S})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm inc}\Delta_{\min})\setminus\mathcal{S}),$
(2.46)
$\displaystyle\leq\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm inc}\Delta_{\min})\setminus\mathcal{S}),$
(2.47)
$\displaystyle=\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}\setminus\mathcal{S})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}^{C}\setminus\mathcal{S}),$ (2.48)
$\displaystyle\leq\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+2\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}\setminus\mathcal{S})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}\cap\mathcal{S})+1,$ (2.49)

where the first inequality follows from Lemma 2.11, the second inequality
follows from Lemma 2.13 and $\Delta_{\min}\leq\overline{\gamma}_{\rm
inc}^{-1}\Delta_{0}$, and the last line follows from Lemma 2.15 and
$\mathcal{VS}\subset\mathcal{S}$. Now we use Lemma 2.12 with
$\Delta_{\min}\leq\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})\overline{\gamma}_{\rm inc}^{-1}\Delta^{*}(\epsilon)$ to get

$\displaystyle\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min}))$
$\displaystyle\leq\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}\cap\mathcal{S})+1,$ (2.50)
$\displaystyle=\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{1}\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\cap\mathcal{S})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{D}^{C}(\max(\gamma_{C},\gamma_{F},\gamma_{\rm
dec})^{-1}\overline{\gamma}_{\rm
inc}\Delta_{\min})\cap\mathcal{L}\cap\mathcal{S})+1,$ (2.51)
$\displaystyle\leq\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{1}\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))+\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+1,$
(2.52)
$\displaystyle\leq(C_{1}+2)\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))+1,$
(2.53)

where the last line follows from Lemma 2.11.

Separately, we use Lemma 2.15, and apply Lemma 2.12 with
$\Delta_{\min}\leq\Delta^{*}(\epsilon)$ to get

$\displaystyle\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))$
$\displaystyle=\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\setminus\mathcal{VS})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}^{C}\setminus\mathcal{VS}),$
(2.54)
$\displaystyle\leq\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\setminus\mathcal{VS})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L})+1,$
(2.55)
$\displaystyle=\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+2\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\setminus\mathcal{VS})$
$\displaystyle\qquad+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\cap\mathcal{VS})+1,$
(2.56)
$\displaystyle=\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{L}\cap\mathcal{VS})+1,$
(2.57) $\displaystyle\leq
2\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+1.$ (2.58)

We then get

$\displaystyle\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))$
$\displaystyle\leq 2\\#(\mathcal{A}\cap\mathcal{D}^{C}(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min})\cap\mathcal{VS})$
$\displaystyle\qquad+2\\#(\mathcal{A}\cap\mathcal{D}(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min})\cap\mathcal{D}^{C}(\Delta_{\min})\cap\mathcal{VS})+1,$
(2.59) $\displaystyle\leq 2\\#(\mathcal{D}^{C}(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min})\cap\mathcal{VS})$
$\displaystyle\qquad+2\\#(\mathcal{A}\cap\mathcal{D}(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min})\cap\mathcal{VS})+1,$ (2.60) $\displaystyle\leq
2C_{3}\\#(\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})+2\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)+1,$ (2.61)
$\displaystyle=2C_{3}\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})+2C_{3}\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})$
$\displaystyle\qquad+2\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)+1,$ (2.62) $\displaystyle\leq
2C_{3}\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))+2C_{3}\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})$
$\displaystyle\qquad+2\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)+1,$ (2.63)

where the third inequality follows from Lemma 2.11 and Lemma 2.14 with

$\Delta_{\min}\leq\overline{\gamma}_{\rm inc}^{-1}\Delta_{0}\leq\gamma_{\rm
inc}^{-1}\Delta_{0}\leq\min(\Delta_{0},\gamma_{\rm
inc}^{-1}\Delta_{\max})\leq\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-1}\min(\Delta_{0},\gamma_{\rm inc}^{-1}\Delta_{\max}).$
(2.64)

Since $C_{3}\in(0,1/2)$, we can rearrange (2.63) to conclude that

$\displaystyle\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))$
$\displaystyle\qquad\qquad\leq\frac{1}{1-2C_{3}}\left[2C_{3}\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})+2\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)+1\right].$ (2.65)

Now, we combine (2.53) and (2.65) to get

$\displaystyle\\#(\mathcal{A})$
$\displaystyle=\\#(\mathcal{A}\cap\mathcal{D}(\Delta_{\min}))+\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min})),$
(2.66)
$\displaystyle\leq(C_{1}+2)\phi(\Delta_{\min},\epsilon)+C_{1}\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+C_{2}+2\\#(\mathcal{A}\cap\mathcal{D}^{C}(\Delta_{\min}))+1,$
(2.67)
$\displaystyle\leq(C_{1}+2)\phi(\Delta_{\min},\epsilon)+\frac{4\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)}{1-2C_{3}}+C_{1}\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})$
$\displaystyle\qquad+C_{2}+\frac{2}{1-2C_{3}}+1+\frac{4C_{3}}{1-2C_{3}}\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS}),$
(2.68)
$\displaystyle\leq(C_{1}+2)\phi(\Delta_{\min},\epsilon)+\frac{4\phi(\gamma_{\rm
inc}^{-1}\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})\Delta_{\min},\epsilon)}{1-2C_{3}}+C_{2}+\frac{2}{1-2C_{3}}+1$
$\displaystyle\qquad+\max\left(C_{1},\frac{4C_{3}}{1-2C_{3}}\right)\left[\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})\right].$
(2.69)

Since $\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S}$ and
$\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS}$ are
disjoint subsets of $\mathcal{A}^{C}$, we have

$\\#(\mathcal{A}^{C}\cap\mathcal{D}(\Delta_{\min})\cap\mathcal{S})+\\#(\mathcal{A}^{C}\cap\mathcal{D}^{C}(\Delta_{\min})\setminus\mathcal{VS})\leq\\#(\mathcal{A}^{C})=(K+1)-\\#(\mathcal{A}).$
(2.70)

Substituting this into (2.69) and rearranging, we get the desired result. That
$C_{4}>0$ follows from $C_{1}>0$ and $C_{3}\in(0,1/2)$. ∎

###  2.4 Overall Complexity Bound

The key remaining step is to compare $\\#(\mathcal{A})$ with $K$. Since each
event “$Q_{k}$ is well aligned” is effectively an independent Bernoulli trial
with success probability at least $1-\delta_{S}$, we derive the below result
based on a concentration bound for Bernoulli trials [23, Lemma 2.1].

######  Lemma 2.17.

Suppose Assumptions 2.2, 2.3, 2.4 and 2.9 hold. Then we have

$\mathbb{P}\left[\\#(\mathcal{A})+1\leq(1-\delta_{S})(1-\delta)(K+1)\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4},$ (2.71)

for all $\delta\in(0,1)$.

###### Proof.

The CHECK_MODEL=FALSE case of this proof has a general framework based on [42,
Lemma 4.5]—also followed in [16]—with a probabilistic argument from [23, Lemma
2.1].

First, we consider only the subsequence of iterations
$\mathcal{K}_{1}:=\\{k_{0},\ldots,k_{J}\\}\subset\\{0,\ldots,K\\}$ when
$Q_{k}$ is resampled (i.e. where CHECK_MODEL=FALSE, so $Q_{k}\neq Q_{k-1}$).
For convenience, we define $\mathcal{A}_{1}:=\mathcal{A}\cap\mathcal{K}_{1}$
and $\mathcal{A}^{C}_{1}:=\mathcal{A}^{C}\cap\mathcal{K}_{1}$.

Let $T_{k_{j}}$ be the indicator function for the event “$Q_{k_{j}}$ is well-
aligned”, and so $\\#(\mathcal{A}_{1})=\sum_{j=0}^{J}T_{k_{j}}$. Since
$T_{k_{j}}\in\\{0,1\\}$, and denoting
$p_{k_{j}}:=\mathbb{P}\left[T_{k_{j}}=1|\,\bm{x}_{k_{j}}\right]$, for any
$t>0$ we have

$\mathbb{E}\left[e^{-t(T_{k_{j}}-p_{k_{j}})}|\,\bm{x}_{k_{j}}\right]=p_{k_{j}}e^{-t(1-p_{k_{j}})}+(1-p_{k_{j}})e^{tp_{k_{j}}}=e^{tp_{k_{j}}+\log(1-p_{k_{j}}+p_{k_{j}}e^{-t})}\leq
e^{t^{2}p_{k_{j}}/2},$ (2.72)

where the inequality from the identity $px+\log(1-p+pe^{-x})\leq px^{2}/2$,
for all $p\in[0,1]$ and $x\geq 0$, shown in [23, Lemma 2.1].

Using the tower property of conditional expectations and the fact that, since
$k_{j}\in\mathcal{K}_{1}$, $T_{k_{j}}$ only depends on $\bm{x}_{k_{j}}$ and
not any previous iteration, we then get

$\displaystyle\mathbb{E}\left[e^{-t(\\#(\mathcal{A}_{1})-\sum_{j=0}^{J}p_{k_{j}})}\right]$
$\displaystyle=\mathbb{E}\left[e^{-t\sum_{j=0}^{J}(T_{k_{j}}-p_{k_{j}})}\right],$
(2.73)
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[e^{-t\sum_{j=0}^{J}(T_{k_{j}}-p_{k_{j}})}|\,Q_{0},\ldots,Q_{k_{J}-1},\bm{x}_{0},\ldots,\bm{x}_{k_{J}}\right]\right],$
(2.74)
$\displaystyle=\mathbb{E}\left[e^{-t\sum_{j=0}^{J-1}(T_{k_{j}}-p_{k_{j}})}\mathbb{E}\left[e^{-t(T_{k_{J}}-p_{k_{J}})}|\,Q_{0},\ldots,Q_{k_{J}-1},\bm{x}_{0},\ldots,\bm{x}_{k_{J}}\right]\right],$
(2.75)
$\displaystyle=\mathbb{E}\left[e^{-t\sum_{j=0}^{J-1}(T_{k_{j}}-p_{k_{j}})}\mathbb{E}\left[e^{-t(T_{k_{J}}-p_{k_{J}})}|\,\bm{x}_{k_{J}}\right]\right],$
(2.76) $\displaystyle\leq
e^{t^{2}p_{k_{J}}/2}\mathbb{E}\left[e^{-t\sum_{j=0}^{J-1}(T_{k_{j}}-p_{k_{j}})}\right],$
(2.77) $\displaystyle\leq e^{t^{2}(\sum_{j=0}^{J}p_{k_{j}})/2},$ (2.78)

where the second-last line follows from (2.72) and the last line follows by
induction. This means that

$\displaystyle\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq\sum_{j=0}^{J}p_{k_{j}}-\lambda\right]$
$\displaystyle=\mathbb{P}\left[e^{-t\left(\\#(\mathcal{A}_{1})-\sum_{j=0}^{J}p_{k_{j}}\right)}>e^{t\lambda}\right],$
(2.79) $\displaystyle\leq
e^{-t\lambda}\mathbb{E}\left[e^{-t\left(\\#(\mathcal{A}_{1})-\sum_{j=0}^{J}p_{k_{j}}\right)}\right],$
(2.80) $\displaystyle\leq
e^{t^{2}\left(\sum_{j=0}^{J}p_{k_{j}}\right)/2-t\lambda},$ (2.81)

where the inequalities follow from Markov’s inequality and (2.78)
respectively. Taking $t=\lambda/\sum_{j=0}^{J}p_{k_{j}}$, we get

$\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq\sum_{j=0}^{J}p_{k_{j}}-\lambda\right]\leq
e^{-\lambda^{2}/\left(2\sum_{j=0}^{J}p_{k_{j}}\right)}.$ (2.82)

Finally, we take $\lambda=\delta\sum_{j=0}^{J}p_{k_{j}}$ for some
$\delta\in(0,1)$ and note that $p_{k_{j}}\geq(1-\delta_{S})$ (from Assumption
2.9), to conclude

$\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq(1-\delta)(1-\delta_{S})(J+1)\right]\leq\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq(1-\delta)\sum_{j=0}^{J}p_{k_{j}}\right]\leq
e^{-\delta^{2}\left(\sum_{j=0}^{J}p_{k_{j}}\right)/2},$ (2.83)

or equivalently, using the partition
$\mathcal{K}_{1}=\mathcal{A}_{1}\cup\mathcal{A}_{1}^{C}$,

$\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq(1-\delta)(1-\delta_{S})[\\#(\mathcal{A}_{1})+\\#(\mathcal{A}_{1}^{C})]\right]\leq
e^{-\delta^{2}(1-\delta_{S})[\\#(\mathcal{A}_{1})+\\#(\mathcal{A}_{1}^{C})]/2}.$
(2.84)

Now we must consider the iterations for which CHECK_MODEL=TRUE (so
$Q_{k}=Q_{k-1}$), which we denote $\mathcal{K}_{1}^{C}$. The algorithm ensures
that if $k\in\mathcal{K}_{1}^{C}$, then $k+1\in\mathcal{K}_{1}$ (unless we are
in the last iteration we consider, $k=K$). Futher, the algorithm guarantees
that if $k\in\mathcal{K}_{1}^{C}$, then $k>0$ and $k\in\mathcal{A}$ if and
only if $k-1\in\mathcal{A}$. These are the key implications of RSDFO that we
will now use.

Firstly, we have $\\#(\mathcal{K}_{1}^{C})\leq\\#(\mathcal{K}_{1})+1$, and so

$K+1=\\#(\mathcal{K}_{1})+\\#(\mathcal{K}_{1}^{C})\leq
2[\\#(\mathcal{A}_{1})+\\#(\mathcal{A}_{1}^{C})]+1,$ (2.85)

which means (2.84) becomes

$\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq(1-\delta)(1-\delta_{S})[\\#(\mathcal{A}_{1})+\\#(\mathcal{A}_{1}^{C})]\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4}.$ (2.86)

Setting $\alpha:=\delta+\delta_{S}+\delta\delta_{S}$, we have
$(1-\delta)(1-\delta_{S})=1-\alpha$, and so

$\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq\frac{1-\alpha}{\alpha}\\#(\mathcal{A}_{1}^{C})\right]=\mathbb{P}\left[\\#(\mathcal{A}_{1})\leq(1-\alpha)[\\#(\mathcal{A}_{1})+\\#(\mathcal{A}_{1}^{C})]\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4}.$ (2.87)

Secondly, we have
$\\#(\mathcal{K}_{1}^{C}\cap\mathcal{A}^{C})\leq\\#(\mathcal{A}_{1}^{C})+1$,
and so $\\#(\mathcal{A}^{C})\leq 2\\#(\mathcal{A}_{1}^{C})+1$. This and
$\mathcal{A}_{1}\subset\mathcal{A}$ give

$\mathbb{P}\left[\\#(\mathcal{A})\leq\frac{1-\alpha}{2\alpha}[\\#(\mathcal{A}^{C})-1]\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4}.$ (2.88)

We then note that $K+1=\\#(\mathcal{A})+\\#(\mathcal{A}^{C})$, and so

$\displaystyle\mathbb{P}\left[\\#(\mathcal{A})\leq\frac{1-\alpha}{2\alpha}[K+1-\\#(\mathcal{A})-1]\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4},$ (2.89)
$\displaystyle\mathbb{P}\left[\\#(\mathcal{A})+\frac{1-\alpha}{1+\alpha}\leq\frac{1-\alpha}{1+\alpha}(K+1)\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4},$ (2.90)
$\displaystyle\mathbb{P}\left[\\#(\mathcal{A})+1\leq(1-\alpha)(K+1)\right]\leq
e^{-\delta^{2}(1-\delta_{S})K/4},$ (2.91)

since $\alpha>0$. ∎

######  Theorem 2.18.

Suppose Assumptions 2.2, 2.3, 2.4 and 2.9 hold, and we have $\beta_{F}\leq
c_{2}$, $\delta_{S}<1/(1+C_{4})$ for $C_{4}$ defined in Lemma 2.16, and
$\gamma_{\rm inc}>\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-2}$. Then for any $\epsilon>0$ and

$k\geq\frac{2(\psi(\epsilon)+1)}{1-\delta_{S}-C_{4}/(1+C_{4})},$ (2.92)

we have

$\mathbb{P}\left[\min_{j\leq k}\|\nabla f(\bm{x}_{j})\|\leq\epsilon\right]\geq
1-\exp\left(-k\frac{(1-\delta_{S}-C_{4}/(1+C_{4}))^{2}}{16(1-\delta_{S})}\right).$
(2.93)

Alternatively, if $K_{\epsilon}:=\min\\{k:\|\nabla
f(\bm{x}_{k})\|\leq\epsilon\\}$ for any $\epsilon>0$, then

$\mathbb{P}\left[K_{\epsilon}\leq\left\lceil\frac{2(\psi(\epsilon)+1)}{1-\delta_{S}-C_{4}/(1+C_{4})}\right\rceil\right]\geq
1-\exp\left(-\frac{(\psi(\epsilon)+1)[1-\delta_{S}-C_{4}/(1+C_{4})]}{8(1-\delta_{S})}\right),$
(2.94)

where $\psi(\epsilon)$ is defined in Lemma 2.16.

###### Proof.

First, fix some arbitrary $k\geq 0$. Let $\epsilon_{k}:=\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|$ and $A_{k}$ be the number of well-aligned iterations in
$\\{0,\ldots,k\\}$. If $\epsilon_{k}>0$, from Lemma 2.16, we have

$A_{k}\leq\psi(\epsilon_{k})+\frac{C_{4}}{1+C_{4}}(k+1).$ (2.95)

For any $\delta>0$ such that

$\delta<1-\frac{C_{4}}{(1+C_{4})(1-\delta_{S})},$ (2.96)

we have $(1-\delta_{S})(1-\delta)>C_{4}/(1+C_{4})$, and so we can compute

$\displaystyle\mathbb{P}\left[\psi(\epsilon_{k})\leq\left[(1-\delta_{S})(1-\delta)-\frac{C_{4}}{1+C_{4}}\right](k+1)-1\right]$
$\displaystyle\leq\mathbb{P}\left[A_{k}\leq(1-\delta_{S})(1-\delta)(k+1)\right],$
(2.97) $\displaystyle\leq e^{-\delta^{2}(1-\delta_{S})k/4},$ (2.98)

using Lemma 2.17. Defining

$\delta:=\frac{1}{2}\left[1-\frac{C_{4}}{(1+C_{4})(1-\delta_{S})}\right],$
(2.99)

we have

$(1-\delta_{S})(1-\delta)=\frac{1}{2}\left[1-\delta_{S}+\frac{C_{4}}{1+C_{4}}\right]>\frac{C_{4}}{1+C_{4}},$
(2.100)

since $1-\delta_{S}>C_{4}/(1+C_{4})$ from our assumption on $\delta_{S}$.
Hence we get

$\mathbb{P}\left[\psi(\epsilon_{k})\leq\frac{1}{2}\left(1-\delta_{S}-\frac{C_{4}}{1+C_{4}}\right)(k+1)-1\right]\leq
e^{-k[1-\delta_{S}-C_{4}/(1+C_{4})]^{2}/\left[16(1-\delta_{S})\right]},$
(2.101)

and we note that this result is still holds if $\epsilon_{k}=0$, as
$\lim_{\epsilon\to 0}\psi(\epsilon)=\infty$.

Now we fix $\epsilon>0$ and choose $k$ satisfying (2.92). We use the fact that
$\psi(\cdot)$ is non-increasing to get

$\displaystyle\mathbb{P}\left[\epsilon_{k}\geq\epsilon\right]$
$\displaystyle\leq\mathbb{P}\left[\psi(\epsilon_{k})\leq\psi(\epsilon)\right],$
(2.102)
$\displaystyle\leq\mathbb{P}\left[\psi(\epsilon_{k})\leq\frac{1}{2}(1-\delta_{S}-C_{4}/(1+C_{4}))k-1\right],$
(2.103)
$\displaystyle\leq\mathbb{P}\left[\psi(\epsilon_{k})\leq\frac{1}{2}(1-\delta_{S}-C_{4}/(1+C_{4}))(k+1)-1\right],$
(2.104)

and (2.93) follows. Lastly, we fix

$k=\left\lceil\frac{2(\psi(\epsilon)+1)}{1-\delta_{S}-C_{4}/(1+C_{4})}\right\rceil,$
(2.105)

and we use (2.93) and the definition of $K_{\epsilon}$ to get

$\displaystyle\mathbb{P}\left[K_{\epsilon}\geq k\right]$
$\displaystyle=\mathbb{P}\left[\epsilon_{k}\geq\epsilon\right],$ (2.106)
$\displaystyle\leq
e^{-k[1-\delta_{S}-C_{4}/(1+C_{4})]^{2}/[16(1-\delta_{S})]},$ (2.107)
$\displaystyle\leq\exp\left(-\frac{(\psi(\epsilon)+1)[1-\delta_{S}-C_{4}/(1+C_{4})]}{8(1-\delta_{S})}\right),$
(2.108)

and we get (2.94). ∎

######  Corollary 2.19.

Suppose the assumptions of Theorem 2.18 hold. Then for $k\geq k_{0}$ for some
$k_{0}$, we have

$\mathbb{P}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|\leq\frac{C}{\sqrt{k}}\right]\geq 1-e^{-ck},$ (2.109)

for some constants $c,C>0$. Alternatively, for $\epsilon\in(0,\epsilon_{0})$
for some $\epsilon_{0}$, we have

$\mathbb{P}\left[K_{\epsilon}\leq\widetilde{C}\epsilon^{-2}\right]\geq
1-e^{-\widetilde{c}\epsilon^{-2}},$ (2.110)

for constants $\widetilde{c},\widetilde{C}>0$.

###### Proof.

For $\epsilon$ sufficiently small, both $\epsilon_{g}(\epsilon)$ and
$\Delta_{\min}(\epsilon)$ are equal to a multiple of $\epsilon$, and so
$\psi(\epsilon)=\alpha_{1}\epsilon^{-2}+\alpha_{2}=\Theta(\epsilon^{-2})$, for
some constants $\alpha_{1},\alpha_{2}>0$.

Therefore for $k$ sufficiently large, the choice

$\epsilon=\sqrt{\frac{2\alpha_{1}}{(1-\delta_{S}-C_{4}/(1+C_{4}))k-2-2\alpha_{2}}}=\Theta(k^{-1/2}),$
(2.111)

is sufficiently small that
$\psi(\epsilon)=\alpha_{1}\epsilon^{-2}+\alpha_{2}$, and gives (2.92) with
equality. The first result then follows from (2.93).

The second result follows immediately from
$\psi(\epsilon)=\Theta(\epsilon^{-2})$ and (2.94). ∎

######  Remark 2.20.

All the above analysis holds with minimal modifications if we replace the
trust-region mechanisms in RSDFO with more standard trust-region updating
mechanisms. This includes, for example, having no safety step (i.e.
$\beta_{F}=0$), and replacing (2.8) with

$\bm{x}_{k+1}=\begin{cases}\bm{x}_{k}+\bm{s}_{k},&\rho_{k}\geq\eta,\\\
\bm{x}_{k},&\rho_{k}<\eta,\end{cases}\quad\text{and}\quad\Delta_{k+1}=\begin{cases}\min(\gamma_{\rm
inc}\Delta_{k},\Delta_{\rm max}),&\rho_{k}\geq\eta,\\\ \gamma_{\rm
dec}\Delta_{k},&\rho_{k}<\eta,\end{cases}$ (2.112)

for some $\eta\in(0,1)$. The corresponding requirement on the trust-region
updating parameters to prove a version of Theorem 2.18 is simply $\gamma_{\rm
inc}>\gamma_{\rm dec}^{-2}$ (provided we also set $\gamma_{C}=\gamma_{\rm
dec}$).

###  2.5 Remarks on Complexity Bound

Our final complexity bounds for RSDFO in Corollary 2.19 are comparable to
probabilistic direct search [42, Corollary 4.9]. They also match—in the order
of $\epsilon$—the standard bounds for (full space) model-based DFO methods for
general objective [71, 33] and nonlinear least-squares [19] problems.

Following [42], we may also derive complexity bounds on the expected first-
order optimality measure (of $\mathcal{O}(k^{-1/2})$) and the expected worst-
case complexity (of $\mathcal{O}(\epsilon^{-2})$ iterations) for RSDFO.

######  Theorem 2.21.

Suppose the assumptions of Theorem 2.18 hold. Then for $k\geq k_{0}$, the
iterates of RSDFO satisfy

$\mathbb{E}\left[\min_{j\leq k}\|\nabla f(\bm{x}_{j})\|\right]\leq
Ck^{-1/2}+\|\nabla f(\bm{x}_{0})\|e^{-ck},$ (2.113)

for $c,C>0$ from (2.109), and for $\epsilon\in(0,\epsilon_{0})$ we have

$\mathbb{E}\left[K_{\epsilon}\right]\leq\widetilde{C}_{1}\epsilon^{-2}+\frac{1}{\widetilde{c}_{1}},$
(2.114)

for constants $\widetilde{c}_{1},\widetilde{C}_{1}>0$. Here, $k_{0}$ and
$\epsilon_{0}$ are the same as in Corollary 2.19.

###### Proof.

First, for $k\geq k_{0}$ define the random variable $H_{k}$ as

$H_{k}:=\begin{cases}Ck^{-1/2},&\text{if $\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|\leq Ck^{-1/2}$},\\\ \|\nabla
f(\bm{x}_{0})\|&\text{otherwise}.\end{cases}$ (2.115)

Then since $\min_{j\leq k}\|\nabla f(\bm{x}_{j})\|\leq H_{k}$, we get

$\mathbb{E}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|\right]\leq\mathbb{E}\left[H_{k}\right]\leq Ck^{-1/2}+\|\nabla
f(\bm{x}_{0})\|\>\mathbb{P}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|>Ck^{-1/2}\right],$ (2.116)

and we get the first result by applying Corollary 2.19.

Next, if $\epsilon\in(0,\epsilon_{0})$ then

$k\geq
k_{0}(\epsilon):=\frac{2(\psi(\epsilon)+1)}{1-\delta_{S}-C_{4}/(1+C_{4})}=\Theta(\epsilon^{-2}),$
(2.117)

and so from Theorem 2.18 we have

$\mathbb{P}\left[K_{\epsilon}\leq k\right]=\mathbb{P}\left[\min_{j\leq
k}\|\nabla f(\bm{x}_{j})\|\leq\epsilon\right]\geq 1-e^{-\widetilde{c}_{1}k},$
(2.118)

where
$\widetilde{c}_{1}:=(1-\delta_{S}-C_{4}/(1+C_{4}))^{2}/[16(1-\delta_{S})]$. We
use the identity
$\mathbb{E}\left[X\right]=\int_{0}^{\infty}\mathbb{P}\left[X>t\right]dt$ for
non-negative random variables $X$ (e.g. [69, eqn. (1.9)]) to get

$\mathbb{E}\left[K_{\epsilon}\right]\leq
k_{0}(\epsilon)+\int_{k_{0}(\epsilon)}^{\infty}\mathbb{P}\left[K_{\epsilon}>t\right]dt\leq
k_{0}(\epsilon)+\sum_{k=k_{0}(\epsilon)}^{\infty}e^{-\widetilde{c}_{1}k}=k_{0}(\epsilon)+\frac{e^{-\widetilde{c}_{1}k_{0}(\epsilon)}}{1-e^{-\widetilde{c}_{1}}},$
(2.119)

where $\widetilde{C}_{1}$ comes from $k_{0}(\epsilon)=\Theta(\epsilon^{-2})$,
which concludes our proof. ∎

Furthermore, we also get almost-sure convergence of $\liminf$ type, similar to
[26, Theorem 10.12] in the deterministic case.

######  Theorem 2.22.

Suppose the assumptions of Theorem 2.18 hold. Then the iterates of RSDFO
satisfy $\inf_{k\geq 0}\|\nabla f(\bm{x}_{k})\|=0$ almost surely.

###### Proof.

From Theorem 2.18, for any $\epsilon>0$ we have

$\lim_{k\to\infty}\mathbb{P}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|>\epsilon\right]=0.$ (2.120)

However, $\mathbb{P}\left[\inf_{k\geq 0}\|\nabla
f(\bm{x}_{k})\|>\epsilon\right]\leq\mathbb{P}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|>\epsilon\right]$ for all $k$, and so

$\mathbb{P}\left[\inf_{k\geq 0}\|\nabla f(\bm{x}_{k})\|>\epsilon\right]=0.$
(2.121)

The result follows from the union bound applied to any sequence
$\epsilon_{n}=n^{-1}$, for example. ∎

In particular, if $\|\nabla f(\bm{x}_{k})\|>0$ for all $k$, then Theorem 2.22
implies $\liminf_{k\to\infty}\|\nabla f(\bm{x}_{k})\|=0$ almost surely.

###  2.6 Selecting a Subspace Dimension

We now specify how to generate our subspaces $Q_{k}$ to be probabilistically
well-aligned and uniformly bounded (Assumption 2.9). These requirements are
quite weak, and so there are several possible approaches for constructing
$Q_{k}$. First, we discuss the case where $Q_{k}$ is chosen to be a random
matrix with orthonormal columns, and show that we need to choose our subspace
dimension $p\sim\sqrt{n}$. This is related to how we ultimately select $Q_{k}$
in the practical implementation DFBGN (Section 4). However, we then show that
by instead taking $Q_{k}$ to be a Johnson-Lindenstrauss transform, we can
choose a value of $p$ independent of the ambient dimension $n$.

####  2.6.1 Random Orthogonal Basis

First, we consider the case where $Q_{k}$ is a randomly generated matrix with
orthonormal columns. We have the below result, a consequence of [54, Theorem
9].

######  Theorem 2.23.

Suppose the columns of $Q_{k}\in\mathbb{R}^{n\times p}$ form an orthonormal
basis for a randomly generated $p$-dimensional subspace of $\mathbb{R}^{n}$.
Then for any fixed vector $\bm{v}$,

$\mathbb{P}\left[\|Q_{k}^{T}\bm{v}\|>\left(\frac{p}{n}-\theta\sqrt{\frac{p}{n}}\right)\|\bm{v}\|\right]>1-3\exp(-p\theta^{2}/64),$
(2.122)

for all $\theta\in(0,1)$.

For some $\delta_{S}\in(0,1)$, if we set

$\theta=\frac{8\sqrt{\log(3/\delta_{S})}}{\sqrt{p}},$ (2.123)

then we get

$\mathbb{P}\left[\|Q_{k}^{T}\nabla
f(\bm{x}_{k})\|>\left(\frac{p}{n}-\theta\sqrt{\frac{p}{n}}\right)\|\nabla
f(\bm{x}_{k})\|\right]>1-\delta_{S}.$ (2.124)

Therefore Assumption 2.9 is achieved provided we set

$p\geq\alpha_{Q}n+\sqrt{64\log(3/\delta_{S})n\,}.$ (2.125)

In particular, since our theory holds if we take $\alpha_{Q}>0$ arbitrarily
small, our only limitation on the subspace dimension $p$ is
$\delta_{S}<1/(1+C_{4})$ from Theorem 2.18, yielding the minimum requirement

$p>\sqrt{64\log(3+3C_{4})n\,},$ (2.126)

where $C_{4}$ depends only on our choice of trust-region algorithm parameters.
The boundedness condition Assumption 2.9(b) holds with $Q_{\max}=1$
automatically. This gives us considerable scope to use very small subspace
dimensions, which means our subspace approach has a strong chance of providing
a substantial reduction in linear algebra costs.

####  2.6.2 Johnson-Lindenstrauss Embeddings

We can improve on the requirement (2.126) on $p$ by using $Q_{k}$ with non-
orthonormal columns. Specifically, we take $Q_{k}$ to be a Johnson-
Lindenstrauss transform (JLT) [74]. The application of these techniques to
random subspace optimization algorithms follows [15, 16].

######  Definition 2.24.

A random matrix $S\in\mathbb{R}^{p\times n}$ is an $(\epsilon,\delta)$-JLT if,
for any point $\bm{v}\in\mathbb{R}^{n}$, we have

$\displaystyle\mathbb{P}\left[(1-\epsilon)\|\bm{v}\|^{2}\leq\|S\bm{v}\|^{2}\leq(1+\epsilon)\|\bm{v}\|^{2}\right]\geq
1-\delta.$ (2.127)

There have been many different approaches for constructing
$(\epsilon,\delta)$-JLT matrices proposed. Two common examples are:

  * •

If $S$ is a random Gaussian matrix with independent entries $S_{i,j}\sim
N(0,1/p)$ and $p=\Omega(\epsilon^{-2}|\log\delta|)$, then $S$ is an
$(\epsilon,\delta)$-JLT (see [12, Theorem 2.13], for example).

  * •

We say that $S$ is an $s$-hashing matrix if it has exactly $s$ nonzero entries
per column (indices sampled independently), which take values $\pm 1/\sqrt{s}$
selected independently with probability 1/2. If $S$ is an $s$-hashing matrix
with $s=\Theta(\epsilon^{-1}|\log\delta|)$ and
$p=\Omega(\epsilon^{-2}|\log\delta|)$, then $S$ is an $(\epsilon,\delta)$-JLT
[48].

By taking $\bm{v}=\nabla f(\bm{x}_{k})$ in iteration $k$, and noting
$(1-\epsilon)^{2}\leq 1-\epsilon$ for all $\epsilon\in(0,1)$, we have that
Assumption 2.9(a) holds if we take $Q_{k}=S^{T}$, where $S$ is any
$(1-\alpha_{Q},\delta_{S})$-JLT. That is, Assumption 2.9(a) is satisfied using
either of the constructions above and
$p=\Omega((1-\alpha_{Q})^{-2}|\log\delta_{S}|)$. Importantly, these
constructions allow us to choose a subspace dimension $p$ which has _no
dependence on the ambient dimension_ $n$ (i.e. $p=\mathcal{O}(1)$ as
$n\to\infty$), an improvement on (2.126). We note that the requirement
$\delta_{S}<1/(1+C_{4})$ in Theorem 2.18 yields a very mild dependence of $p$
on the choice of trust-region updating parameters.

We conclude by noting that the uniform boundedness property Assumption 2.9(b)
is trivial if $S$ is a hashing matrix. If $S$ is Gaussian, by Bernstein’s
inequality we can choose $Q_{\max}$ large enough that
$\mathbb{P}\left[\|S\|>Q_{\max}\right]$ is small. Then, we generate our final
$Q_{k}$ by sampling $S$ independently until $\|S\|\leq Q_{\max}$ holds (which
almost-surely takes finite time). A union bound argument then gives us
Assumption 2.9(a) for this choice of $Q_{k}$, and so Assumption 2.9 is
completely satisfied with the same (asymptotic) requirements on $p$.

##  3 Random Subspace Nonlinear Least-Squares Method

We now describe how RSDFO (Algorithm 1) can be specialized to the
unconstrained nonlinear least-squares problem

$\min_{\bm{x}\in\mathbb{R}^{n}}f(\bm{x}):=\frac{1}{2}\|\bm{r}(\bm{x})\|^{2}=\frac{1}{2}\sum_{i=1}^{m}r_{i}(\bm{x})^{2},$
(3.1)

where $\bm{r}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is given by
$\bm{r}(\bm{x}):=[r_{1}(\bm{x}),\ldots,r_{m}(\bm{x})]^{T}$. We assume that
$\bm{r}$ is differentiable, but that access to the Jacobian
$J:\mathbb{R}^{n}\to\mathbb{R}^{m\times n}$ is not possible. In addition, we
typically assume that $m\geq n$ (regression), but everything here also applies
to the case $m<n$ (inverse problems). We now introduce the algorithm RSDFO-GN
(Randomized Subspace DFO with Gauss-Newton), which is a randomized subspace
version of a model-based DFO variant of the Gauss-Newton method [18].

Following the construction from [18], we assume that we have selected the
$p$-dimensional search space $\mathcal{Y}_{k}$ defined by
$Q_{k}\in\mathbb{R}^{n\times p}$ (as in RSDFO above). Then, we suppose that we
have evaluated $\bm{r}$ at $p+1$ points
$Y_{k}:=\\{\bm{x}_{k},\bm{y}_{1},\ldots,\bm{y}_{n}\\}\subset\mathcal{Y}_{k}$
(which typically are all close to $\bm{x}_{k}$). Since
$\bm{y}_{t}\in\mathcal{Y}_{k}$ for each $t=1,\ldots,p$, from (2.3) we have
$\bm{y}_{t}=\bm{x}_{k}+Q_{k}\hat{\bm{s}}_{t}$ for some
$\hat{\bm{s}}_{t}\in\mathbb{R}^{p}$.

Given this interpolation set, we first wish to construct a local subspace
linear model for $\bm{r}$:

$\displaystyle\bm{r}(\bm{x}_{k}+Q_{k}\hat{\bm{s}})\approx\hat{\bm{m}}_{k}(\hat{\bm{s}})=\bm{r}(\bm{x}_{k})+\hat{J}_{k}\hat{\bm{s}}.$
(3.2)

To do this, we choose the approximate subspace Jacobian
$\hat{J}_{k}\in\mathbb{R}^{m\times p}$ by requiring that $\hat{\bm{m}}_{k}$
interpolate $\bm{r}$ at our interpolation points $Y_{k}$. That is, we impose

$\displaystyle\hat{\bm{m}}_{k}(\hat{\bm{s}}_{t})=\bm{r}(\bm{y}_{t}),\qquad\forall
t=1,\ldots,p,$ (3.3)

which yields the $p\times p$ linear system (with $m$ right-hand sides)

$\displaystyle\hat{W}_{k}\hat{J}_{k}^{T}:=\begin{bmatrix}\hat{\bm{s}}_{1}^{T}\\\
\vdots\\\
\hat{\bm{s}}_{p}^{T}\end{bmatrix}\hat{J}_{k}^{T}=\begin{bmatrix}(\bm{r}(\bm{y}_{1})-\bm{r}(\bm{x}_{k}))^{T}\\\
\vdots\\\ (\bm{r}(\bm{y}_{p})-\bm{r}(\bm{x}_{k}))^{T}\end{bmatrix}.$ (3.4)

Our linear subspace model $\hat{\bm{m}}_{k}$ (3.2) naturally yields a local
subspace quadratic model for $f$, as in the classical Gauss-Newton method,
namely (c.f. (2.4)),

$f(\bm{x}_{k}+Q_{k}\hat{\bm{s}})\approx\hat{m}_{k}(\hat{\bm{s}}):=\frac{1}{2}\|\hat{\bm{m}}_{k}(\hat{\bm{s}})\|^{2}=f(\bm{x}_{k})+\hat{\bm{g}}_{k}^{T}\hat{\bm{s}}+\frac{1}{2}\hat{\bm{s}}^{T}\hat{H}_{k}\hat{\bm{s}},$
(3.5)

where $\hat{\bm{g}}_{k}:=\hat{J}_{k}^{T}\bm{r}(\bm{x}_{k})$ and
$\hat{H}_{k}:=\hat{J}_{k}^{T}\hat{J}_{k}$.

###  3.1 Constructing $Q_{k}$-Fully Linear Models

We now describe how we can achieve $Q_{k}$-fully linear models of the form
(3.5) in RSDFO-GN.

As in [18], we will need to define the Lagrange polynomials and
$\Lambda$-poisedness of an interpolation set. Given our interpolation set
$Y_{k}$ lies inside $\mathcal{Y}_{k}$, we consider the (low-dimensional)
Lagrange polynomials associated with $Y_{k}$. These are the linear functions
$\hat{\ell}_{0},\ldots,\hat{\ell}_{p}:\mathbb{R}^{p}\to\mathbb{R}$, defined by
the interpolation conditions

$\displaystyle\hat{\ell}_{t}(\hat{\bm{s}}_{t^{\prime}})=\delta_{t,t^{\prime}},\qquad\forall
t,t^{\prime}=0,\ldots,p,$ (3.6)

with the convention $\hat{\bm{s}}_{0}=\bm{0}$ corresponding to the
interpolation point $\bm{x}_{k}$. The Lagrange polynomials exist and are
unique whenever $\hat{W}_{k}$ (3.4) is invertible, which we typically ensure
through judicious updating of $Y_{k}$ at each iteration.

######  Definition 3.1.

For any $\Lambda>0$, the set $Y_{k}$ is $\Lambda$-poised in the
$p$-dimensional ball $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ if

$\max_{t=0,\ldots,p}\>\max_{\|\hat{\bm{s}}\|\leq\Delta_{k}}|\hat{\ell}_{t}(\hat{\bm{s}})|\leq\Lambda.$
(3.7)

Note that since $\hat{\ell}_{0}(\bm{0})=1$, for the set $Y_{k}$ to be
$\Lambda$-poised we require $\Lambda\geq 1$. In general, a larger $\Lambda$
indicates that $Y_{k}$ has “worse” geometry, which leads to a less accurate
approximation for $f$. This notion of $\Lambda$-poisedness (in a subspace) is
sufficient to construct $Q_{k}$-fully linear models (3.5) for $f$.

######  Lemma 3.2.

Suppose Assumption 2.9(b) holds, $J(\bm{x})$ is Lipschitz continuous, and
$\bm{r}$ and $J$ are uniformly bounded above in $\cup_{k\geq
0}B(\bm{x}_{k},\Delta_{\max})$. If $Y_{k}\subset
B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ and $Y_{k}$ is $\Lambda$-poised
in $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$, then $\hat{m}_{k}$ (3.5) is
a $Q_{k}$-fully linear model for $f$, with $\kappa_{\rm ef},\kappa_{\rm
eg}=\mathcal{O}(p^{2}\Lambda^{2})$.

###### Proof.

Consider the low-dimensional functions
$\hat{\bm{r}}:\mathbb{R}^{p}\to\mathbb{R}^{m}$ and
$\hat{f}:\mathbb{R}^{p}\to\mathbb{R}$ given by
$\hat{\bm{r}}_{k}(\hat{\bm{s}}):=\bm{r}(\bm{x}_{k}+Q_{k}\hat{\bm{s}})$ and
$\hat{f}(\hat{\bm{s}}):=\frac{1}{2}\|\hat{\bm{r}}(\hat{\bm{s}})\|^{2}$
respectively. We note that $\hat{\bm{r}}$ is continuously differentiable with
Jacobian $\hat{J}(\hat{\bm{s}})=J(\bm{x}_{k}+Q_{k}\hat{\bm{s}})Q_{k}$. Then
since $\|Q_{k}\|\leq Q_{\max}$ from Assumption 2.9(b), it is straightforward
to show that both $\hat{\bm{r}}$ and $\hat{J}$ are uniformly bounded above and
$\hat{J}$ is Lipschitz continuous (with a Lipschitz constant $Q_{\max}^{2}$
times larger than for $J(\bm{x})$).

We can then consider $\hat{\bm{m}_{k}}$ (3.2) and $\hat{m}_{k}$ (3.5) to be
interpolation models for $\hat{r}$ and $\hat{f}$ in the low-dimensional ball
$B(\bm{0},\Delta_{k})\subset\mathbb{R}^{p}$. From [18, Lemma 3.3], we conclude
that $\hat{m}_{k}$ is a fully linear model for $\hat{f}$ with constants
$\kappa_{\rm ef},\kappa_{\rm eg}=\mathcal{O}(p^{2}\Lambda^{2})$. The
$Q_{k}$-fully linear property follows immediately from this, noting that
$\nabla\hat{f}_{k}(\hat{\bm{s}})=Q_{k}^{T}\nabla
f(\bm{x}_{k}+Q_{k}\hat{\bm{s}})$. ∎

Given this result, the procedures in [26, Chapter 6] allow us to check and/or
guarantee the $\Lambda$-poisedness of an interpolation set, and we have met
all the requirements needed to fully specify RSDFO-GN.

Lastly, we note that underdetermined linear interpolation, where (3.4) is
underdetermined and solved in a minimal norm sense, has been recently shown to
yield a property similar to $Q_{k}$-full linearity [47, Theorem 3.6].

##### Complete RSDFO-GN Algorithm

A complete statement of RSDFO-GN is given in Algorithm 2. This exactly follows
RSDFO (Algorithm 1), but where we ask that the interpolation set satsifies the
conditions: $Y_{k}\subset B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ and
$Y_{k}$ is $\Lambda$-poised in $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$.
From Lemma 3.2, this is sufficient to guarantee $Q_{k}$-full linearity of
$\hat{m}_{k}$.

1:Starting point $\bm{x}_{0}\in\mathbb{R}^{n}$, initial trust region radius
$\Delta_{0}>0$, and subspace dimension $p\in\\{1,\ldots,n\\}$.

2:Parameters: maximum trust-region radius $\Delta_{\rm max}\geq\Delta_{0}$,
trust-region radius scalings $0<\gamma_{\rm dec}<1<\gamma_{\rm
inc}\leq\overline{\gamma}_{\rm inc}$, criticality constants
$\epsilon_{C},\mu>0$ and trust-region scaling $0<\gamma_{C}<1$, safety step
threshold $\beta_{F}>0$ and trust-region scaling $0<\gamma_{F}<1$, acceptance
thresholds $0<\eta_{1}\leq\eta_{2}<1$, and poisedness constant $\Lambda>1$.

3:Set flag CHECK_MODEL=FALSE.

4:for $k=0,1,2,\ldots$ do

5: if CHECK_MODEL=TRUE then

6: Set $Q_{k}=Q_{k-1}$.

7: Construct an interpolation set $Y_{k}\subset
B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ which is $Y_{k}$ is
$\Lambda$-poised in $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$.

8: Build the reduced model $\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ (3.5) by
solving (3.4).

9: else

10: Define a subspace by randomly sampling $Q_{k}\in\mathbb{R}^{n\times p}$.

11: Construct a reduced model $\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ (3.5)
by solving (3.4), where the interpolation points $Y_{k}\subset\mathcal{Y}_{k}$
need not be contained in $B(\bm{x}_{k},\Delta_{k})$ or be $\Lambda$-poised in
$B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$.

12: end if

13: Follow lines 12 to 24 of RSDFO (Algorithm 1), but replace every instance
of checking $Q_{k}$-full linearity of $\hat{m}_{k}$ in
$B(\bm{x}_{k},\Delta_{k})$ with checking that $Y_{k}\subset
B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ and $Y_{k}$ is $\Lambda$-poised
in $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$.

14:end for

Algorithm 2 RSDFO-GN (Randomized Subspace Derivative-Free Optimization with
Gauss-Newton) for solving (3.1).

###  3.2 Complexity Analysis for RSDFO-GN

We are now in a position to specialize our complexity analysis for RSDFO to
RSDFO-GN. For this, we need to impose a smoothness assumption on $\bm{r}$.

######  Assumption 3.3.

The extended level set
$\mathcal{L}:=\\{\bm{y}\in\mathbb{R}^{n}:\|\bm{y}-\bm{x}\|\leq\Delta_{\max}\>\text{for
some}\>f(\bm{x})\leq f(\bm{x}_{0})\\}$ is bounded, $\bm{r}$ is continuously
differentiable, and the Jacobian $J$ is Lipschitz continuous on $\mathcal{L}$.

This smoothness requirement allows us to immediately apply the complexity
analysis for RSDFO, yielding the following result.

######  Corollary 3.4.

Suppose Assumptions 3.3, 2.3, 2.4 and 2.9 hold, and we have $\beta_{F}\leq
c_{2}$, $\delta_{S}<1/(1+C_{4})$ for $C_{4}$ defined in Lemma 2.16, and
$\gamma_{\rm inc}>\min(\gamma_{C},\gamma_{F},\gamma_{\rm
dec},\beta_{F})^{-2}$. Then for the iterates generated by RSDFO-GN and $k$
sufficiently large,

$\mathbb{P}\left[\min_{j\leq k}\|\nabla
f(\bm{x}_{j})\|\leq\frac{C}{\sqrt{k}}\right]\geq 1-e^{-ck},$ (3.8)

for some constants $c,C>0$. Alternatively, for $\epsilon\in(0,\epsilon_{0})$
for some $\epsilon_{0}$, we have

$\mathbb{P}\left[K_{\epsilon}\leq\widetilde{C}\epsilon^{-2}\right]\geq
1-e^{-\widetilde{c}\epsilon^{-2}},$ (3.9)

for constants $\widetilde{c},\widetilde{C}>0$.

###### Proof.

Assumption 3.3 implies that both $\bm{r}$ and $J$ are uniformly bounded above
on $\mathcal{L}$, which is sufficient for Lemma 3.2 to hold. Hence, whenever
we check/ensure that $Y_{k}\subset
B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ and $Y_{k}$ is $\Lambda$-poised
in $B(\bm{x}_{k},\Delta_{k})\cap\mathcal{Y}_{k}$ we are checking/guaranteeing
that $\hat{m}_{k}$ is $Q_{k}$-fully linear in $B(\bm{x}_{k},\Delta_{k})$. In
addition, from [18, Lemma 3.2] and taking $f_{\rm low}=0$, we have that
Assumption 2.2 is satisfied. Therefore the result follows directly from
Corollary 2.19. ∎

##### Bounds on Objective Evaluations

In each iteration of RSDFO, we require at most $p+1$ objective evaluations: at
most $p$ to form the model $\hat{m}_{k}$—and this holds regardless of whether
we need $\hat{m}_{k}$ to be $Q_{k}$-fully linear or not—and one evaluation for
$\bm{x}_{k}+\bm{s}_{k}$. Hence the above bounds on $K_{\epsilon}$ also hold
for the number of objective evaluations required to first achieve $\|\nabla
f(\bm{x}_{k})\|\leq\epsilon$, up to a constant factor of $p+1$.

##### Interpolation Models for RSDFO

Using existing techniques for constructing $\Lambda$-poised interpolation sets
for quadratic interpolation (as outlined in [26]), the specific model
construction ideas presented here can also be applied to general objective
problems, and thus provide a concrete implementation of RSDFO.

###  3.3 Linear Algebra Cost of RSDFO-GN

In RSDFO-GN, the interpolation linear system (3.4) is solved in two steps,
namely: factorize the interpolation matrix $\hat{W}_{k}$, then back-solve for
each right-hand side. Thus, the cost of the linear algebra is:

  1. 1.

Model construction costs $\mathcal{O}(p^{3})$ to compute the factorization of
$\hat{W}_{k}$, and $\mathcal{O}(mp^{2})$ for the back-substitution solves with
$m$ right-hand sides; and

  2. 2.

Lagrange polynomial construction costs $\mathcal{O}(p^{3})$ in total, due to
one backsolve for each of the $p+1$ polynomials (using the pre-existing
factorization of $\hat{W}_{k}$).

By updating the factorization or $\hat{W}_{k}^{-1}$ directly (e.g. via the
Sherman-Morrison formula), we can replace the $\mathcal{O}(p^{3})$
factorization cost with a $\mathcal{O}(p^{2})$ updating cost (c.f. [62]).
However, the dominant $\mathcal{O}(mp^{2})$ model construction cost remains,
and in practice we have observed that the factorization needs to be recomputed
from scratch to avoid the accumulation of rounding errors.

In the case of a full-space method where $p=n$ such as in [18], these costs
becomes $\mathcal{O}(n^{3})$ for the factorization (or $\mathcal{O}(n^{2})$ if
Sherman-Morrison is used) plus $\mathcal{O}(mn^{2})$ for the back-solves. When
$n$ grows large, this linear algebra cost rapidly dominates the total runtime
of these algorithms and limits the efficiency of full-space methods. This
issue is discussed in more detail, with numerical results, in [65, Chapter
7.2].

In light of this discussion, we now turn our attention to building an
implementation of RSDFO-GN that has both strong performance (in terms of
objective evaluations) and low linear algebra cost.

##  4 DFBGN: An Efficient Implementation of RSDFO-GN

An important tenet of DFO is that objective evaluations are often expensive,
and so algorithms should be efficient in reusing information, hence limiting
the total objective evaluations required to achieve a given decrease. However
because we require our model to sit within our active space $\mathcal{Y}_{k}$,
we do not have a natural process by which to reuse evaluations between
iterations, when the space changes. We dedicate this section to outlining an
implementation of RSDFO-GN, which we call DFBGN (Derivative-Free Block Gauss-
Newton). DFBGN is designed to be efficient in its objective queries while
still only building low-dimensional models, and hence is also efficient in
terms of linear algebra. Specifically, we design DFBGN to achieve two aims:

  * •

Low computational cost: we want our implementation to have a per-iteration
linear algebra cost which is linear in the ambient dimension;

  * •

Efficient use of objective evaluations: our implementation should follow the
principles of other DFO methods and make progress with few objective
evaluations. In particular, we hope that, when run with ‘full-space models’
(i.e. $p=n$), our implementation should have (close to) state-of-the-art
performance.

We will assess the second point in Section 5 by comparison with DFO-LS [14] an
open-source model-based DFO Gauss-Newton solver which explores the full space
(i.e. $p=n$).

######  Remark 4.1.

As discussed in [14], DFO-LS has a mechanism to build a model with fewer than
$n+1$ interpolation points. However, in that context we modify the model so
that it varies over the whole space $\mathbb{R}^{n}$, which enables the
interpolation set to grow to the usual $n+1$ points and yield a full-
dimensional model. There, the goal is to make progress with very few
evaluations, but here our goal is scalability, so we keep our model low-
dimensional throughout and instead change the subspace at each iteration.

###  4.1 Subspace Interpolation Models

Similar to Section 3, we assume that, at iteration $k$, our interpolation set
has $p+1$ points
$\\{\bm{x}_{k},\bm{y}_{1},\ldots,\bm{y}_{p}\\}\subset\mathbb{R}^{n}$ with
$1\leq p\leq n$. However, we assume that these points are already given, and
use them to determine the space $\mathcal{Y}_{k}$ (as defined by $Q_{k}$).
That is, given

$\displaystyle W_{k}:=\begin{bmatrix}(\bm{y}_{1}-\bm{x}_{k})^{T}\\\ \vdots\\\
(\bm{y}_{p}-\bm{x}_{k})^{T}\end{bmatrix}\in\mathbb{R}^{p\times n},$ (4.1)

we compute the QR factorization

$W_{k}^{T}=Q_{k}R_{k},$ (4.2)

where $Q_{k}\in\mathbb{R}^{n\times p}$ has orthonormal columns and
$R_{k}\in\mathbb{R}^{p\times p}$ is upper triangular—and invertible provided
$W_{k}^{T}$ is full rank, which we guarantee by judicious replacement of
interpolation points. This gives us the $Q_{k}$ that defines $\mathcal{Y}_{k}$
via (2.3)—in this case $Q_{k}$ is has orthonormal columns—and in this way all
our interpolation points are in $\mathcal{Y}_{k}$.

Since each $\bm{y}_{t}\in\mathcal{Y}_{k}$, from (4.2) we have
$\bm{y}_{t}=\bm{x}_{k}+Q_{k}\hat{\bm{s}}_{t}$, where $\hat{\bm{s}}_{t}$ is the
$t$-th column of $R_{k}$. Hence we have $\hat{W}_{k}=R_{k}^{T}$ in (3.4) and
so $\hat{\bm{m}}_{k}$ (3.2) is given by solving

$\displaystyle
R_{k}^{T}\hat{J}_{k}^{T}=\begin{bmatrix}(\bm{r}(\bm{y}_{1})-\bm{r}(\bm{x}_{k}))^{T}\\\
\vdots\\\ (\bm{r}(\bm{y}_{p})-\bm{r}(\bm{x}_{k}))^{T}\end{bmatrix},$ (4.3)

via forward substitution, since $R_{k}^{T}$ is lower triangular. This
ultimately gives us our local model $\hat{m}_{k}$ via (3.5).

We reiterate that compared to RSDFO-GN, we have used the interpolation set
$Y_{k}$ to determine both $Q_{k}$ and $\hat{m}_{k}$, rather than first
sampling $Q_{k}$, then finding interpolation points
$Y_{k}\subset\mathcal{Y}_{k}$ with which to construct $\hat{m}_{k}$. This
difference is crucial in allowing the reuse of interpolation points between
iterations, and hence lowering the objective evaluation requirements of model
construction.

######  Remark 4.2.

As discussed in [65, Chapter 7.3], we can equivalently recover this
construction by asking for a full-space model
$\bm{m}_{k}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ given by
$\bm{m}_{k}(\bm{s})=\bm{r}(\bm{x}_{k})+J_{k}\bm{s}$ such that the
interpolation conditions
$\bm{m}_{k}(\bm{y}_{t}-\bm{x}_{k})=\bm{r}(\bm{y}_{t})$ are satisfied and
$J_{k}$ has minimal Frobenius norm.

###  4.2 Complete DFBGN Algorithm

A complete statement of DFBGN is given in Algorithm 3. Compared to RSDFO-GN,
we include specific steps to manage the interpolation set, which in turn
dictates the choice of subspace $\mathcal{Y}_{k}$. Specifically, one issue
with our approach is that our new iterate $\bm{x}_{k}+\bm{s}_{k}$ is in
$\mathcal{Y}_{k}$, so if we were to simply add $\bm{x}_{k}+\bm{s}_{k}$ into
the interpolation set, $\mathcal{Y}_{k}$ would not change across iterations,
and we will never explore the whole space. On the other hand, unlike RSDFO and
RSDFO-GN we do not want to completely resample $Q_{k}$ as this would require
too many objective evaluations. Instead, in DFBGN we delete a subset of points
from the interpolation set and add new directions orthogonal to the existing
directions, which ensures that $Q_{k+1}\neq Q_{k}$ in every iteration.555By
contrast, the optional growing mechanism in DFO-LS (Remark 4.1) is designed
such that $\bm{x}_{k}+\bm{s}_{k}$ is not in $\mathcal{Y}_{k}$, and so the
search space is automatically expanded at every iteration. However, this
requires an expensive SVD of $J_{k}\in\mathbb{R}^{m\times n}$ at every
iteration, and so is not suitable for our large-scale setting.

We also note that DFBGN does not include some important algorithmic features
present in RSDFO-GN, DFO-LS or other model-based DFO methods, and hence is
quite simple to state. These features are not necessary for a variety of
reasons, which we now outline.

1:Starting point $\bm{x}_{0}\in\mathbb{R}^{n}$, initial trust region radius
$\Delta_{0}>0$, subspace dimension $p\in\\{1,\ldots,n\\}$, and number of
points to drop at each iteration $p_{\rm drop}\in\\{1,\ldots,p\\}$.

2:Parameters: maximum and minimum trust-region radii $\Delta_{\rm
max}\geq\Delta_{0}>\Delta_{\rm end}>0$, trust-region radius scalings
$0<\gamma_{\rm dec}<1<\gamma_{\rm inc}\leq\overline{\gamma}_{\rm inc}$, and
acceptance thresholds $0<\eta_{1}\leq\eta_{2}<1$.

3:Select random orthonormal directions
$\bm{d}_{1},\ldots,\bm{d}_{p}\in\mathbb{R}^{n}$ using Algorithm 5, and build
initial interpolation set
$Y_{0}:=\\{\bm{x}_{0},\bm{x}_{0}+\Delta_{0}\bm{d}_{1},\ldots,\bm{x}_{0}+\Delta_{0}\bm{d}_{p}\\}$.

4:for $k=0,1,2,\ldots$ do

5: Given $\bm{x}_{k}$ and $Y_{k}$, solve (4.3) to build subspace models
$\hat{\bm{m}}_{k}:\mathbb{R}^{p}\to\mathbb{R}^{m}$ (3.2) and
$\hat{m}_{k}:\mathbb{R}^{p}\to\mathbb{R}$ (3.5).

6: Approximately solve the subspace trust-region subproblem in
$\mathbb{R}^{p}$ (2.5) and calculate the step
$\bm{s}_{k}=Q_{k}\hat{\bm{s}}_{k}\in\mathbb{R}^{n}$.

7: Evaluate $\bm{r}(\bm{x}_{k}+\bm{s}_{k})$ and calculate ratio $\rho_{k}$
(2.7).

8: Accept/reject step and update trust region radius: set $\bm{x}_{k+1}$ and
$\Delta_{k+1}$ as per (2.8).

9: if $\Delta_{k+1}\leq\Delta_{\rm end}$, terminate.

10: if $p<n$ then

11: Set $Y_{k+1}^{\rm init}=Y_{k}\cup\\{\bm{x}_{k}+\bm{s}_{k}\\}$.

12: Remove $\min(\max(p_{\rm drop},2),p)$ points from $Y_{k+1}^{\rm init}$
(without removing $\bm{x}_{k+1}$) using Algorithm 4.

13: else

14: Set $Y_{k+1}^{\rm
init}=Y_{k}\cup\\{\bm{x}_{k}+\bm{s}_{k}\\}\setminus\\{\bm{y}\\}$ for some
$\bm{y}\in Y_{k}\setminus\\{\bm{x}_{k+1}\\}$.

15: Remove $\min(\max(p_{\rm drop},1),p)$ points from $Y_{k+1}^{\rm init}$
(without removing $\bm{x}_{k+1}$) using Algorithm 4.

16: end if

17: Let $q:=p+1-|Y_{k+1}^{\rm init}|$, and generate random orthonormal vectors
$\\{\bm{d}_{1},\ldots,\bm{d}_{q}\\}$ that are also orthogonal to
$\\{\bm{y}-\bm{x}_{k+1}:\bm{y}\in Y_{k+1}^{\rm
init}\setminus\\{\bm{x}_{k+1}\\}\\}$, using Algorithm 5.

18: Set $Y_{k+1}=Y_{k+1}^{\rm
init}\cup\\{\bm{x}_{k+1}+\Delta_{k+1}\bm{d}_{1},\ldots,\bm{x}_{k+1}+\Delta_{k+1}\bm{d}_{q}\\}$.

19:end for

Algorithm 3 DFBGN: Derivative-Free Block Gauss-Newton for solving (3.1).

##### No Criticality and Safety Steps

Compared to RSDFO-GN, the implementation of DFBGN does not include criticality
(which is also the case in DFO-LS) or safety steps. These steps ultimately
function to ensure we do not have $\hat{\bm{g}}_{k}\ll\Delta_{k}$. In DFBGN,
we ensure $\Delta_{k}$ does not get too large compared to $\|\bm{s}_{k}\|$
through (2.8), while $\|\bm{s}_{k}\|$ is linked to $\|\hat{\bm{g}}_{k}\|$
through Lemma 2.5. If $\|\bm{s}_{k}\|$ is much smaller than $\Delta_{k}$ and
our step produces a poor objective decrease, then we will set
$\Delta_{k+1}\leftarrow\|\bm{s}_{k}\|$ for the next iteration. Although Lemma
2.5 allows $\|\bm{s}_{k}\|$ to be large even if $\|\bm{g}_{k}\|$ is small, in
practice we do not observe $\Delta_{k}\gg\|\bm{g}_{k}\|$ without DFBGN setting
$\Delta_{k+1}\leftarrow\|\bm{s}_{k}\|$ after a small number of iterations.

##### No Model-Improving Steps

An important feature of model-based DFO methods are model-improving
procedures, which change the interpolation set to ensure $\Lambda$-poisedness
(Definition 3.1), or equivalently ensure that the local model for $f$ is fully
linear. In RSDFO-GN for instance, model-improvement is performed when
CHECK_MODEL=TRUE, whereas in [26, Algorithm 10.1] there are dedicated model-
improvement phases.

Instead, DFBGN ensures accurate interpolation models via a geometry-aware
(i.e. $\Lambda$-poisedness aware) process for deleting interpolation points at
each iteration, where they are replaced by new points in directions (from
$\bm{x}_{k+1}$) which are orthogonal to $Q_{k}$ and selected at random. The
process for deleting interpolation points—and choosing a suitable number of
points to remove, $p_{\rm drop}$—at each iteration are considered in Sections
4.3 and 4.4 respectively. The process for generating new interpolation points,
Algorithm 5, is outlined in Section 4.3.

A downside of our approach is that the new orthogonal directions are not
chosen by minimizing a model for the objective (i.e. not attempting to reduce
the objective), as we have no information about how the objective varies
outside $\mathcal{Y}_{k}$. This is the fundamental trade-off between a
subspace approach and standard methods (such as DFO-LS); we can reduce the
linear algebra cost, but must spend objective evaluations to change the search
space between iterations.

##### Linear Algebra Cost of DFBGN

As in Section 3.3, our approach in DFBGN yields substantial reductions in the
required linear algebra costs compared to DFO-LS:

  * •

Model construction costs $\mathcal{O}(np^{2})$ for the factorization (4.2) and
$\mathcal{O}(mp^{2})$ for back-substitution solves (4.3), rather than
$\mathcal{O}(n^{3})$ and $\mathcal{O}(mn^{2})$ respectively for DFO-LS; and

  * •

Lagrange polynomial construction costs $\mathcal{O}(p^{3})$ rather than
$\mathcal{O}(n^{3})$.666Here, we use the existing factorization (4.2) and
solve with $p+1$ right-hand sides, as in Section 3.3.

As well as these reductions, we also get a smaller trust-region subproblem
(2.5)—in $\mathbb{R}^{p}$ rather than $\mathbb{R}^{n}$—and smaller memory
requirements for storing the model Jacobian: we only store $\hat{J}_{k}$ and
$Q_{k}$, requiring $\mathcal{O}((m+n)p)$ memory rather than $\mathcal{O}(mn)$
for storing the full $m\times n$ Jacobian. However, in (2.5), we do have the
extra cost of projecting $\hat{\bm{s}}_{k}\in\mathbb{R}^{p}$ into the full
space $\mathbb{R}^{n}$, which requires a multiplication by $Q_{k}$, costing
$\mathcal{O}(np)$. In addition to the reduced linear algebra costs, the
smaller interpolation set means we have a lower evaluation cost to construct
the initial model of $p+1$ evaluations (rather than $n+1$).

No particular choice of $p$ is needed for this method, and anything from $p=1$
(i.e. coordinate search) to $p=n$ (i.e. full space search) is allowed.
However, unsurprisingly, we shall see that larger values of $p$ give better
performance in terms of evaluations, except for the very low-budget phase,
where smaller values of $p$ benefit from a lower initialization cost. Hence,
we expect that our approach with small $p$ is useful when the
$\mathcal{O}(mn^{2}+n^{3})$ per-iteration linear algebra cost of DFO-LS is too
great, and reducing the linear algebra cost is worth (possibly) needing more
objective evaluations to achieve a given accuracy. As a result, $p$ should in
general be set as large as possible, given the linear algebra costs the user
is willing to bear.

In Table 1, we compare the linear algebra costs of DFO-LS and DFBGN. The
overall per-iteration cost of DFO-LS is $\mathcal{O}(mn^{2}+n^{3})$ and the
cost of DFBGN is $\mathcal{O}(mp^{2}+np^{2}+p^{3})$, depending on the choice
of $p\in\\{1,\ldots,n\\}$. The key benefit is that our dependency on the
underlying problem dimension $n$ decreases from cubic in DFO-LS to linear in
DFBGN (provided $p\ll n$). We also note that both methods have linear cost in
the number of residuals $m$, but with a factor that is significantly smaller
in DFBGN than in DFO-LS—$\mathcal{O}(p^{2})$ compared to $\mathcal{O}(n^{2})$.

Algorithm phase | DFO-LS | DFBGN | Comment
---|---|---|---
Form $\hat{\bm{m}}_{k}$ (3.2) | $\mathcal{O}(n^{3}+mn^{2})$ | $\mathcal{O}(np^{2}+mp^{2})$ | Factorization plus linear solves
Form $\hat{m}_{k}$ (3.5) | $\mathcal{O}(mn^{2})$ | $\mathcal{O}(mp^{2})$ |  Form $\hat{J}_{k}^{T}\hat{J}_{k}$
Trust-region subproblem |  $\mathcal{O}(n^{2})$–$\mathcal{O}(n^{3})$ |  $\mathcal{O}(p^{2})$–$\mathcal{O}(p^{3})$ | Depending on # CG iterations*
Calculate $\bm{s}_{k}=Q_{k}\hat{\bm{s}}_{k}$ | — | $\mathcal{O}(np)$ |
Form new step $\bm{x}_{k}+\bm{s}_{k}$ | $\mathcal{O}(n)$ | $\mathcal{O}(n)$ |
Choose point to replace | $\mathcal{O}(n^{3})$ | — | Compute Lagrange polynomials
Model improvement | $\mathcal{O}(n^{3})$ | — | Recompute Lagrange polynomials
Choose points to remove | — | $\mathcal{O}(p^{3}+np)$ |  See Algorithm 4
Generate new directions | — | $\mathcal{O}(np^{2})$ |  See Algorithm 5
Total | $\mathcal{O}(mn^{2}+n^{3})$ | $\mathcal{O}(mp^{2}+np^{2}+p^{3})$ |
Table 1: Comparison of per-iteration linear algebra costs of DFO-LS and DFBGN
(Algorithm 3) with subspace dimension $p\in\\{1,\ldots,n\\}$. *Note that the
trust-region subproblem is solved using a truncated CG method [25, Chapter
7.5.1] originally from [63] in both DFO-LS and DFBGN.

######  Remark 4.3.

In every iteration we must compute the QR factorization (4.2). However, we
note, similar to [18, Section 4.2], that adding, removing and changing
interpolation points all induce simple changes to $\hat{W}_{k}^{T}$ (adding or
removing columns, and low-rank updates). This means that (4.2) can be computed
with cost $\mathcal{O}(np)$ per iteration using the updating methods in [37,
Section 12.5]. In our implementation, however, we do not do this, as we find
that these updates introduce errors777Leading to $\hat{W}_{k}^{T}\neq
Q_{k}R_{k}$, not relating to $Q_{k}$ orthogonal or $R_{k}$ upper-triangular.
that accumulate at every iteration and reduce the accuracy of the resulting
interpolation models. To maintain the numerical performance of our method, we
need to recompute (4.2) from scratch regularly (e.g. every 10 iterations), and
so would not see the $\mathcal{O}(np)$ per-iteration cost, on average.

######  Remark 4.4.

The default parameter choices for DFBGN are the same as DFO-LS, namely:
$\Delta_{\rm max}=10^{10}$, $\Delta_{\rm end}=10^{-8}$, $\gamma_{\rm
dec}=0.5$, $\gamma_{\rm inc}=2$, $\overline{\gamma}_{\rm inc}=4$,
$\eta_{1}=0.1$, and $\eta_{2}=0.7$. DFBGN also uses the same default choice
$\Delta_{0}=0.1\max(\|\bm{x}_{0}\|_{\infty},1)$. The default choice of $p_{\rm
drop}$ is discussed in Section 4.4.

##### Adaptive Choice of $\bm{p}$

One approach that we have considered is to allow $p$ to vary between
iterations of DFBGN, rather than being constant throughout. Instead of adding
$p_{\rm drop}$ new points at the end of each iteration (line 17), we implement
a variable $p$ by adding at least one new point to the interpolation set,
continuing until some criterion is met. This criterion is designed to allow
$p$ small when such a $p$ allows us to make reasonable progress, but to grow
$p$ up to $p\approx n$ when necessary.

We have tested several possible criteria—comparing some combination of model
gradient and Hessian, trust-region radius, trust-region step length, and
predicted decrease from the trust-region step—and found the most effective to
be comparing the model gradient and Hessian with the trust-region radius.
Specifically, we continue adding new directions until (c.f. Lemma 2.10 and
[18, Lemma 3.22])

$\frac{\|\bm{g}_{k}\|}{\max(\|H_{k}\|,1)}\geq\alpha\Delta_{k},$ (4.4)

for some $\alpha>0$ (we use $\alpha=0.2(n-p)/n$ for an interpolation set with
$p+1$ points). However, our numerical testing has shown that DFBGN with $p$
fixed outperforms this approach for all budget and accuracy levels, on both
medium- and large-scale problems, and so we do not consider it further here.
We delegate further study of this approach to future work, to see if
alternative adaptive choices for $p$ can be beneficial.

###  4.3 Interpolation Set Management

We now provide more details about how we manage the interpolation set in
DFBGN. Specifically, we discuss how points are chosen for removal from
$Y_{k}$, and how new interpolations points are calculated.

####  4.3.1 Geometry Management

In the description of DFBGN, there are no explicit mechanisms to ensure that
the interpolation set is well-poised. DFBGN ensures that the interpolation set
has good geometry through two mechanisms:

  * •

We use a geometry-aware mechanism for removing points, based on [63, 18],
which requires the computation of Lagrange polynomials. This mechanism is
given in Algorithm 4, and is called in lines 12 and 15 of DFBGN, as well as to
select a point to replace in line 14; and

  * •

Adding new directions that are orthogonal to existing directions, and of
length $\Delta_{k}$, means adding these new points never causes the
interpolation set to have poor poisedness.

Together, these two mechanisms mean that any points causing poor poisedness
are quickly removed, and replaced by high-quality interpolation points
(orthogonal to existing directions, and within distance $\Delta_{k}$ of the
current iterate).

1:Interpolation set $\\{\bm{x}_{k+1},\bm{y}_{1},\ldots,\bm{y}_{p}\\}$ with
current iterate $\bm{x}_{k+1}$, trust-region radius $\Delta_{k+1}>0$ and
number of points to remove $p_{\rm drop}\in\\{1,\ldots,p\\}$.

2:Compute the (linear) Lagrange polynomials for
$\\{\bm{x}_{k+1},\bm{y}_{1},\ldots,\bm{y}_{p}\\}$ in the same way as (4.3).

3:For $t=1,\ldots,p$ (i.e. all interpolation points except $\bm{x}_{k+1}$),
compute  $\theta_{t}:=\max_{\bm{x}\in
B(\bm{x}_{k+1},\Delta_{k+1})}|\ell_{t}(\bm{x})|\cdot\max\left(\frac{\|\bm{y}_{t}-\bm{x}_{k+1}\|^{4}}{\Delta_{k+1}^{4}},1\right).$
(4.5)

4:Remove the $p_{\rm drop}$ interpolation points with the largest values of
$\theta_{t}$.

Algorithm 4 Mechanism for removing points from the interpolation set in DFBGN.

The linear algebra cost of Algorithm 4 is $\mathcal{O}(p^{3})$ to compute $p$
Lagrange polynomials with cost $\mathcal{O}(p^{2})$ each (since we already
have a factorization of $\hat{W}_{k}^{T}$). Then for each $t$ we must evaluate
$\theta_{t}$ (4.5), with cost $\mathcal{O}(p)$ to maximize $\ell_{t}(\bm{x})$
(since $\ell_{t}$ is linear and varies only in directions $\mathcal{Y}_{k}$),
and $\mathcal{O}(n)$ to calculate $\|\bm{y}_{t}-\bm{x}_{k+1}\|$. This gives a
total cost of $\mathcal{O}(p^{3}+np)$.888We could instead compute
$\|\bm{y}_{t}-\bm{x}_{k+1}\|$ by taking the norm of the $t$-th column of
$R_{k}$, provided we have the factorization (4.2), for cost $\mathcal{O}(p)$
for each $t$. This does not affect the overall conclusion of Table 1.

(a) DFBGN with $p=n/10$

(b) DFBGN with $p=n$

Figure 1: Performance profiles (in evaluations) for DFBGN when $p=n/10$ and
$p=n$, comparing removing points with the geometry-aware Algorithm 4 and
without Lagrange polynomials (by distance to the current iterate). We use
accuracy level $\tau=10^{-5}$, and results are an average of 10 runs, each
with a budget of $100(n+1)$ evaluations. The problem collection is (CR). See
Section 5.1 for details on the testing framework.

##### Alternative Point Removal Mechanism

Instead of Algorithm 4, we could have used a simpler mechanism for removing
points, such as removing the points furthest from the current iterate (with
total cost999As above, if we have (4.2), we could calculate all distances to
$\bm{x}_{k+1}$ using columns of $R_{k}$, with total cost $\mathcal{O}(p^{2})$.
$\mathcal{O}(np)$). However, this leads to a substantial performance penalty.
In Figure 1, we compare these two approaches for selecting points to be
removed, namely Algorithm 4 and distance to $\bm{x}_{k+1}$, on the (CR) test
set with $p=n/10$ and $p=n$ (using the default value of $p_{\rm drop}$, as
detailed in Section 4.4). For more details on the numerical testing framework,
see Section 5.1 below. We see that the geometry-aware criterion (4.5) gives
substantially better performance than the cheaper criterion.

####  4.3.2 Generation of New Directions

We now detail how new directions $\bm{d}_{1},\ldots,\bm{d}_{q}$ are created in
line 17 of DFBGN (Algorithm 3). The same approach is suitable for generating
the initial directions $\bm{s}_{1},\ldots,\bm{s}_{p}$ in line 3 of DFBGN,
using $\widetilde{A}=A$ below (i.e. no $Q$ required).

1:Orthonormal basis for current subspace $Q\in\mathbb{R}^{n\times p_{1}}$
(optional), number of new directions $q\leq n-p_{1}$.

2:Generate $A\in\mathbb{R}^{n\times q}$ with i.i.d. standard normal entries.

3:If $Q$ is specified, calculate $\widetilde{A}=A-QQ^{T}A$, otherwise set
$\widetilde{A}=A$.

4:Perform the QR factorization $\widetilde{Q}\widetilde{R}=\widetilde{A}$ and
return $\bm{d}_{1},\ldots,\bm{d}_{q}$ as the columns of $\widetilde{Q}$.

Algorithm 5 Mechanism for generating new directions in DFBGN.

Suppose our current subspace is defined by the orthonormal columns of
$Q\in\mathbb{R}^{n\times p_{1}}$, and we wish to generate $q$ new orthonormal
vectors that are also orthogonal to the columns of $Q$ (with $p_{1}+q\leq n$).
When called in line 17 of DFBGN, we will have $p_{1}=p-p_{\rm drop}$ and
$q=p_{\rm drop}$. We use the approach in Algorithm 5. From the QR
factorization, the columns of $\widetilde{Q}$ are orthonormal, and if
$\widetilde{A}$ is full rank (which occurs with probability 1; see Lemma 4.5
below) then we also have
$\operatorname{col}(\widetilde{Q})=\operatorname{col}(\widetilde{A})$. So, to
confirm the columns of $\widetilde{Q}$ are orthogonal to $Q$, we only need to
check that the columns of $\widetilde{A}$ are orthogonal to $Q$. Let
$\bm{\widetilde{a}}_{i}$ be the $i$-th column of $\widetilde{A}$ and
$\bm{q}_{j}$ be the $j$-th column of $Q$. Then, if $\bm{a}_{i}$ is the $i$-th
column of $A$, we have

$\bm{\widetilde{a}}_{i}^{T}\bm{q}_{j}=\bm{a}_{i}^{T}(I-QQ^{T})\bm{q}_{j}=\bm{a}_{i}^{T}(\bm{q}_{j}-Q\bm{e}_{j})=0,$
(4.6)

as required.

The cost of Algorithm 5 is $\mathcal{O}(nq)$ to generate $A$,
$\mathcal{O}(np_{1}q)$ to form $\widetilde{A}$ and $\mathcal{O}(nq^{2})$ for
the QR factorization. Since $p_{1},q\leq p$ (since $p_{1}$ is the number of
directions remaining in the interpolation set and $q$ is the number of new
directions to be added), the whole process has cost at most
$\mathcal{O}(np^{2})$. This bound is tight, up to constant factors, as we
could take $p_{1}=q=p/2$, for instance.

######  Lemma 4.5.

The matrix $\widetilde{A}$ has full column rank with probability 1.

###### Proof.

Let $\bm{a}_{i}$ and $\bm{\widetilde{a}}_{i}$ be the $i$-th columns of $A$ and
$\widetilde{A}$ respectively. From [29, Proposition 7.1], $A$ has full column
rank with probability 1, and each $\bm{a}_{i}\notin\operatorname{col}(Q)$ with
probability 1. Now suppose we have constants $c_{1},\ldots,c_{q}$ so that
$\sum_{i=1}^{q}c_{i}\bm{\widetilde{a}}_{i}=\bm{0}$. Then since
$\widetilde{\bm{a}}_{i}=\bm{a}_{i}-QQ^{T}\bm{a}_{i}$, we have

$\sum_{i=1}^{q}c_{i}\bm{a}_{i}=\sum_{i=1}^{q}c_{i}QQ^{T}\bm{a}_{i}.$ (4.7)

The right-hand side is in $\operatorname{col}(Q)$, so since
$\bm{a}_{i}\notin\operatorname{col}(Q)$, we must have
$\sum_{i=1}^{q}c_{i}\bm{a}_{i}=\bm{0}$. Thus $c_{1}=\cdots=c_{q}=0$ since $A$
has full column rank, and so $\widetilde{A}$ has full column rank. ∎

###  4.4 Selecting an Appropriate Value of $\bm{p_{\rm drop}}$

An important component of DFBGN that we have not yet specified is how many
points to remove from the interpolation set at each iteration, $p_{\rm
drop}\in\\{1,\ldots,p\\}$.

On one hand, a large $p_{\rm drop}$ enables us to change the subspace by a
large amount between iterations, ensuring we explore the whole of
$\mathbb{R}^{n}$ quickly, rather than searching in unproductive subspaces for
many iterations. However, a small $p_{\rm drop}$ means we require few
objective evaluations per iteration, and so are more likely to use our
evaluation budget efficiently. We consider two choices of $p_{\rm drop}$,
aimed at each of these possible benefits: $p_{\rm drop}=p/10$ to change
subspaces quickly, and $p_{\rm drop}=1$ (the minimum possible value) to use
few objective evaluations.

Another approach that we consider is a compromise between these two choices.
We note that having $p_{\rm drop}=1$ is useful to make progress with few
evaluations, so we use this value while we are making progress—we consider
this to occur when we have a successful iteration (i.e.
$\rho_{k}\geq\eta_{1}$). When we are not progressing (i.e. unsuccessful steps
with $\rho_{k}<\eta_{1}$), we use the larger value $p_{\rm drop}=p/10$.

(a) DFBGN with $p=n/2$

(b) DFBGN with $p=n$

Figure 2: Performance profiles (in evaluations) comparing different choices of
$p_{\rm drop}$, for DFBGN with $p=n/2$ and $p=n$, with accuracy level
$\tau=10^{-5}$. The choice ‘$p_{\rm drop}$ mixed’ uses $p_{\rm drop}=1$ for
successful iterations and $p_{\rm drop}=p/10$ for unsuccessful iterations.
Results an average of 10 runs, each with a budget of $100(n+1)$ evaluations.
The problem collection is (CR). See Section 5.1 for details on the testing
framework.

We compare these three approaches on the (CR) problem collection with a budget
of $100(n+1)$ evaluations in Figure 2. Since these different choices of
$p_{\rm drop}$ are similar when $p$ is small, we show results for subspace
dimensions $p=n/2$ and $p=n$. We first see that, even with $p=n/2$, the three
approaches all perform similarly. However, for $p=n$ the compromise choice
$p_{\rm drop}\in\\{1,p/10\\}$ performs better than the two constant-$p$
approaches. In addition, $p_{\rm drop}=1$ outperforms $p_{\rm drop}=p/10$ for
small performance ratios, but is less robust and solves fewer problems
overall.

Given these results, in DFBGN we use the compromise choice as the default
mechanism: $p_{\rm drop}=1$ on successful iterations and $p_{\rm drop}=p/10$
on unsuccessful iterations.

##### Relationship to model-improvement phases

The CHECK_MODEL flag in RSDFO-GN is important for ensuring we do not reduce
$\Delta_{k}$ too quickly without first ensuring the quality of the
interpolation model.101010This is related to ensuring $\Delta_{k}$ does not
get too small compared to $\|\hat{\bm{g}}_{k}\|$ via Lemma 2.6. For a similar
purpose, DFO-LS incorporates a second trust-region radius which also is
involved with ensuring $\Delta_{k}$ does not decrease too quickly [14]. In
DFBGN, as described in Section 4.3.1, we maintain the geometry of the
interpolation set by replacing poorly-located points with orthogonal
directions around the current iterate; in practice this ensures the quality of
the interpolation set. However, the choice of $p_{\rm drop}$ has a large
impact on causing $\Delta_{k}$ to shrink too quickly.

In many cases, DFBGN may reach a point where its model is not accurate and we
start to have unsuccessful iterations. To fix this (and continue making
progress), we need to introduce several new interpolation points to produce a
high-quality model. If $p_{\rm drop}$ is small, this may take many
unsuccessful iterations, causing $\Delta_{k}$ to decrease quickly.

(a) drcavty1, $p_{\rm drop}=1$

(b) luksan13, $p_{\rm drop}=1$

Figure 3: Comparison of $\Delta_{k}$, $\|\hat{\bm{g}}_{k}\|$ (left $y$-axis)
and normalized objective value (right $y$-axis) for DFBGN with $p=n$ and
$p_{\rm drop}=1$. The two problems are taken from the (CR) collection. See
Section 5.1 for details on the testing framework.

The result of having $p_{\rm drop}$ small is seen in Figure 3. Here, we show
$\Delta_{k}$, $\|\hat{\bm{g}}_{k}\|$ and $f(\bm{x}_{k})$ for DFBGN with $p=n$
and $p_{\rm drop}=1$ for two problems from the (CR) collection. Both problems
show that $\Delta_{k}$ can quickly shrink to be much smaller than
$\|\hat{\bm{g}}_{k}\|$ before reaching optimality. In the case of drcavty1, we
see multiple instances where, after several unsuccessful iterations, we
recover a high-quality model and continue making progress (causing
$\Delta_{k}$ to increase again); this manifests itself as large oscillations
in $\Delta_{k}$ with comparatively little change in $\|\hat{\bm{g}}_{k}\|$.
Ultimately, as we terminate on $\Delta_{k}\leq\Delta_{\rm end}=10^{-8}$, DFBGN
quits without solving the problem (reaching accuracy $\tau\approx 6\times
10^{-3}$). A more extreme version of this behaviour is seen for problem
luksan13, where we terminate on small $\Delta_{k}$ in the first sequence of
unsuccessful iterations—DFBGN does not allow enough time to recover a high-
quality model and terminates after achieving accuracy $\tau\approx 0.3$.

(a) drcavty1, $p_{\rm drop}$ mixed

(b) luksan13, $p_{\rm drop}$ mixed

Figure 4: Comparison of $\Delta_{k}$, $\|\hat{\bm{g}}_{k}\|$ (left $y$-axis)
and normalized objective value (right $y$-axis) for DFBGN with $p=n$ and the
default $p_{\rm drop}\in\\{1,p/10\\}$. The two problems are taken from the
(CR) collection. See Section 5.1 for details on the testing framework.

This effect is mitigated by our default choice of $p_{\rm
drop}\in\\{1,p/10\\}$. By using a larger $p_{\rm drop}$ on unsuccessful
iterations, when our model is performing poorly, our interpolation set is
changed quickly. This results in DFBGN recovering a high-quality model after a
smaller decrease in $\Delta_{k}$. To demonstrate this, in Figure 4 we show the
results of DFBGN with this $p_{\rm drop}$ for the same problems as Figure 3
above. In both cases, we still see oscillations in $\Delta_{k}$, but their
magnitude is substantially reduced—it takes fewer iterations to get successful
steps, and $\Delta_{k}$ stays well above $\Delta_{\rm end}$. This leads to
both problems being solved to high accuracy.

In Figure 4, we also see that, as we approach the solution,
$\|\hat{\bm{g}}_{k}\|$ and $\Delta_{k}$ decrease at the same rate, as we would
hope. For drcavty1 after iteration 150, we also see the phenomenon described
above, where $\Delta_{k}$ can become much larger than $\|\hat{\bm{g}}_{k}\|$
due to many successful iterations, before an unsuccessful iteration with
$\|\bm{s}_{k}\|$ small means that $\Delta_{k}$ returns to the level of
$\|\hat{\bm{g}}_{k}\|$ regularly.

##### Alternative Mechanism for Recovering High-Quality Models

An alternative approach for avoiding unnecessary decreases in $\Delta_{k}$
while the interpolation model quality is improved is to simply decrease
$\Delta_{k}$ more slowly on unsuccessful iterations. This corresponds to
setting $\gamma_{\rm dec}$ to be closer to 1, which is the default choice of
DFO-LS for noisy problems (see [14, Section 3.1]), and aligns with our
theoretical requirements on the trust-region parameters (Theorem 2.18).

In Figure 5, we compare the DFBGN default choices, of $p_{\rm
drop}\in\\{1,p/10\\}$ and $\gamma_{\rm dec}=0.5$, with $p_{\rm drop}=1$ and
$\gamma_{\rm dec}\in\\{0.5,0.98\\}$ on the (CR) problem collection. For small
values of $p$ (where the different choices of $p_{\rm drop}$ are essentially
identical), the choice of $\gamma_{\rm dec}$ has almost no impact on the
performance of DFBGN. For larger values of $p$, using $\gamma_{\rm dec}=0.98$
with $p_{\rm drop}=1$ performs comparably well to the DFBGN default
($\gamma_{\rm dec}=0.5$ with $p_{\rm drop}\in\\{1,p/10\\}$). However, we opt
for keeping $\gamma_{\rm dec}=0.5$, to allow us to use the larger value for
noisy problems (just as in DFO-LS), and to reduce the risk of overfitting our
trust-region parameters to a particular problem collection.

(a) DFBGN with $p=n/100$

(b) DFBGN with $p=n/10$

(c) DFBGN with $p=n/2$

(d) DFBGN with $p=n$

Figure 5: Performance profiles (in evaluations) comparing different choices of
$p_{\rm drop}$ and $\gamma_{\rm dec}$ for DFBGN, at accuracy level
$\tau=10^{-5}$. The choice ‘$p_{\rm drop}$ mixed’ uses $p_{\rm drop}=1$ for
successful iterations and $p_{\rm drop}=p/10$ for unsuccessful iterations.
Results an average of 10 runs, each with a budget of $100(n+1)$ evaluations.
The problem collection is (CR). See Section 5.1 for details on the testing
framework.

##  5 Numerical Results

In this section we compare the performance of DFBGN (Algorithm 3) to that of
DFO-LS. We note that that DFO-LS has been shown to have state-of-the-art
performance compared to other solvers in [14]. As described in Section 4.2,
the implementation of DFBGN is based on the decision to reduce the linear
algebra cost of the algorithm at the expense of more objective evaluations per
iteration. However, we still maintain the goal of DFBGN achieving (close to)
state-of-the-art performance when it is run as a ‘full space’ method (i.e.
$p=n$). Here, we will investigate this tradeoff in practice.

###  5.1 Testing Framework

In our testing, we will compare a Python implementation of DFBGN (Algorithm 3)
against DFO-LS version 1.0.2 (also implemented in Python). The implementation
of DFBGN is available on Github.111111See
https://github.com/numericalalgorithmsgroup/dfbgn. Results here use version
0.1. We will consider both the standard version of DFO-LS, and one where we
use a reduced initialization cost of $n/100$ evaluations (c.f. Remark 4.1).
This will allow us to compare both the overall performance of DFBGN and its
performance with small budgets (since DFBGN also has a reduced initialization
cost of $p+1$ evaluations). We compare these against DFBGN with the choices
$p\in\\{n/100,n/10,n/2,n\\}$ and the adaptive choice of $p_{\rm
drop}\in\\{1,p/10\\}$ (Section 4.4). All default settings are used for both
solvers, and since both are randomized (DFO-LS uses random initial directions
only, and DFBGN is randomized through Algorithm 5), we run 10 instances of
each problem under all solver configurations.

##### Test Problems

We will consider two collections of nonlinear least-squares test problems,
both taken from the CUTEst collection [38]. The first, denoted (CR), is a
collection of 60 medium-scale problems (with $25\leq n\leq 110$ and $n\leq
m\leq 400$). Full details of the (CR) collection may be found in [18, Table
3]. The second, denoted (CR-large), is a collection of 28 large-scale problems
(with $1000\leq n\leq 5000$ and $n\leq m\leq 9998$). This collection is a
subset of problems from (CR), with their dimension increased substantially.
Full details of the (CR-large) collection are given in Appendix B. Note that
the 12 hour runtime limit was only relevant for (CR-large) in all cases.

##### Measuring Solver Performance

For every problem, we allow all solvers a budget of at most $100(n+1)$
objective evaluations (i.e. evaluations of the full vector $\bm{r}(\bm{x})$).
This dimension-dependent choice may be understood as equivalent to 100
evaluations of $\bm{r}(\bm{x})$ and the Jacobian $J(\bm{x})$ via finite
differencing. However, given the importance of linear algebra cost to our
comparisons, we allow each solver a maximum runtime of 12 hours for each
instance of each problem.121212Since all problems are implemented in Fortran
via CUTEst, the cost of objective evaluations for this testing is minimal. For
each solver $S$, each problem instance $P$, and accuracy level $\tau\in(0,1)$,
we calculate

$\displaystyle N(S,P,\tau)$ $\displaystyle:=\text{\\# evaluations of
$\bm{r}(\bm{x})$ required to find a point $\bm{x}$ with}$
$\displaystyle\qquad\qquad f(\bm{x})\leq
f(\bm{x}^{*})+\tau(f(\bm{x}_{0})-f(\bm{x}^{*})),$ (5.1)

where $f(\bm{x}^{*})$ is an estimate of the minimum of $f$ as listed in [18,
Table 3] for (CR) and Appendix B for (CR-large). If this objective decrease is
not achieved by a solver before its budget or runtime limit is hit, we set
$N(S,P,\tau)=\infty$. We then compare solver performances on a problem
collection $\mathcal{P}$ by plotting either data profiles [55]