# A matrix algebra approach to approximate Hessians

Warren Hare Department of Mathematics, University of British Columbia,
Okanagan Campus, Kelowna, B.C. V1V 1V7, Canada. Research partially supported
by NSERC of Canada Discovery Grant 2018-03865<EMAIL_ADDRESS>ORCID
0000-0002-4240-3903 Gabriel Jarry–Bolduc Department of Mathematics, University
of British Columbia, Okanagan Campus, Kelowna, B.C. V1V 1V7, Canada. Research
partially supported by Natural Sciences and Engineering Research Council
(NSERC) of Canada Discovery Grant 2018-03865<EMAIL_ADDRESS>Chayne
Planiden School of Mathematics and Applied Statistics, University of
Wollongong, Wollongong, NSW, 2500, Australia. Research supported by University
of Wollongong<EMAIL_ADDRESS>ORCID 0000-0002-0412-8445

###### Abstract

This work presents a novel matrix-based method for constructing an
approximation Hessian using only function evaluations. The method requires
less computational power than interpolation-based methods and is easy to
implement in matrix-based programming languages such as MATLAB. As only
function evaluations are required, the method is suitable for use in
derivative-free algorithms.

For reasonably structured sample sets, the method is proven to create an
order-$1$ accurate approximation of the full Hessian. Under more specialized
structures, the method is proved to yield order-$2$ accuracy. The undetermined
case, where the number of sample points is less than required for full
interpolation, is studied and error bounds are developed for the resulting
partial Hessians.

Keywords: Hessian approximation; derivative-free optimization; order-$N$
accuracy; generalized simplex gradient;generalized simplex Hessian;
generalized centered simplex Hessian.

AMS subject classification: primary 65K10, 90C56.

##  1 Introduction

Derivative-free optimization (DFO) is the study of finding the minimum value
and minimizers of a function without using gradients or higher-order
derivative information [3]. DFO is gaining in popularity in recent years and
is useful in cases where gradients are unavailable, computationally expensive,
or simply difficult to obtain [3, 17]. DFO methods have been used in both
smooth [9, 14, 22, 34, 39, 40, 43, 45] and nonsmooth [4, 7, 29, 31, 33, 36]
optimization, and most commonly take the form of either direct-search [1, 2,
6, 8, 23] or model-based [30, 34, 35, 44, 45] methods, while some use a blend
of both [5, 20]. Model-based DFO methods use numerical analysis techniques to
approximate gradients and Hessians in a manner that has controllable error
bounds. For example, a technique to approximate gradients is the simplex
gradient. This method consists of taking a set of $n+1$ sample points in
$\operatorname{\mathbb{R}}^{n}$ that are appropriately spaced, called a
simplex, and using them to build a linear interpolation function that
approximates the objective function locally, then calculating the gradient of
the affine function. The simplex gradient has been shown to have a nicely
bounded error, for functions whose gradients are locally Lipschitz continuous
[10].

One recent line of research explores methods of improving or generalizing the
simplex gradient for its use in DFO [15, 16, 18, 21, 25, 28, 38, 41]. Working
from ideas in [16], in [25] the simplex gradient was generalized so as not to
require exactly $n+1$ points; an error-controlled approximation can now be
found using any finite number of properly-spaced points. In [28], a similar
approach was used to examine the generalized centered simplex gradient, an
approximation that uses twice as many points, but results in an improvement on
the error control from order-$1$ to order-$2$ [24].

In a similar vein, researchers have also explored methods to approximate full
Hessians or partial Hessians. In [21], the authors outline an idea for a
“simplex Hessian” that is constructed via quadratic interpolation through
$(n+1)(n+2)/2$ well-poised sample points. They further suggested that if only
a portion of the Hessian were desired (say the diagonal component), then fewer
points could be used. These ideas were formalized in [16] through quadratic
interpolation and analyzed through the use of Lagrange interpolating
polynomials. Obtaining an approximation of the diagonal component of a Hessian
is also discussed in [19, 32].

In this work, we continue the development of these tools by introducing the
_generalized simplex Hessian (GSH)_ and the _generalized centered simplex
Hessian (GCSH)_. The GSH is closely linked to the Hessian of a quadratic
interpolation function. We will see that in certain situations, both
approaches yield the same result.

This paper can be viewed as an extension of the work related to simplex
Hessians introduced in [17, 21]. We introduce an explicit formula based on
matrix algebra concepts to approximate the Hessian of a function $f$. This
formula provides several advantages. First, the formula provides an approach
that is well-defined regardless of the number of sample points utilized.
Indeed, as long as the matrices used to build the sample set of points are
nonempty, the GSH and the GCSH are well-defined. In particular, the GSH
provides an explicit formula to approximate the Hessian even when the
quadratic interpolation function does not exist or is not unique. However, we
note that when the matrices used in the computation of the GSH have a specific
structure, the GSH is equivalent to a _forward-finite-difference
approximation_ of the Hessian.

The GSH provides an accurate approximation of the Hessian of an objective
function $f:\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$ under
reasonable assumptions. The technique uses a point of interest $x^{0}$ and
sets of directions, stored in a set of matrices $\\{S,T_{1},\ldots,T_{m}\\}$
whose columns are vectors that are added to the point of interest to build the
sample points. Defining $\Delta_{S}$ as the radius of $S$ and $\Delta_{T}$ as
the maximum radius of the $T_{i}$, we prove that if $S$ and all $T_{i}$ are
full row rank and the Hessian $\nabla^{2}f$ exists and is Lipschitz
continuous, then the GSH is an accurate estimate of $\nabla^{2}f$ to within a
multiple of $\Delta_{u}:=\max\\{\Delta_{S},\Delta_{T}\\}$. In terms of
order-$N$ accuracy [24], the GSH provides order-$1$ accurate approximations of
the full Hessian (see Definition 2.11 herein). Furthermore, the GCSH provides
order-2 accurate approximations of the full Hessian. Error bounds for the
_underdetermined case_ (see Definition 2.9) are presented showing that the GSH
and its centered version, the GCSH, are order-1 and order-2 accurate
approximation of the appropriate partial Hessian, respectively.

The error bounds presented in this paper share similarities to those presented
in [15, 16]. In those papers, the authors provide error bounds for the Hessian
of a fully quadratic model [15, Theorem 3], the Hessian of a quadratic
regression model [16, Theorem 3.2], and the Hessian of a quadratic
undetermined model [16, Theorem 5.12].

A second benefit of the GSH formula is that by employing a matrix structure in
the definition, the Hessian approximation is extremely easy to implement in
any matrix-based programming language. Moreover, after function evaluations
are complete, the matrix structure requires less computational effort (in
flops) than quadratic interpolation. Indeed, when $S$ and all $T_{i}$ are
square matrices, then the GSH requires order $O(n^{4})$ flops to compute. This
is an improvement to the order $O(n^{6})$ flops required to find the Hessian
of the quadratic interpolation function using a set of sample points that is
poised for quadratic interpolation [15] (see also [17, §6.2]).

The remainder of this paper is organized as follows. Section 2 contains a
description of notation and some needed definitions, including those of the
generalized simplex gradient, the GSH and the GCSH. Section 3 clarifies the
relation between the GSH and the GCSH. Section 4 presents several error bounds
for the GSH and the GCSH. The error bounds are divided into two categories:
the situation where all matrices of directions $T_{i}$ are equal, and the
situation where all matrices $T_{i}$ are not necessarily equal. For boths
situations, the error bound is valid for the underdetermined case. In Section
5, we establish the minimal poised set for GSH computation and show that it is
well-poised for quadratic interpolation. Finally, Section 6 contains
concluding remarks and recommends areas of future research in this vein.

##  2 Preliminaries

Throughout this work, we use the standard notation found in [42]. The domain
of a function $f$ is denoted by $\operatorname{dom}f$. The transpose of a
matrix $A$ is denoted by $A^{\top}$. We work in finite-dimensional space
$\operatorname{\mathbb{R}}^{n}$ with inner product
$x^{\top}y=\sum_{i=1}^{n}x_{i}y_{i}$ and induced norm
$\|x\|=\sqrt{x^{\top}x}$. The identity matrix in
$\operatorname{\mathbb{R}}^{n\times n}$ is denoted by $\operatorname{Id}_{n}$.
We use $e^{i}_{n}\in\operatorname{\mathbb{R}}^{n}$ for $i\in\\{1,\dots,n\\}$
to denote the standard unit basis vectors in $\operatorname{\mathbb{R}}^{n}$,
i.e. the $i$th column of $\operatorname{Id}_{n}.$ The matrix
$\operatorname{Diag}(v)=\operatorname{Diag}(v_{1},\dots,v_{n})=:M\in\operatorname{\mathbb{R}}^{n\times
n},$ where $v\in\operatorname{\mathbb{R}}^{n},$ represents the diagonal matrix
with diagonal entries $M_{i,i}=v_{i}$ for all $i\in\\{1,\dots,n\\}.$ The zero
vector in $\operatorname{\mathbb{R}}^{n}$ is denoted $\mathbf{0}$ and the zero
matrix in $\operatorname{\mathbb{R}}^{n\times m}$ is denoted
$\mathbf{0}_{n\times m}$. The entry in the $i$th row and $j$th column of a
matrix $A$ is denoted $A_{i,j}.$

Given a matrix $A\in\operatorname{\mathbb{R}}^{n\times m},$ we use the induced
matrix norm

$\displaystyle\|A\|=\|A\|_{2}=\max\\{\|Ax\|_{2}\,:\,\|x\|_{2}=1\\}$

and the Frobenius norm

$\displaystyle\|A\|_{F}=\left(\sum_{i=1}^{n}\sum_{j=1}^{m}A_{i,j}^{2}\right)^{\frac{1}{2}}.$

We denote by $B(x^{0},\Delta)$ and $\overline{B}(x^{0},\Delta)$ the open and
closed balls, respectively, centered about $x^{0}$ with radius $\Delta$. The
_Minkowski sum_ of two sets of vectors $A$ and $B$ is denoted by $A\oplus B.$
That is

$\displaystyle A\oplus B=\\{a+b:a\in A,b\in B\\}.$

We define a quadratic function
$Q:\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$ to be a function
of the form $Q(x)=\alpha_{0}+\alpha^{\top}x+\frac{1}{2}x^{\top}\mathbf{H}x$
where
$\alpha_{0}\in\operatorname{\mathbb{R}},\alpha\in\operatorname{\mathbb{R}}^{n}$
and $\mathbf{H}=\mathbf{H}^{\top}\in\operatorname{\mathbb{R}}^{n\times n}.$ An
affine function
$\mathcal{L}:\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$ is
defined to be any function that can be written in the form
$\mathcal{L}(x)=\alpha_{0}+\alpha^{\top}x$. Note, affine functions and
constant functions $C(x)=\alpha_{0}$ are also considered quadratic functions,
with $\mathbf{H}=\mathbf{0}_{n\times n}$.

Next, we introduce fundamental definitions and notation that will be used
throughout this paper.

######  Definition 2.1 (Poised for quadratic interpolation).

_[15, 16, 17]_ The set of distinct points
$\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{m}\\}\subset\operatorname{\mathbb{R}}^{n}$
with $m=\frac{1}{2}(n+1)(n+2)-1$ is poised for quadratic interpolation if the
system

$\displaystyle\alpha_{0}+\alpha^{\top}y^{i}+\frac{1}{2}(y^{i})^{\top}\mathbf{H}y^{i}=0,\quad
i\in\\{0,1,\dots,m\\},$ (1)

has a unique solution for
$\alpha_{0}\in\operatorname{\mathbb{R}},\alpha\in\operatorname{\mathbb{R}}^{n},$
and $\mathbf{H}=\mathbf{H}^{\top}\in\operatorname{\mathbb{R}}^{n\times n}.$

When a set of sample points $\mathcal{Y}$ is poised for quadratic
interpolation, it means that there exists a unique quadratic model passing
through all the sample points in $\mathcal{Y}.$ Therefore, the Hessian of the
quadratic model defined through $\mathcal{Y}$ is unique. Note that a set of
sample points containing $(n+1)(n+2)/2$ distinct sample points is not
necessarily poised for quadratic interpolation. The set of sample points is
poised for quadratic interpolation if and only the matrix associated to the
linear system (1) is square and full rank. For instance, the set of sample
points

$\left\\{\begin{bmatrix}1\\\ 0\end{bmatrix},\begin{bmatrix}0\\\
1\end{bmatrix},\begin{bmatrix}-1\\\ 0\end{bmatrix},\begin{bmatrix}0\\\
-1\end{bmatrix},\begin{bmatrix}\frac{1}{\sqrt{2}}\\\
\frac{1}{\sqrt{2}}\end{bmatrix},\begin{bmatrix}-\frac{1}{\sqrt{2}}\\\
-\frac{1}{\sqrt{2}}\end{bmatrix}\right\\}\subset\operatorname{\mathbb{R}}^{2}$

is not poised for quadratic interpolation since the rank of the matrix
associated to the system (1) is $5<6$. The set

$\left\\{\begin{bmatrix}0\\\ 0\end{bmatrix},\begin{bmatrix}1\\\
0\end{bmatrix},\begin{bmatrix}0\\\ 1\end{bmatrix},\begin{bmatrix}-1\\\
0\end{bmatrix},\begin{bmatrix}0\\\
-1\end{bmatrix},\begin{bmatrix}\frac{1}{\sqrt{2}}\\\
-\frac{1}{\sqrt{2}}\end{bmatrix}\right\\}\subset\operatorname{\mathbb{R}}^{2}$

is poised for quadratic interpolation since the rank of the matrix associated
to the system (1) is 6. When the number of distinct sample points in
$\mathcal{Y}$ is fewer than $(n+1)(n+2)/2,$ then the quadratic interpolation
model is no longer unique. This case is referred as the underdetermined case
[17, Chapter 5]. In this case, a quadratic model may be defined through the
minimum Frobenius norm problem [17, Chapter 5] (see also [16, Section 5]).
When the number of distinct sample points in $\mathcal{Y}$ is greater than
$(n+1)(n+2)/2$ and there exists no quadratic function that passes through all
points in $\mathcal{Y},$ the quadratic model may be determined through a least
square regression problem [17, Chapter 4] (see also [16, Sections 2,3 & 4]).

######  Definition 2.2 (Quadratic interpolation function).

_[3, Definition 9.9]_ Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{m}\\}\subset\operatorname{dom}f$
with $m=\frac{1}{2}(n+1)(n+2)-1$ be poised for quadratic interpolation. Then
the quadratic interpolation function of $f$ over $\mathcal{Y}$ is

$\displaystyle Q_{f}(\operatorname{\mathcal{Y}})(x)$
$\displaystyle=\alpha_{0}+\alpha^{\top}x+\frac{1}{2}x^{\top}\mathbf{H}x,$

where $(\alpha_{0},\alpha,\mathbf{H}=\mathbf{H}^{\top})$ is the unique
solution to

$\displaystyle\alpha_{0}+\alpha^{\top}y^{i}+\frac{1}{2}(y^{i})^{\top}\mathbf{H}y^{i}$
$\displaystyle=f(y^{i}),\quad i\in\\{0,1,\dots,m\\}.$

In the next definition, we introduce key notation used in the construction of
the GSH. Within, we write a set of vectors in matrix form, by which we mean
that each column of the matrix is a vector in the set.

######  Definition 2.3 (GSH notation).

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the _point of interest_. Let

$\displaystyle S$
$\displaystyle=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}~{}~{}\mbox{and}$ $\displaystyle T_{i}$
$\displaystyle=\begin{bmatrix}t^{1}_{i}&t^{2}_{i}&\cdots&t^{k}_{i}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
k_{i}},i\in\\{1,\ldots,m\\}$

be sets of directions contained in $\operatorname{\mathbb{R}}^{n}$, written in
matrix form. Define

$T_{1:m}=\\{T_{1},\ldots,T_{m}\\}.$

Assume that $x^{0}\oplus T_{i},x^{0}+s^{i},x^{0}+s^{i}\oplus T_{i}$ are
contained in $\operatorname{dom}f$ for all $i~{}\in\\{1,\ldots,m\\}$. Define

$\Delta_{S}=\max\limits_{i\,\in\\{1,\ldots,m\\}}\|s^{i}\|,\quad\Delta_{T_{i}}=\max\limits_{j\,\in\\{1,\ldots,k_{i}\\}}\|t^{j}_{i}\|,\quad\Delta_{T}=\max\limits_{i\in\\{1,\ldots,m\\}}\Delta_{T_{i}},$
$\displaystyle\widehat{S}=\frac{1}{\Delta_{S}}S,$ (2)

and

$\displaystyle\delta_{s}\left(x^{0};T_{i}\right)$
$\displaystyle=\left[\begin{array}[]{c}f(x^{0}+t^{1}_{i})-f(x^{0})\\\
f(x^{0}+t^{2}_{i})-f(x^{0})\\\ \vdots\\\
f(x^{0}+t^{k}_{i})-f(x^{0})\end{array}\right]\in\operatorname{\mathbb{R}}^{k_{i}}.$

In this paper, all matrices involved are assumed to have non-null rank. This
ensures that the radius of the matrix is nonzero and hence, (2) is always
well-defined. Note that in the previous definition,
$m\in\operatorname{\mathbb{N}}$ and $k_{i}$ can be any positive integer for
all $i\in\\{1,\dots,m\\}$.

Recall that for nonsquare matrices, a generalization of the matrix inverse is
the pseudoinverse. The most well-known type of matrix pseudoinverse is the
Moore–Penrose pseudoinverse.

######  Definition 2.4 (Moore–Penrose pseudoinverse).

Let $A\in\operatorname{\mathbb{R}}^{n\times m}$. The _Moore–Penrose
pseudoinverse_ of $A$, denoted by $A^{\dagger}$, is the unique matrix in
$\operatorname{\mathbb{R}}^{m\times n}$ that satisfies the following four
equations:

  1. (i)

$AA^{\dagger}A=A,$

  2. (ii)

$A^{\dagger}AA^{\dagger}=A^{\dagger},$

  3. (iii)

$(AA^{\dagger})^{\top}=AA^{\dagger},$

  4. (iv)

$(A^{\dagger}A)^{\top}=A^{\dagger}A.$

Note that given $A\in\operatorname{\mathbb{R}}^{n\times m},$ there exists a
unique Moore–Penrose pseudoinverse
$A^{\dagger}\in\operatorname{\mathbb{R}}^{m\times n}.$ The following two
properties hold.

  1. (i)

If $A$ has full column rank $m$, then $A^{\dagger}$ is a left-inverse of $A$,
that is $A^{\dagger}A=\operatorname{Id}_{m}.$ In this case,

$A^{\dagger}=(A^{\top}A)^{-1}A^{\top}.$ (3)

  2. (ii)

If $A$ has full row rank $n$, then $A^{\dagger}$ is a right-inverse of $A$,
that is $AA^{\dagger}=\operatorname{Id}_{n}.$ In this case,
$A^{\dagger}=A^{\top}(AA^{\top})^{-1}.$

Before introducing the GSH, it is valuable to recall the definition of the
generalized simplex gradient.

######  Definition 2.5 (Generalized simplex gradient).

_[25, Definition 2]_ Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S\in\operatorname{\mathbb{R}}^{n\times m}$ with $x^{0}\oplus
S\subset\operatorname{dom}f.$ The _generalized simplex gradient_ of $f$ at
$x^{0}$ over $S$ is denoted by $\nabla_{s}f(x^{0};S)$ and defined by

$\nabla_{s}f(x^{0};S)=(S^{\top})^{\dagger}\delta_{s}(x^{0};S).$

The Moore–Penrose pseudoinverse allows for $S$ to contain any finite number of
columns, rather than being restricted to $m=n$. When
$S\in\operatorname{\mathbb{R}}^{n\times n}$ and $\operatorname{rank}S=n$, then
the Moore–Penrose pseudoinverse of $S^{\top}$ is the inverse of $S^{\top}$ and
we recover the definition of the simplex gradient [3, Definition 9.5]. Note
that the generalized simplex gradient is a generalization of the familiar
_forward-finite-difference approximation_ of the gradient [12, Section 4.1].
Indeed, if $S=H\operatorname{Id}_{n},$ where
$H=\operatorname{Diag}\begin{bmatrix}h_{1}&h_{2}&\dots&h_{n}\end{bmatrix},h_{i}>0,$
then we recover the forward-finite-difference approximation of the gradient.

We now introduce the two key definitions for this paper, the generalized
simplex Hessian and the generalized centered simplex Hessian. It requires
finite sets of directions and a point of interest at which the Hessian is
approximated. Similarly to the generalized simplex gradient, the GSH involves
the Moore–Penrose pseudoinverse and a difference matrix. In this case, the
difference matrix consists of the difference between generalized simplex
gradients.

For ease of notation, we use the definition $T_{1:m}=\\{T_{1},\ldots,T_{m}\\}$
where convenient.

######  Definition 2.6 (Generalized simplex Hessian).

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and $T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}$ with $x^{0}\oplus
T_{i},x^{0}\oplus S,x^{0}+s^{i}\oplus T_{i}$ contained in
$\operatorname{dom}f$ for all $i\in\\{1,\dots,m\\}.$ The _generalized simplex
Hessian_ of $f$ at $x^{0}$ over $S$ and $T_{1:m}$ is denoted by
$\nabla^{2}_{s}f(x^{0};S,T_{1:m})$ and defined by

$\nabla^{2}_{s}f(x^{0};S,T_{1:m})=(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m}),$
(4)

where

$\delta^{2}_{s}(x^{0};S;T_{1:m})=\left[\begin{array}[]{c}(\nabla_{s}f(x^{0}+s^{1};T_{1})-\nabla_{s}f(x^{0};T_{1}))^{\top}\\\
(\nabla_{s}f(x^{0}+s^{2};T_{2})-\nabla_{s}f(x^{0};T_{2}))^{\top}\\\ \vdots\\\
(\nabla_{s}f(x^{0}+s^{m};T_{m})-\nabla_{s}f(x^{0};T_{m}))^{\top}\end{array}\right]\in\operatorname{\mathbb{R}}^{m\times
n}.$

In the case $T_{1}=T_{2}=\cdots=T_{m}$, we use $\overline{T}:=T_{i}$ to
simplify notation and write $\nabla^{2}_{s}f(x^{0};S;\overline{T})$ to
emphasize the special case.

We remark that the $s$ in $\nabla^{2}_{s}f,\delta_{s}$ and $\delta^{2}_{s}$ is
for simplex. This becomes important as we now introduce the generalized
centered simplex Hessian. The ‘square’ in $\delta^{2}_{s}$ emphasizes the fact
that we are taking differences of first-order objects (hence getting a second-
order object).

######  Definition 2.7 (Generalized centered simplex Hessian).

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S\in\operatorname{\mathbb{R}}^{n\times m}$ and
$T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}$ with $x^{0}\oplus S\oplus
T_{i},x^{0}\oplus(-S)\oplus(-T_{i}),x^{0}\oplus(\pm S),$ and $x^{0}\oplus(\pm
T_{i})$ contained in $\operatorname{dom}f$ for all $i\in\\{1,\dots,m\\}.$ The
generalized centered simplex Hessian of $f$ at $x^{0}$ over $S$ and $T_{1:m}$
is denoted by $\nabla^{2}_{c}f(x^{0};S,T_{1:m})$ and defined by

$\nabla^{2}_{c}f(x^{0};S,T_{1:m})=\frac{1}{2}\left(\nabla_{s}^{2}f(x^{0};S,T_{1:m})+\nabla_{s}^{2}f(x^{0};-S,-T_{1:m})\right).$
(5)

Note that the $c$ in $\nabla^{2}_{c}f$ is for centered.

The GCSH can be viewed as a generalization of the centered-finite difference
approximation of the Hessian, also called midpoint approximation [12, Section
4.1]. Similarly, The GSH can be viewed as a generalization of the forward-
finite-difference approximation of the Hessian [12, Section 4.1]. Indeed, when
$S=T_{1}=\dots=T_{m}=H\operatorname{Id}_{n},$ where
$H=\operatorname{Diag}\begin{bmatrix}h_{1}&h_{2}&\dots&h_{n}\end{bmatrix},h_{i}>0,$
then the GSH (GCSH) is equivalent to a forward-finite-difference approximation
(centered-finite-difference approximation) of the second-order derivatives.
Moreover, the GSH can be viewed as the DFO analog of the Hessian approximation
given by Taylor’s Theorem in the derivative-based context, which approximates
the gradient by differences of true gradients (see [37, Section 8.1]). An
advantage of the GSH (GCSH) over forward-finite-difference approximations
(centered-finite-difference approximations) is that it is less restrictive. In
particular, it is well-defined as long as the matrices $S$ and
$T_{1},\dots,T_{m}$ are nonempty.

######  Remark 2.8.

Consider the case where $S$ and $T_{i}$ are square matrices, i.e., $m=n$ and
$k_{i}=n$ for all $i\in\\{1,\dots,m\\}.$ The computation of the GSH begins by
the construction of $n+1$ generalized simplex gradients. Setting aside
function evaluations, computing a simplex gradient is dominated by computing a
matrix inverse, which is $O(n^{3})$ flops. Thus the construction of
$\delta^{2}_{s}(x^{0};S;T_{1:m})$ requires $O(n^{4})$ flops. The inverse of
$S$ is computed ($O(n^{3})$ flops) and multiplied by
$\delta^{2}_{s}(x^{0};S;T_{1:m})$, $O(n^{3})$ flops. Hence, setting aside
function evaluations, computing the GSH requires $O(n^{4})$ flops.

Conversely, setting aside function evaluations, quadratic interpolation
requires $O(n^{6})$ flops to compute.

######  Definition 2.9.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and $T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}$ with $x^{0}\oplus
S\oplus T_{i},x^{0}\oplus(-S)\oplus(-T_{i}),x^{0}\oplus(\pm S),$ and
$x^{0}\oplus(\pm T_{i})$ contained in $\operatorname{dom}f$ for all
$i\in\\{1,\dots,m\\}.$ Assume that all matrices are non-null rank. We define
the following four cases to characterize the matrix
$S\in\operatorname{\mathbb{R}}^{n\times m}$ and the set of matrices
$T_{1:m}:=\\{T_{1},\dots,T_{m}\\}.$

  * •

Underdetermined: the GSH (GCSH) is said to be $S$-underdetermined if $S$ is
non-square and full column rank. We say that it is $T_{1:m}$-underdetermined
if all matrices in the set $T_{1:m}$ are full column rank and at least one
matrix is non-square.

  * •

Determined: the GSH (GCSH) is said to be $S$-determined if $S$ is square and
full rank. It is $T_{1:m}$-derdetermined if all matrices in the set $T_{1:m}$
are square and full rank.

  * •

Overdetermined: the GSH (GCSH) is said to be $S$-overdetermined if $S$ is non-
square and full row rank. It is $T_{1:m}$-overderdetermined if all matrices in
the set $T_{1:m}$ are full row rank at least one is non-square.

  * •

Nondetermined: the GSH (GCSH) is said to be $S$-nondetermined if it is not in
any of the previous three cases. It is $T_{1:m}$-nonderdetermined if the set
$T_{1:m}$ is not in any of the previous three cases.

Note that the definition of an $S$-underdetermined GSH (GCSH) implies that
$\operatorname{span}S\neq\operatorname{\mathbb{R}}^{n}$, which is true if and
only if $SS^{\dagger}\neq\operatorname{Id}_{n}$. Similarly, the definition of
a $T_{1:m}$-underdetermined GSH (GCSH) implies that
$\operatorname{span}T_{i}\neq\operatorname{\mathbb{R}}^{n}\text{ for some
}\bar{i}\in\\{1,\dots,m\\}$, which is true if and only if
$T_{i}T_{i}^{\dagger}\neq\operatorname{Id}_{n}\text{ for }\bar{i}$.

Defining underdetermined, determined, overdetermined and nondetermined as
above creates 16 different cases to classify a GSH (GCSH), all of which are
investigated in Section 4.

The following proposition presents a bound for the absolute error between the
generalized simplex gradient $\nabla_{s}f(x^{0};S)$ and the true gradient
$\nabla f(x^{0}).$ Note that the assumption that $f\in\mathcal{C}^{2}$ on
$B(x^{0},\overline{\Delta})$ implies that $\nabla f$ is Lipschitz continuous
on $\overline{B}(x^{0},\Delta)$ where $\Delta<\overline{\Delta}.$

We remark that results akin to Proposition 2.10 have been proven several times
using a variety of different techniques [21, 41, 25]. In [15, 16] (see also
[17] and [21, Theorem 3.1]), simplex gradients were analyzed using Lagrange
interpolating polynomials. In [41], (see also [3, 25]), error bounds for
simplex gradients were analyzed using matrix condition numbers. In Proposition
2.10, we present the error bound in this second format, as further results
will examine generalized simplex Hessians in the same style.

######  Proposition 2.10 (Error bound for the generalized simplex gradient).

Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ with radius $\Delta_{S}>0$. Let $\widehat{S}=S/\Delta_{S}$ and let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
be $\mathcal{C}^{2}$ on $B(x^{0},\overline{\Delta})$ where $x^{0}$ is the
point of interest and $\overline{\Delta}>\Delta_{S}.$ Denote by $L_{\nabla
f}\geq 0$ the Lipschitz constant of $\nabla f$ on
$\overline{B}(x^{0},\Delta_{S}).$ Then

$\displaystyle\|\operatorname{Proj}_{S}\nabla_{s}f(x^{0};S)-\operatorname{Proj}_{S}\nabla
f(x^{0})\|$
$\displaystyle=\|\nabla_{s}f(x^{0};S)-\operatorname{Proj}_{S}\nabla
f(x^{0})\|$ $\displaystyle\leq\frac{\sqrt{m}}{2}L_{\nabla
f}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\Delta_{S},$

where $\operatorname{Proj}_{S}(\cdot)=(S^{\top})^{\dagger}S^{\top}(\cdot).$

The error bound presented in Proposition 2.10 only requires that $S$ has non-
null rank. Therefore, it covers all possible cases: $S$-underdetermined,
$S$-determined case, $S$-overdetermined and $S$-nondetermined. In Section 4,
error bounds with a similar format than the one in Proposition 2.10 are
developed for the GSH and the GCSH. To do this, we next define the projection
onto a pair of subspaces.

Given matrices $S\in\operatorname{\mathbb{R}}^{n\times m}$ and
$T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}$, the projection of the
matrix $H\in\operatorname{\mathbb{R}}^{n\times n}$ onto $S$ and $T_{1:m}$ is
denoted by $\operatorname{Proj}_{S,T_{1:m}}H$ and defined by

$\displaystyle\operatorname{Proj}_{S,T_{1:m}}H$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}HT_{i}T_{i}^{\dagger}.$

In the case where $T_{1}=T_{2}=\dots=T_{m}=:\overline{T},$ the projection of
$H$ onto $S$ and $\overline{T}$ is denoted by
$\operatorname{Proj}_{S,\overline{T}}H,$ and reduces to

$\displaystyle\operatorname{Proj}_{S,\overline{T}}H$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}H\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}\left(\sum_{i=1}^{m}e^{i}_{m}(e^{i}_{m})^{\top}\right)S^{\top}H\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}{\operatorname{Id}}_{m}S^{\top}H\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}H\overline{T}\overline{T}^{\dagger}.$

We conclude this section by defining _order-N accuracy_ as it used in this
paper several times to describe the quality of the Hessian approximation
techniques developed.

######  Definition 2.11 (Order-$N$ accuracy).

Let $P\in\operatorname{\mathbb{R}}^{n\times n}$. Given $f\in\mathcal{C}^{2}$,
$x^{0}\in\operatorname{dom}f$ and $\overline{\Delta}>0,$ we say that
$\\{\tilde{f}_{\Delta}\\}_{\Delta\in(0,\overline{\Delta}]}$ is a class of
models of $f$ at $x^{0}$ parametrized by $\Delta$ that provides order-N
accuracy of $P$ at $x^{0}$ if there exists a scalar $\kappa(x^{0})$ such that,
given any $\Delta\in(0,\overline{\Delta}],$ the model $\tilde{f}_{\Delta}$
satisfies

$\|P-\nabla^{2}\tilde{f}_{\Delta}(x^{0})\|\leq\kappa(x^{0})\Delta^{N}.$

When defining error bounds in Section 4, the matrix $P$ in the previous
definition will be equal to
$\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$. We may refer to
$\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$ as a _partial Hessian_.
Note that when $S$ is full row rank and $T_{i}$ is full row rank for all $i,$
then $\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$ is equal to the true
Hessian $\nabla^{2}f(x^{0}).$ We will see that we can obtain order-1 accuracy,
or order-2 accuracy, of the partial Hessian
$\operatorname{Proj}_{S,T_{1:m}}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$ under the
appropriate assumptions. First, the relation between the GSH and the GCSH is
clarified in the next section.

##  3 The relation between the GSH and the GCSH

In this section, we investigate the relation between the GSH and the GCSH. We
begin by providing a formula to compute the GCSH in a format similar to the
formula for the GSH (4).

In Lemma 3.1 we introduce the notation $\delta_{c}^{2}$. This is a centered
version of $\delta^{2}_{s}$; hence, the $s$ for simplex is replaced with $c$
for centered.

######  Lemma 3.1.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and $T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}$ with $x^{0}\oplus
S\oplus T_{i},x^{0}\oplus(-S)\oplus(-T_{i}),x^{0}\oplus(\pm S),x^{0}\oplus(\pm
T_{i})$ contained in $\operatorname{dom}f$ for all $i\in\\{1,\dots,m\\}.$ Then

$\displaystyle\nabla_{c}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=(S^{\top})^{\dagger}\delta_{c}^{2}(x^{0};S,T_{1:m})$

where

$\delta_{c}^{2}(x^{0};S,T_{1:m})=\frac{1}{2}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};T_{1})+\delta_{s}(x^{0}-s^{1};-T_{1})-\delta_{s}(x^{0};T_{1})-\delta_{s}(x^{0};-T_{1})\right)^{\top}T_{1}^{\dagger}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};T_{m})+\delta_{s}(x^{0}-s^{m};-T_{m})-\delta_{s}(x^{0};T_{m})-\delta_{s}(x^{0};-T_{m})\right)^{\top}T_{m}^{\dagger}\end{bmatrix}.$

###### Proof.

We have

$\displaystyle\nabla_{c}^{2}f(x^{0};S,T_{1:m})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\nabla_{s}^{2}f(x^{0};S,T_{1:m})+\nabla_{s}^{2}f(x^{0};-S,-T_{1:m})\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}(S^{\top})^{\dagger}\left(\delta^{2}_{s}(x^{0};S,T_{1;m})-\delta^{2}_{s}(x^{0};-S,-T_{1:m})\right)$
$\displaystyle=$
$\displaystyle(S^{\top})^{\dagger}\frac{1}{2}\begin{bmatrix}\left(\nabla_{s}f(x^{0}+s^{1};T_{1})-\nabla_{s}f(x^{0};T_{1})\right)^{\top}-\left(\nabla_{s}f(x^{0}-s^{1};-T_{1})-\nabla_{s}f(x^{0};-T_{1})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}+s^{m};T_{m})-\nabla_{s}f(x^{0};T_{m})\right)^{\top}-\left(\nabla_{s}f(x^{0}-s^{m};-T_{m})-\nabla_{s}f(x^{0};-T_{m})\right)^{\top}\end{bmatrix}$
$\displaystyle=$
$\displaystyle(S^{\top})^{\dagger}\frac{1}{2}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};T_{1})-\delta_{s}(x^{0};T_{1})\right)^{\top}T_{1}^{\dagger}-\left(\left(\delta_{s}(x^{0}-s^{1};-T_{1})-\delta_{s}(x^{0};-T_{1})\right)^{\top}(-T_{1})^{\dagger}\right)\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};T_{m})-\delta_{s}(x^{0};T_{m})\right)^{\top}T_{m}^{\dagger}-\left(\left(\delta_{s}(x^{0}-s^{m};-T_{m})-\delta_{s}(x^{0};-T_{m})\right)^{\top}(-T_{m})^{\dagger}\right)\end{bmatrix}$
$\displaystyle=$
$\displaystyle(S^{\top})^{\dagger}\frac{1}{2}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};T_{1})+\delta_{s}(x^{0}-s^{1};-T_{1})-\delta_{s}(x^{0};T_{1})-\delta_{s}(x^{0};-T_{1})\right)^{\top}T_{1}^{\dagger}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};T_{m})+\delta_{s}(x^{0}-s^{m};-T_{m})-\delta_{s}(x^{0};T_{m})-\delta_{s}(x^{0};-T_{m})\right)^{\top}T_{m}^{\dagger}\end{bmatrix}$
$\displaystyle=$
$\displaystyle(S^{\top})^{\dagger}\delta_{c}^{2}(x^{0};S,T_{1:m}).$

∎

Notice that

$\displaystyle\left[\delta_{s}(x^{0}+s^{i};T_{i})+\delta_{s}(x^{0}-s^{i};-T_{i})-\delta_{s}(x^{0};T_{i})-\delta_{s}(x^{0};-T_{i})\right]_{j_{i}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}f(x^{0}+s^{i}+t_{i}^{j_{i}})-f(x^{0}+s^{i})+f(x^{0}-s^{i}-t_{i}^{j_{i}})-f(x^{0}-s^{i})$
$\displaystyle-f(x^{0}+t_{i}^{j_{i}})+f(x^{0})-f(x^{0}-t_{i}^{j_{i}})+f(x^{0})\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}f(x^{0}+s^{i}+t_{i}^{j_{i}})+f(x^{0}-s^{i}-t_{i}^{j_{i}})$
$\displaystyle-f(x^{0}+s^{i})-f(x^{0}-s^{i})-f(x^{0}+t_{i}^{j_{i}})-f(x^{0}-t_{i}^{j_{i}})+2f(x^{0})\big{)}$
(6)

for $j_{i}\in\\{1,\dots,k_{i}\\}$ and $i\in\\{1,\dots,m\\}.$ We will use (6)
in Theorem 4.3.

Next, we show that the GCSH is actually a special case of the GSH where the
matrices of directions have a specific form. This result is similar to the
result in [28, Proposition 2.10], which showed that the generalized centered
gradient is equal to the generalized simplex gradient when the matrix of
directions used to compute the generalized simplex gradient has the form
$A=\begin{bmatrix}S&-S\end{bmatrix}$ for some matrix
$S\in\operatorname{\mathbb{R}}^{n\times m}.$ We shall use the following lemma.

######  Lemma 3.2.

Let $A=\begin{bmatrix}S&-S\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
2m}$ for some matrix $S\in\operatorname{\mathbb{R}}^{n\times m}.$ Then

$A^{\dagger}=\frac{1}{2}\begin{bmatrix}S^{\dagger}\\\
-S^{\dagger}\end{bmatrix}.$

###### Proof.

This follows from basic algebra. ∎

######  Proposition 3.3.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$A=\begin{bmatrix}a^{1}&a^{2}&\cdots&a^{2m}\end{bmatrix}=\begin{bmatrix}S&-S\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
2m}$ for some
$S=\begin{bmatrix}s^{1}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and let $B_{i}=T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}},$ and
$B_{m+i}=-T_{i}$ for all $i\in\\{1,\dots,m\\}.$ Suppose that $x^{0}\oplus
S\oplus T_{i},x^{0}\oplus(-S)\oplus(-T_{i}),x^{0}\oplus(\pm S),x^{0}\oplus(\pm
T_{i})$ are contained in $\operatorname{dom}f$ for all $i\in\\{1,\dots,m\\}.$
Then

$\nabla_{s}^{2}f(x^{0};A,B_{1:2m})=\nabla_{c}^{2}f(x^{0};S,T_{1:m}).$

###### Proof.

We have

$\displaystyle\nabla_{s}^{2}f(x^{0};A,B_{1:2m})$
$\displaystyle=\begin{bmatrix}S^{\top}\\\
-S^{\top}\end{bmatrix}^{\dagger}\begin{bmatrix}\left(\nabla_{s}f(x^{0}+a^{1};B_{1})-\nabla_{s}f(x^{0};B_{1})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}+a^{2m};B_{2m})-\nabla_{s}f(x^{0};B_{2m})\right)^{\top}\end{bmatrix}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}(S^{\top})^{\dagger}&-(S^{\top})^{\dagger}\end{bmatrix}\begin{bmatrix}\left(\nabla_{s}f(x^{0}+a^{1};B_{1})-\nabla_{s}f(x^{0};B_{1})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}+a^{2m};B_{2m})-\nabla_{s}f(x^{0};B_{2m})\right)^{\top}\end{bmatrix}$
$\displaystyle=\frac{1}{2}(S^{\top})^{\dagger}\begin{bmatrix}\left(\nabla_{s}f(x^{0}+s^{1};T_{1})-\nabla_{s}f(x^{0};T_{1})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}+s^{m};T_{m})-\nabla_{s}f(x^{0};T_{m})\right)^{\top}\end{bmatrix}+\frac{1}{2}(-S^{\top})^{\dagger}\begin{bmatrix}\left(\nabla_{s}f(x^{0}-s^{1};-T_{1})-\nabla_{s}f(x^{0};-T_{1})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}-s^{m};-T_{m})-\nabla_{s}f(x^{0};-T_{m})\right)^{\top}\end{bmatrix}$
$\displaystyle=\frac{1}{2}\left(\nabla_{s}^{2}f(x^{0};S,T_{1:m})+\nabla_{s}^{2}f(x^{0};-S,-T_{1:m})\right)$
$\displaystyle=\nabla_{c}^{2}f(x^{0};S,T_{1:m}).$

∎

We are now ready to examine the error bounds for the GSH and the GCSH.

##  4 Error bounds

In this section, we provide error bounds for the GSH and the GCSH. These error
bounds are separated in two categories: the general case where all matrices
$T_{i}$ are not necessarily equal, and the special case where all matrices
$T_{i}=:\overline{T}$ are equal. The error bounds for the GSH show an order-1
accurate approximation of the full Hessian $\nabla^{2}f(x^{0})$, or an order-1
accurate approximation of the partial Hessian
$\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$ if the GSH is $S$\- or
$T_{1:m}$-underdetermined. Similarly, the error bounds for the GCSH show an
order-2 accurate approximation of the full Hessian $\nabla^{2}f(x^{0})$, or an
order-2 accurate approximation of the partial Hessian
$\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$ if the GCSH is $S$\- or
$T_{1:m}$-underdetermined.

Before introducing the error bounds, we note that in certain situations, the
projection of the GSH (GCSH) onto $S$ and $T_{1:m}$ is equal to the GSH
(GCSH).

######  Proposition 4.1.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and $x^{0}\in\operatorname{dom}f.$ Let
$S=\begin{bmatrix}s^{1}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and $T_{i}\in\operatorname{\mathbb{R}}^{n\times k_{i}}.$ Assume that
$x^{0}\oplus(\pm S),x^{0}\oplus(\pm T_{i}),x^{0}\oplus S\oplus T_{i}$ and
$x^{0}\oplus-S\oplus-T_{i}$ are contained in $\operatorname{dom}f$ for all
$i.$ Then the following hold.

  1. (i)

If $S$ is full column rank or $T_{i}$ is full row rank for all
$i\in\\{1,\dots,m\\}$, then

$\displaystyle\operatorname{Proj}_{S,T_{1:m}}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=\nabla_{s}^{2}f(x^{0};S,T_{1:m})$ (7)

and

$\displaystyle\operatorname{Proj}_{S,T_{1:m}}\nabla_{c}^{2}f(x^{0};S,T_{1:m})=\nabla_{c}^{2}f(x^{0};S,T_{1:m}).$
(8)

  2. (ii)

If $T_{1}=T_{2}=\dots=T_{m}=:\overline{T},$ then

$\displaystyle\operatorname{Proj}_{S,\overline{T}}\nabla_{s}^{2}f(x^{0};S,\overline{T})$
$\displaystyle=\nabla_{s}^{2}f(x^{0};S,\overline{T})$ (9)

and

$\displaystyle\operatorname{Proj}_{S,\overline{T}}\nabla_{c}^{2}f(x^{0};S,\overline{T})$
$\displaystyle=\nabla_{c}^{2}f(x^{0};S,\overline{T}).$ (10)

###### Proof.

Suppose $S$ is full column rank. Since
$T_{i}^{\dagger}T_{i}T_{i}^{\dagger}=T_{i}^{\dagger}$ by Property (ii) of the
Moore–Penrose pseudoinverse, we obtain

$\displaystyle\operatorname{Proj}_{S,T_{1:m}}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}\nabla_{s}^{2}f(x^{0};S,T_{1:m})T_{i}T_{i}^{\dagger}$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m})T_{i}T_{i}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}\sum_{i=1}^{m}\operatorname{Diag}(e^{i}_{m})\delta^{2}_{s}(x^{0};S,T_{1:m})T_{i}T_{i}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m})$
$\displaystyle=\nabla_{s}^{2}f(x^{0};S,T_{1:m}).$

Equation (8) follows analogously.

Now suppose that $T_{i}$ is full row rank for all $i\in\\{1,\dots,m\\}.$ Then
we know $T_{i}T_{i}^{\dagger}=\operatorname{Id}_{n}$ for all
$i\in\\{1,\dots,m\\}.$ We obtain

$\displaystyle\operatorname{Proj}_{S,T_{1:m}}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}\nabla_{s}^{2}f(x^{0};S,T_{1:m})T_{i}T_{i}^{\dagger}$
$\displaystyle=\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=(S^{\top})^{\dagger}\left(\sum_{i=1}^{m}e^{i}_{m}(e^{i}_{m})^{\top}\right)S^{\top}\nabla_{s}^{2}f(x^{0};S,T_{1:m})$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m})$
$\displaystyle=(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m})=\nabla_{s}^{2}f(x^{0};S,T_{1:m}).$

Equation (8) follows analogously.

Finally, suppose that $T_{1}=T_{2}=\dots=T_{m}=:\overline{T}$. We have

$\displaystyle\operatorname{Proj}_{S,\overline{T}}\nabla_{s}^{2}f(x^{0};S,\overline{T})$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}\nabla_{s}^{2}f(x^{0};S,\overline{T})\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,\overline{T})\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}(S^{\top})^{\dagger}\begin{bmatrix}\left(\nabla_{s}f(x^{0}+s^{1};\overline{T})-\nabla_{s}f(x^{0};\overline{T})\right)^{\top}\\\
\vdots\\\
\left(\nabla_{s}f(x^{0}+s^{m};\overline{T})-\nabla_{s}f(x^{0};\overline{T})\right)^{\top}\end{bmatrix}\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}(S^{\top})^{\dagger}\begin{bmatrix}\left((\overline{T}^{\top})^{\dagger}\delta_{s}(x^{0}+s^{1};\overline{T})-(\overline{T}^{\top})^{\dagger}\delta_{s}(x^{0};\overline{T})\right)^{\top}\\\
\vdots\\\
\left((\overline{T}^{\top})^{\dagger}\delta_{s}(x^{0}+s^{m};\overline{T})-(\overline{T}^{\top})^{\dagger}\delta_{s}(x^{0};\overline{T})\right)^{\top}\end{bmatrix}\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}S^{\top}(S^{\top})^{\dagger}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};\overline{T})-\delta_{s}(x^{0};\overline{T})\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};\overline{T})-\delta_{s}(x^{0};\overline{T})\right)^{\top}\end{bmatrix}\overline{T}^{\dagger}\overline{T}\overline{T}^{\dagger}$
$\displaystyle=(S^{\top})^{\dagger}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};\overline{T})-\delta_{s}(x^{0};\overline{T})\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};\overline{T})-\delta_{s}(x^{0};\overline{T})\right)^{\top}\end{bmatrix}\overline{T}^{\dagger}$

by Property (ii) of the Moore-Penrose pseudoinverse. Thus,
$\operatorname{Proj}_{S,\overline{T}}\nabla_{s}^{2}f(x^{0};S,\overline{T})=\nabla_{s}^{2}f(x^{0};S,\overline{T}).$
Equation (10) follows analogously.

∎

In the following theorem and the remainder of this paper, we use the notation

$\displaystyle\widehat{T}_{i}$ $\displaystyle=T_{i}/\Delta_{T_{i}},\quad
i\in\\{1,\dots,m\\},$ $\displaystyle\Delta_{u}$
$\displaystyle=\max\\{\Delta_{S},\Delta_{T_{1}},\dots,\Delta_{T_{m}}\\},$
$\displaystyle\Delta_{l}$
$\displaystyle=\min\\{\Delta_{S},\Delta_{T_{1}},\dots,\Delta_{T_{m}}\\},$
$\displaystyle\Delta_{T}$
$\displaystyle=\max\\{\Delta_{T_{1}},\dots,\Delta_{T_{m}}\\}$
$\displaystyle\widehat{T}$ $\displaystyle=\widehat{T}_{i}\quad\text{such
that}\quad\left\|\widehat{T}_{i}^{\dagger}\right\|\quad\text{is maximal,}\quad
i\in\\{1,\dots,m\\},$ $\displaystyle k$
$\displaystyle=\max\\{k_{1},\dots,k_{m}\\},$ $\displaystyle H$
$\displaystyle=\nabla^{2}f(x^{0}).$

######  Theorem 4.2 (Error bounds for the GSH).

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
be $\mathcal{C}^{3}$ on $B(x^{0},\overline{\Delta})$ where
$x^{0}\in\operatorname{dom}f$ is the point of interest and
$\overline{\Delta}>0$. Denote by $L_{\nabla^{2}f}\geq 0$ the Lipschitz
constant of $\nabla^{2}f$ on $\overline{B}(x^{0},\overline{\Delta})$
respectively. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and
$T_{i}=\begin{bmatrix}t^{1}_{i}&t_{i}^{2}&\cdots&t_{i}^{k_{i}}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
k_{i}}$ for all $i\in\\{1,\dots,m\\}.$ Assume that
$B(x^{0},\Delta_{T_{i}})\subset B(x^{0},\overline{\Delta})$ and
$B(x^{0}+s^{i},\Delta_{T_{i}})\subset B(x^{0},\overline{\Delta})$ for all
$i\in\\{1,\dots m\\}.$

  1. (i)

If $S$ is full column rank or $T_{i}$ is full row rank for all
$i\in\\{1,\dots,m\\}$, then

$\displaystyle\|\operatorname{Proj}_{S,T_{1:m}}\nabla_{s}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})\|$
$\displaystyle=\|\nabla_{s}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})\|$
$\displaystyle\leq
4m\sqrt{k}L_{\nabla^{2}f}\|(\widehat{S}^{\top})^{\dagger}\|\|\widehat{T}^{\dagger}\|\left(\frac{\Delta_{u}}{\Delta_{l}}\right)^{2}\Delta_{u}.$
(11)

  2. (ii)

If $T_{1}=T_{2}=\dots=T_{m}=:\overline{T}$, then

$\displaystyle\|\operatorname{Proj}_{S,\overline{T}}\nabla_{s}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}\nabla^{2}f(x^{0})\|$
$\displaystyle=\|\nabla_{s}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}\nabla^{2}f(x^{0})\|$
$\displaystyle\leq
4\sqrt{mk}L_{\nabla^{2}f}\frac{\Delta_{u}}{\Delta_{l}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\Delta_{u}.$
(12)

###### Proof.

The equalities in Items $(i)$ and $(ii)$ follow from (7) and (9) respectively.
We have

$\displaystyle\|\nabla_{s}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}H\|$
$\displaystyle=\left\|(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,T_{1:m})-\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}HT_{i}T_{i}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\delta^{2}_{s}(x^{0};S,T_{1:m})-\sum_{i=1}^{m}(s^{i})^{\top}HT_{i}T_{i}^{\dagger}\right\|$
$\displaystyle=\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\sum_{i=1}^{m}(e^{i}_{m})^{\top}\delta^{2}_{s}(x^{0};S,T_{1:m})-(s^{i})^{\top}HT_{i}T_{i}^{\dagger}\right\|$
$\displaystyle\leq\sum_{i=1}^{m}\frac{1}{\Delta_{S}\Delta_{T_{i}}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T_{i}}^{\dagger}\right\|\left\|\left(\delta_{s}(x^{0}+s^{i};T_{i})-\delta_{s}(x^{0};T_{i})\right)^{\top}-(s^{i})^{\top}HT_{i}\right\|.$
(13)

Let
$u_{i}=\left(\delta_{s}(x^{0}+s^{i};T_{i})-\delta_{s}(x^{0};T_{i})\right)^{\top}\in\operatorname{\mathbb{R}}^{k_{i}}$
for all $\\{1,\dots,m\\}.$ Next we find a bound for each
$\left[u_{i}-(s^{i})^{\top}HT_{i}\right]_{j_{i}},j_{i}\in\\{1,\dots,k_{i}\\}.$
Note that

$\displaystyle\left[u_{i}-(s^{i})^{\top}HT_{i}\right]_{j_{i}}$
$\displaystyle=f(x^{0}+s^{i}+t_{i}^{j_{i}})-f(x^{0}+s^{i})-f(x^{0}+t_{i}^{j_{i}})+f(x^{0})-(s^{i})^{\top}Ht_{i}^{j_{i}}.$

Each function value $f(x^{0}+s^{i}+t_{i}^{j_{i}}),f(x^{0}+s^{i}),$ and
$f(x^{0}+t_{i}^{j_{i}}),$ can be written as a second-order Taylor expansion
about $x^{0}$ plus a remainder term $R_{2}(\cdot).$ Hence, we obtain

$\displaystyle\left|f(x^{0}+s^{i}+t_{i}^{j_{i}})-f(x^{0}+s^{i})-f(x^{0}+t_{i}^{j_{i}}+f(x^{0})-(s^{i})^{\top}HT_{i}^{j_{i}}\right|$
$\displaystyle\leq|R_{2}(x^{0}+s^{i}+t_{i}^{j_{i}})|+|R_{2}(x^{0}+s^{i})|+|R_{2}(x^{0}+t_{i}^{j_{i}})|$
$\displaystyle\leq\
\frac{1}{6}L_{\nabla^{2}f}\|s^{i}+t_{i}^{j_{i}}\|^{3}+\frac{1}{6}L_{\nabla^{2}f}\|s^{i}\|^{3}+\frac{1}{6}L_{\nabla^{2}f}\|t_{i}^{j_{i}}\|^{3}$
$\displaystyle\leq\frac{1}{2}L_{\nabla^{2}f}(\Delta_{S}+\Delta_{T_{i}})^{3}.$
(14)

Plugging (14) in (13), we obtain

$\displaystyle\|\nabla_{s}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}H\|$
$\displaystyle\leq\sum_{i=1}^{m}\frac{1}{\Delta_{S}\Delta_{T_{i}}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T_{i}}^{\dagger}\right\|\frac{\sqrt{k_{i}}}{2}L_{\nabla^{2}f}(\Delta_{S}+\Delta_{T_{i}})^{3}$
$\displaystyle\leq\sum_{i=1}^{m}\frac{1}{\Delta_{l}^{2}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\frac{\sqrt{k}}{2}L_{\nabla^{2}f}(2\Delta_{u})^{3}$
$\displaystyle=4m\sqrt{k}L_{\nabla^{2}f}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\left(\frac{\Delta_{u}}{\Delta_{l}}\right)^{2}\Delta_{u}.$

Now, suppose that $T_{1}=T_{2}=\dots=T_{m}=:\overline{T}.$ We have

$\displaystyle\|\nabla_{s}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}H\|$
$\displaystyle=\left\|(S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,\overline{T})-(S^{\top})^{\dagger}S^{\top}H\overline{T}\overline{T}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\delta^{2}_{s}(x^{0};S,\overline{T})-S^{\top}H\overline{T}\overline{T}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}\Delta_{T}}\|(\widehat{S}^{\top})^{\dagger}\|\widehat{T}^{\dagger}\|\|\|(D_{s}^{2})^{\top}-S^{\top}H\overline{T}\|,$
(15)

where

$(D_{s}^{2})^{\top}=\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};\overline{T}))-\delta_{s}(x^{0};\overline{T})\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};\overline{T})-\delta_{s}(x^{0};\overline{T})\right)^{\top}\end{bmatrix}\in\operatorname{\mathbb{R}}^{m\times
k}.$

Each entry $\left[(D_{s}^{2})^{\top}-S^{\top}H\overline{T}\right]_{i,j}$ has
the previous bound obtained in (14). Since

$\|(D_{s}^{2})^{\top}-S^{\top}H\overline{T}\|\leq\|(D_{s}^{2})^{\top}-S^{\top}H\overline{T}\|_{F},$

using (14) in (15), we get the inequality

$\displaystyle\|\nabla_{s}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}H\|$
$\displaystyle\leq\frac{\sqrt{mk}}{2}L_{\nabla^{2}f}\frac{(\Delta_{S}+\Delta_{T})^{3}}{\Delta_{S}\Delta_{T}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|$
$\displaystyle\leq
4\sqrt{mk}L_{\nabla^{2}f}\frac{\Delta_{u}}{\Delta_{l}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\Delta_{u}.$

∎

The previous error bounds show that the GSH is an order-1 accurate
approximation of a partial Hessian, or an order-1 accurate approximation of
the full Hessian when the GSH is $S$-determined\overdetermined and
$T_{1:m}$-determined\overdetermined (equivalently,
$(S^{\top})^{\dagger}S^{\top}=\operatorname{Id}_{n}$ and
$T_{i}T_{i}^{\dagger}=\operatorname{Id}_{n}$ for all $i$). The ratio
$\frac{\Delta_{u}}{\Delta_{l}}$ suggests that the radii of the matrices should
be taken to be of the same magnitude. If one of the radius is decreased, it
suggests that all radii should be decreased by the same amount. Finally, the
presence of $m$ and $k$ in the error bounds roughly suggests that there is no
advantage of having $S$-overdetermined or $T_{1:m}$-overdetermined. More
details on this topic can be found in [11, 27].

When all matrices $T_{i}$ are equal, note that Item $(ii)$ covers all 16 cases
from Definition 2.9. When all matrices $T_{i}$ are not equal, Item $(i)$ in
the previous theorem covers 12 out of the 16 cases. It does not cover the case
$S$-nondetermined\overdetermined and $T_{1:m}$-nondetermined \underdetermined.
These cases are left as a future research direction and discuss in the
conclusion.

######  Theorem 4.3 (Error bounds for the GCSH).

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
be $\mathcal{C}^{4}$ on $B(x^{0},\overline{\Delta})$ where
$x^{0}\in\operatorname{dom}f$ is the point of interest and
$\overline{\Delta}>0$. Denote by $L_{\nabla^{3}f}$ the Lipschitz constant of
$\nabla^{3}f$ on $\overline{B}(x^{0},\overline{\Delta}).$ Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m},$
$T_{i}=\begin{bmatrix}t_{i}^{1}&t_{i}^{2}&\cdots&t_{i}^{k_{i}}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
k_{i}}$ with $\overline{B}(x^{0}\pm s^{i},\Delta_{T_{i}})\subset
B(x^{0},\overline{\Delta})$ for all $i\in\\{1,\dots,m\\}$.

  1. (i)

If $S$ is full column rank or $T_{i}$ is full row rank for all
$i\in\\{1,\dots,m\\}$, then

$\displaystyle\|\operatorname{Proj}_{S,T_{1:m}}\nabla_{c}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})\|$
$\displaystyle=\|\nabla_{c}^{2}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})\|$
$\displaystyle\leq$ $\displaystyle
2m\sqrt{k}L_{\nabla^{3}f}\left(\frac{\Delta_{u}}{\Delta_{l}}\right)^{2}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\left(\widehat{T}\right)^{\dagger}\right\|\Delta^{2}_{u}.$
(16)

  2. (ii)

If $T_{1}=T_{2}=\dots=T_{m}=:\overline{T}\in\operatorname{\mathbb{R}}^{n\times
k},$ then

$\displaystyle\left\|\operatorname{Proj}_{S,\overline{T}}\nabla_{c}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}\nabla^{2}f(x^{0})\right\|$
$\displaystyle=\left\|\nabla_{c}^{2}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}\nabla^{2}f(x^{0})\right\|$
$\displaystyle\leq
2\sqrt{mk}L_{\nabla^{3}f}\frac{\Delta_{u}}{\Delta_{l}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\left(\widehat{T}\right)^{\dagger}\right\|\Delta^{2}_{u}.$
(17)

###### Proof.

The equalities in Items $(i)$ and $(ii)$ follow from (8) and (10)
respectively. We have

$\displaystyle\left\|\nabla^{2}_{c}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}H\right\|$
$\displaystyle=\left\|(S^{\top})^{\dagger}\delta_{c}^{2}(x^{0};S,T_{1:m})-\sum_{i=1}^{m}(S^{\top})^{\dagger}e^{i}_{m}(e^{i}_{m})^{\top}S^{\top}HT_{i}T_{i}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\delta_{c}^{2}(x^{0};S,T_{1:m})-\sum_{i=1}^{m}(s^{i})^{\top}HT_{i}T_{i}^{\dagger}\right\|$
$\displaystyle=\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\sum_{i=1}^{m}(e^{i}_{m})^{\top}\delta_{c}^{2}(x^{0};S,T_{1:m})-(s^{i})^{\top}HT_{i}T_{i}^{\dagger}\right\|$

Let
$w_{i}=\frac{1}{2}\left(\delta_{s}(x^{0}+s^{i};T_{i})+\delta_{s}(x^{0}-s^{i};-T_{i})-\delta_{s}(x^{0};T_{i})-\delta_{s}(x^{0};-T_{i})\right)^{\top}\in\operatorname{\mathbb{R}}^{k_{i}}.$
Next we find a bound for each entry
$\left[w_{i}-(s^{i})^{\top}HT_{i}\right]_{j_{i}},j_{i}\in\\{1,\dots,k_{i}\\}.$
Using (6), we know that

$\displaystyle\left[w_{i}-(s^{i})^{\top}HT_{i}\right]_{j_{i}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(f(x^{0}+s^{i}+t_{i}^{j_{i}})+f(x^{0}-s^{i}-t_{i}^{j_{i}})-f(x^{0}+s^{i})-f(x^{0}-s^{i})-f(x^{0}+t_{i}^{j_{i}})-f(x^{0}-t_{i}^{j_{i}})+2f(x^{0})\right)$
$\displaystyle-(s^{i})^{\top}Ht_{i}^{j_{i}}.$ (18)

Each function value
$f(x^{0}+s^{i}+t_{i}^{j_{i}}),f(x^{0}-s^{i}-t_{i}^{j_{i}}),f(x^{0}+s^{i}),f(x^{0}-s^{i}),f(x^{0}+t_{i}^{j_{i}})$
and $f(x^{0}-t_{i}^{j_{i}})$ can be written as a third-order Taylor expansion
about $x^{0}$ plus a remainder term $R_{3}(\cdot).$ It follows that

$\displaystyle\medmath{\left|\frac{1}{2}\left(f(x^{0}+s^{i}+t_{i}^{j_{i}})+f(x^{0}-s^{i}-t_{i}^{j_{i}})-f(x^{0}+s^{i})-f(x^{0}-s^{i})-f(x^{0}+t_{i}^{j_{i}})-f(x^{0}-t_{i}^{j_{i}})+2f(x^{0})\right)-(s^{i})^{\top}Ht_{i}^{j_{i}}\right|}$
$\displaystyle\leq\medmath{\frac{1}{2}\left(\|R_{3}(x^{0}+s^{i}+t_{i}^{j_{i}})\|+\|R_{3}(x^{0}-s^{i}-t_{i}^{j_{i}})\|+\|R_{3}(x^{0}+s^{i})\|+\|R_{3}(x^{0}-s^{i})\|+\|R_{3}(x^{0}+t_{i}^{j_{i}}\|+\|R_{3}(x^{0}-t_{i}^{j_{i}})\|\right)}$
$\displaystyle\leq\medmath{\frac{1}{48}\left(L_{\nabla^{3}f}\|s^{i}+t_{i}^{j_{i}}\|^{4}+\frac{1}{24}L_{\nabla^{3}f}\|-s^{i}-t_{i}^{j_{i}}\|^{4}+L_{\nabla^{3}f}\|-s^{i}\|^{4}+L_{\nabla^{3}f}\|s^{i}\|^{4}+L_{\nabla^{3}f}\|t_{i}^{j_{i}}\|+L_{\nabla^{3}f}\|-t_{i}^{j_{i}}\|\right)}$
$\displaystyle\leq\frac{1}{8}L_{\nabla^{3}f}(\Delta_{S}+\Delta_{T_{i}})^{4}.$
(19)

Plugging the bound from (19) in (18), we get

$\displaystyle\|\nabla^{2}_{c}f(x^{0};S,T_{1:m})-\operatorname{Proj}_{S,T_{1:m}}H\|$
$\displaystyle\leq\sum_{i=1}^{m}\frac{1}{\Delta_{S}\Delta_{T_{i}}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T_{i}}^{\dagger}\right\|\sqrt{k_{i}}\frac{1}{8}L_{\nabla^{3}f}(\Delta_{S}+\Delta_{T_{i}})^{4}$
$\displaystyle\leq\sum_{i=1}^{m}\frac{1}{\Delta_{l}^{2}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\sqrt{k}\frac{1}{8}L_{\nabla^{3}f}(2\Delta_{u})^{4}$
$\displaystyle=2m\sqrt{k}L_{\nabla^{3}f}\left(\frac{\Delta_{u}}{\Delta_{l}}\right)^{2}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\left(\widehat{T}\right)^{\dagger}\right\|\Delta^{2}_{u}.$

Now, suppose that $T_{1}=T_{2}=\dots=T_{m}=:\overline{T}.$ We have

$\displaystyle\|\nabla^{2}_{c}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}H\|$
$\displaystyle=\left\|(S^{\top})^{\dagger}\delta_{c}^{2}(x^{0};S,\overline{T})-(S^{\top})^{\dagger}S^{\top}H\overline{T}\overline{T}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\delta_{c}^{2}(x^{0};S,\overline{T})-S^{\top}H\overline{T}\overline{T}^{\dagger}\right\|$
$\displaystyle\leq\frac{1}{\Delta_{S}\Delta_{T}}\|(\widehat{S}^{\top})^{\dagger}\|\widehat{T}^{\dagger}\|\|\|(D_{c}^{2})^{\top}-S^{\top}H\overline{T}\|,$
(20)

where

$(D_{c}^{2})^{\top}=\frac{1}{2}\begin{bmatrix}\left(\delta_{s}(x^{0}+s^{1};\overline{T})+\delta_{s}(x^{0}-s^{1};-\overline{T})-\delta_{s}(x^{0};\overline{T})-\delta_{s}(x^{0};-\overline{T})\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+s^{m};\overline{T})+\delta_{s}(x^{0}-s^{m};-\overline{T})-\delta_{s}(x^{0};\overline{T})-\delta_{s}(x^{0};-\overline{T})\right)^{\top}\end{bmatrix}\in\operatorname{\mathbb{R}}^{m\times
k}.$

Each entry $\left[(D_{c}^{2})^{\top}-S^{\top}H\overline{T}\right]_{i,j}$ is
bounded by (19). Since

$\|(D_{c}^{2})^{\top}-S^{\top}H\overline{T}\|\leq\|(D_{f}^{c})^{\top}-S^{\top}H\overline{T}\|_{F},$

using (19) in (20), we get the inequality

$\displaystyle\|\nabla^{2}_{c}f(x^{0};S,\overline{T})-\operatorname{Proj}_{S,\overline{T}}H\|$
$\displaystyle\leq\frac{\sqrt{mk}}{24}L_{\nabla^{3}f}\frac{(\Delta_{S}+\Delta_{T})^{4}}{\Delta_{S}\Delta_{T}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|$
$\displaystyle\leq\frac{2\sqrt{mk}}{3}L_{\nabla^{3}f}\frac{\Delta_{u}}{\Delta_{l}}\left\|(\widehat{S}^{\top})^{\dagger}\right\|\left\|\widehat{T}^{\dagger}\right\|\Delta_{u}^{2}.$

∎

The main differences between the error bounds presented for the GCSH in
Theorem 4.3 and the GSH in Theorem 4.2 are

  * •

the error bounds in Theorem 4.3 assume that $f\in\mathcal{C}^{4}$ and the
error bounds in Theorem 4.2 assume $f\in\mathcal{C}^{3},$

  * •

the error bounds in Theorem 4.3 involve the Lipschitz constant
$L_{\nabla^{3}f}$ and the error bounds in Theorem 4.2 involve the Lipschitz
constant $L_{\nabla^{2}f},$

  * •

the error bounds in Theorem 4.3 are order-2 accurate and the error bounds in
Theorem 4.2 are order-1 accurate.

Item $(ii)$ of Theorems 4.2 and 4.3 is more restrictive than Item $(i)$, but
in return it can offer additional benefits, such as symmetry. In the following
proposition, we show that the transpose of the GSH (GCSH) of $f$ at $x^{0}$
over $S$ and $\overline{T}$ is the same as the GSH (GCSH) of $f$ at $x^{0}$
over $\overline{T}$ and $S$.

######  Proposition 4.4.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
and let $x^{0}\in\operatorname{dom}f$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{m}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
m}$ and
$\overline{T}=\begin{bmatrix}t^{1}&t^{2}&\cdots&t^{k}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
k}$ with $x^{0}\oplus S\oplus
T_{i},x^{0}\oplus(-S)\oplus(-T_{i}),x^{0}\oplus(\pm S),$ and $x^{0}\oplus(\pm
T_{i})$ contained in $\operatorname{dom}f$ for all $i\in\\{1,\dots,m\\}.$Then

$\left(\nabla^{2}_{s}f(x^{0};S,\overline{T})\right)^{\top}=\nabla^{2}_{s}f(x^{0};\overline{T},S),$

and

$\left(\nabla^{2}_{c}f(x^{0};S,\overline{T})\right)^{\top}=\nabla^{2}_{c}f(x^{0};\overline{T},S).$

###### Proof.

We have

$\displaystyle\left(\nabla_{s}^{2}f(x^{0};S,\overline{T})\right)^{\top}$
$\displaystyle=\left((S^{\top})^{\dagger}\delta^{2}_{s}(x^{0};S,\overline{T})\right)^{\top}$
$\displaystyle=(\delta^{2}_{s}(x^{0};S,\overline{T}))^{\top}S^{\dagger}$
$\displaystyle=\begin{bmatrix}\nabla_{s}f(x^{0}+s^{1};\overline{T})-\nabla_{s}f(x^{0};\overline{T})&\cdots&\nabla_{s}f(x^{0}+s^{m};\overline{T})-\nabla_{s}f(x^{0};\overline{T})\end{bmatrix}S^{\dagger}$
$\displaystyle=(\overline{T}^{\top})^{\dagger}\begin{bmatrix}\delta_{s}(x^{0}+s^{1};\overline{T})-\delta_{s}(x^{0};\overline{T})&\cdots&\delta_{s}(x^{0}+s^{m};\overline{T})-\delta_{s}(x^{0};\overline{T})\end{bmatrix}S^{\dagger}.$

Let
$D_{s}^{2}=\begin{bmatrix}\delta_{s}(x^{0}+s^{1};\overline{T})-\delta_{s}(x^{0};\overline{T})&\cdots&\delta_{s}(x^{0}+s^{m};\overline{T})-\delta_{s}(x^{0};\overline{T})\end{bmatrix}\in\operatorname{\mathbb{R}}^{k\times
m}.$ Then

$\displaystyle[D_{s}^{2}]_{i,j}$
$\displaystyle=f(x^{0}+s^{i}+t^{j})-f(x^{0}+s^{i})-f(x^{0}+t^{j})+f(x^{0}).$

Note that

$\displaystyle D_{s}^{2}$
$\displaystyle=\begin{bmatrix}\left(\delta_{s}(x^{0}+t^{1};S)-\delta_{s}(x^{0};S)\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+t^{k};S)-\delta_{s}(x^{0};S)\right)^{\top}\end{bmatrix}.$

Hence,

$\displaystyle\left(\nabla_{s}^{2}f(x^{0};S,\overline{T})\right)^{\top}$
$\displaystyle=(\overline{T}^{\top})^{\dagger}\begin{bmatrix}\left(\delta_{s}(x^{0}+t^{1};S)-\delta_{s}(x^{0};S)\right)^{\top}\\\
\vdots\\\
\left(\delta_{s}(x^{0}+t^{k};S)-\delta_{s}(x^{0};S)\right)^{\top}\end{bmatrix}S^{\dagger}$
$\displaystyle=(\overline{T}^{\top})^{\dagger}\begin{bmatrix}\left((S^{\dagger})^{\top}\delta_{s}(x^{0}+t^{1};S)-(S^{\dagger})^{\top}\delta_{s}(x^{0};S)\right)^{\top}\\\
\vdots\\\
\left((S^{\dagger})^{\top}\delta_{s}(x^{0}+t^{k};S)-(S^{\dagger})^{\top}\delta_{s}(x^{0};S)\right)^{\top}\end{bmatrix}$
$\displaystyle=(\overline{T}^{\top})^{\dagger}\delta^{2}_{s}(x^{0};\overline{T},S)=\nabla_{s}^{2}f(x^{0};\overline{T},S).$

The second equality follows from the definition of the GCSH (Definition 2.7).
∎

In the next section, we investigate how to obtain an order-$1$ accurate
approximation of the full Hessian with a minimal number of sample points.

##  5 Minimal poised set and quadratic interpolation

In this section, we investigate how to choose the matrices of directions $S$
and $T_{i}$ to obtain an approximation of the full Hessian with a minimal
number of sample points. Reducing the number of distinct points is extremely
valuable, as it will decrease the number of distinct function evaluations
necessary to compute the GSH. For all of this section, we set
$T_{1}=\cdots=T_{m}=:\overline{T}\in\operatorname{\mathbb{R}}^{n\times k}$. We
show that if $S$ and $\overline{T}$ have a specific structure, then the number
of distinct function evaluations necessary to compute the GSH is
$(n+1)(n+2)/2.$ We then explore some results that occur when such a structure
is used. We begin by defining a _set for GSH computation_.

######  Definition 5.1.

The set of all distinct points utilized in the computation of
$\nabla_{s}^{2}f(x^{0};S,\overline{T})$ is said to be the set for GSH
computation and is denoted $\operatorname{\mathcal{S}}(x^{0};S,\overline{T}).$

Note that $\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ contains at most
$(m+1)(k+1)$ distinct points, but can contain fewer points if some of them
overlap. Next, we introduce the definition of _minimal poised set_ for the
GSH.

######  Definition 5.2 (Minimal poised set for the GSH).

Let $x^{0}\in\operatorname{\mathbb{R}}^{n}$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{n}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
n}$ and
$\overline{T}=\begin{bmatrix}t^{1}&t^{2}&\cdots&t^{n}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
n}$. We say that $\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ is a
_minimal poised set for the GSH_ at $x^{0}$ if and only if $S$ and
$\overline{T}$ are full rank and
$\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ contains exactly
$(n+1)(n+2)/2$ distinct points.

This definition requires $S$ and $\overline{T}$ to have exactly $n$ columns
and to be full rank. This is the case where the GSH is $S$-determined and
$\overline{T}$-determined. This implies that $S$ and $\overline{T}$ are of
minimal size to ensure that the GSH is an order-1 accurate approximation of
the full Hessian.

Note that a minimal poised set for the GCSH could be defined in a similar
manner. Obviously, such a set must contain more than $(n+1)(n+2)/2.$ An
investigation of this topic is left as a future research direction. We next
show that it is possible to create a minimal poised set for the GSH.

######  Proposition 5.3.

Let $x^{0}$ be the point of interest. Let
$S=\begin{bmatrix}s^{1}&s^{2}&\cdots&s^{n}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
n}$. Define the set $U_{\ell}$ for each index $\ell\in\\{0,1,\ldots,n\\}$ as

$U_{0}=S$

and

$U_{\ell}=\begin{bmatrix}s^{1}-s^{\ell}&s^{2}-s^{\ell}&\cdots&s^{\ell-1}-s^{\ell}&-s^{\ell}&s^{\ell+1}-s^{\ell}&\cdots&s^{n}-s^{\ell}\end{bmatrix}\in\operatorname{\mathbb{R}}^{n\times
n},\ell\neq 0.$

Then for each $\ell$,
$|\operatorname{\mathcal{S}}(x^{0};S,U_{\ell})|\leq(n+1)(n+2)/2$. Moreover, if
$S$ has full rank, then
$|\operatorname{\mathcal{S}}(x^{0};S,U_{\ell})|=(n+1)(n+2)/2$.

###### Proof.

Without loss of generality, let $x^{0}=\mathbf{0}$. First, suppose
$\ell\in\\{1,2,\dots,n\\}.$ For arbitrary function $f$, consider the matrix
$\delta^{2}_{s}(x^{0};S,\overline{T})$ where $\overline{T}=U_{\ell}.$ The
computation of $\nabla_{s}f(x^{0};U_{\ell})$ evaluates $f$ at the points

$\begin{array}[]{l}~{}~{}~{}\\{\mathbf{0},(s^{1}-s^{\ell}),\dots,(s^{\ell-1}-s^{\ell}),-s^{\ell},(s^{\ell+1}-s^{\ell}),\dots,(s^{n}-s^{\ell})\\}\\\
=\\{\mathbf{0}\\}\cup\\{-s^{\ell}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\left\\{s^{i}-s^{\ell}\right\\}.\end{array}$

For $i\neq\ell$, the computation of $\nabla_{s}f(x^{0}+s^{i};U_{k})$ evaluates
$f$ at the points

$\begin{array}[]{l}~{}~{}~{}\\{s^{i},s^{i}+(s^{1}-s^{\ell}),\dots,s^{i}+(s^{\ell-1}-s^{\ell}),s^{i}-s^{\ell},s^{i}+(s^{\ell+1}-s^{\ell}),\dots,s^{i}+(s^{n}-s^{\ell})\\}\\\
=\\{s^{i}\\}\cup\\{s^{i}-s^{\ell}\\}\cup\bigcup\limits_{\begin{subarray}{c}j=1\\\
j\neq\ell\end{subarray}}^{n}\\{s^{i}+s^{j}-s^{\ell}\\}.\end{array}$

The computation of $\nabla_{s}f(x^{0}+s^{\ell};U_{\ell})$ evaluates $f$ at the
points

$\begin{array}[]{l}~{}~{}~{}\\{s^{\ell},s^{\ell}+(s^{1}-s^{\ell}),\dots,s^{\ell}+(s^{\ell-1}-s^{\ell}),s^{\ell}-s^{\ell},s^{\ell}+(s^{\ell+1}-s^{\ell}),\dots,s^{\ell}+(s^{n}-s^{\ell})\\}\\\
=\\{s^{\ell}\\}\cup\\{\mathbf{0}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\\{s^{i}\\}=\\{\mathbf{0}\\}\cup S.\end{array}$

Thus, $f$ is evaluated at the points

$\begin{array}[]{l}\tiny{\bigg{(}\\{\mathbf{0}\\}\cup\\{-s^{\ell}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\left\\{s^{i}-s^{\ell}\right\\}\bigg{)}\cup\bigg{(}\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\\{s^{i}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\\{s^{i}-s^{\ell}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\bigcup\limits_{\begin{subarray}{c}j\geq i\\\
j\neq\ell\end{subarray}}\\{s^{i}+s^{j}-s^{\ell}\\}\bigg{)}\cup\bigg{(}\\{\mathbf{0}\\}\cup
S\bigg{)}}\\\ =\\{\mathbf{0}\\}\cup
S\cup\\{-s^{\ell}\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\left\\{s^{i}-s^{\ell}\right\\}\cup\bigcup\limits_{\begin{subarray}{c}i=1\\\
i\neq\ell\end{subarray}}^{n}\bigcup\limits_{\begin{subarray}{c}j\geq i\\\
j\neq\ell\end{subarray}}\\{s^{i}+s^{j}-s^{\ell}\\}.\end{array}$ (23)

This is at most $1+n+1+(n-1)+(n-1)(n-2)/2=(n+1)(n+2)/2$ points.

Now, suppose $\ell=0.$ Using a similar process to the above, we find that $f$
is evaluated at the points

$\begin{array}[]{l}\bigg{(}\\{\mathbf{0}\\}\cup\bigcup\limits_{i=1}^{n}\left\\{s^{i}\right\\}\bigg{)}~{}\cup~{}\bigg{(}\bigcup\limits_{i=1}^{n}\\{2s^{i}\\}\cup\bigcup\limits_{i=1}^{n}\bigcup\limits_{j>i}\\{s^{i}+s^{j}\\}\bigg{)}.\end{array}$
(25)

This is at most $(n+1)(n+2)/2$ points.

Finally, if $S$ is full rank, then the four sets in (23) and the four sets in
(25) are disjoint, so we have exactly $(n+1)(n+2)/2$ function evaluations.∎

Using $S=\operatorname{Id}_{n}$ in Proposition 5.3, we can create $n+1$
canonical minimal poised sets for the GSH.

######  Definition 5.4 ($\ell^{\mbox{th}}$-canonical minimal poised set for
the GSH).

Let $x^{0}\in\operatorname{\mathbb{R}}^{n}$ be the point of interest. Let
$S=\operatorname{Id}_{n}$. Fix $\ell\in\\{0,1,\dots,n\\}$. Let

$E_{0}=\operatorname{Id}_{n}$

and

$E_{\ell}=\begin{bmatrix}e^{1}-e^{\ell}&e^{2}-e^{\ell}&\cdots&e^{\ell-1}-e^{\ell}&-e^{\ell}&e^{\ell+1}-e^{\ell}&\cdots&e^{n}-e^{\ell}\end{bmatrix},\ell\neq
0.$

Then $\operatorname{\mathcal{S}}(x^{0};\operatorname{Id}_{n},E_{\ell})$ is
called the $\ell^{\mbox{th}}$-canonical minimal poised set for the GSH at
$x^{0}$.

From Proposition 5.3 and the fact that any matrix $E_{k}$ is full rank in
Definition 5.4, the $\ell^{\mbox{th}}$-canonical minimal poised set for the
GSH is indeed a minimal poised set for the GSH. Henceforth, we use the
notation $\operatorname{\mathcal{M}}(x^{0};S,U_{k})$ to denote a minimal
poised set for the GSH at $x^{0}$ that takes the form constructed in
Proposition 5.3.

Note that the order of the directions in $S$ and $\overline{T}$ is arbitrary.
Thus, it is immediately clear that if
$\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ is a minimal poised set for
the GSH and $P_{1},P_{2}\in\operatorname{\mathbb{R}}^{n\times n}$ are
permutation matrices, then
$\operatorname{\mathcal{S}}(x^{0};SP_{1},\overline{T}P_{2})$ is also a minimal
poised set for the GSH at $x^{0}$. The next proposition expands this idea and
demonstrates how to construct minimal poised sets for the GSH.

######  Proposition 5.5.

Let $x^{0}\in\operatorname{\mathbb{R}}^{n}$ be the point of interest. Let
$S,\overline{T}\in\operatorname{\mathbb{R}}^{n\times n}$. Let
$N\in\operatorname{\mathbb{R}}^{n\times n}$ be an invertible matrix and
$P_{1},P_{2}\in\operatorname{\mathbb{R}}^{n\times n}$ be permutation matrices.
Then $\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ is a minimal poised
set for the GSH at $x^{0}$ if and only if
$\operatorname{\mathcal{S}}(x^{0};NSP_{1},N\overline{T}P_{2})$ is a minimal
poised set for the GSH at $x^{0}$.

###### Proof.

The proof follows trivially from properties of matrices.∎

It follows that $\operatorname{\mathcal{S}}(x^{0};S,\overline{T})$ is a
minimal poised set for the GSH at $x^{0}$ if and only if the set
$\operatorname{\mathcal{S}}(x^{0};\beta SP_{1},\beta\overline{T}P_{2})$ is a
minimal poised set for the GSH at $x^{0}$, where $\beta$ is a nonzero scalar
and $P_{1},P_{2}$ are permutaion matrices.

######  Example 5.6.

Let $x^{0}=(0,0)$. The $2^{\mbox{nd}}$-canonical minimal poised set for the
GSH in $\operatorname{\mathbb{R}}^{2}$ contains the points
$(0,-1),(0,0),(0,1),(1,-1),(1,0)$ and $(2,-1)$. In this case,
$S=\\{e^{1},e^{2}\\}$ and $\overline{T}=\\{(e^{1}-e^{2}),-e^{2}\\}$. Figure 1
illustrates this set.

Figure 1: The $2^{\mbox{nd}}$-canonical minimal poised set for the GSH at
$x^{0}$ in $\operatorname{\mathbb{R}}^{2}$.

The points in $\\{(0,0)\\}\cup\\{(0,0)\oplus S\\}=\\{0,e^{1},e^{2}\\}$ are
represented with solid black borders. These are the base points where simplex
gradients will be computed. The lines represent the vectors corresponding to
$\overline{T}$ emanate from $\\{(0,0)\\}\cup\\{(0,0)\oplus S\\}$. The points
in $\left(\\{(0,0)\\}\cup\\{(0,0)\oplus S\\}\right)+\overline{T}$ are
represented with cyan cores. These are the points used to construct the
simplex gradients. Notice the points $(0,0)$ and $(1,0)$ have both black
borders and cyan cores. These are the common points that allow the number of
function evaluations to be reduced to $(n+1)(n+2)/2=6$.

We next demonstrate that every minimal poised set for the GSH of the form
$\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})$ is poised for quadratic
interpolation. We then show that the converse is not true; it is possible to
construct a set that is poised for quadratic interpolation, but does not take
the form of $\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})$.

######  Proposition 5.7.

Let $S\in\operatorname{\mathbb{R}}^{n\times n}$ be full rank. Select
$\ell\in\\{0,1,\ldots,n\\}$ and define $\overline{T}:=U_{\ell}$ as in
Proposition 5.3. Then $\operatorname{\mathcal{M}}(x^{0};S,\overline{T})$ is
poised for quadratic interpolation.

###### Proof.

Noting that $\operatorname{\mathcal{M}}(x^{0};S,\overline{T})$ is poised for
quadratic interpolation if and only if the set
$\operatorname{\mathcal{M}}(x^{0};S,\overline{T})\oplus-x^{0}$ is poised for
quadratic interpolation, we assume without loss of generality that
$x^{0}=\mathbf{0}$. Case I: $\overline{T}=U_{\ell}$ for
$\ell\in\\{1,2,\ldots,n\\}$. Without loss of generality, by Proposition 5.5,
assume that $\ell=n$. The points contained in
$\operatorname{\mathcal{M}}(x^{0};S,U_{n})$ are

$\displaystyle\\{\mathbf{0}\\}\cup
S\cup\\{-s^{n}\\}\cup\bigcup\limits_{i=1}^{n-1}\left\\{s^{i}-s^{n}\right\\}\cup\bigcup\limits_{i=1}^{n-1}\bigcup\limits_{\begin{subarray}{c}j\geq
i\\\ j\neq n\end{subarray}}\\{s^{i}+s^{j}-s^{n}\\}$ $\displaystyle=$
$\displaystyle\\{\mathbf{0}\\}\cup
S\cup\\{-s^{n}\\}\cup\bigcup\limits_{i=1}^{n-1}\left\\{s^{i}-s^{n}\right\\}\cup\\{2s^{i}-s^{n}\\}_{i=1}^{n-1}\cup\bigcup\limits_{i=1}^{n-1}\bigcup\limits_{\begin{subarray}{c}j>i\\\
j\neq n\end{subarray}}\left\\{s^{i}+{s^{j}}-s^{n}\right\\}.$

We show that using this set, the only solution to (1) is the trivial solution.
Considering the point $\mathbf{0}$, we obtain

$\displaystyle\alpha_{0}$ $\displaystyle=0.$

Considering the points $s^{i}$ for $i\in\\{1,2,\dots,n\\}$ and noting that for
all $i\in\\{1,2,\dots,n\\},$ $s^{i}=Se^{i}$, we obtain

$\displaystyle\overline{\alpha}^{\top}e^{i}+\frac{1}{2}(e^{i})^{\top}\widehat{\mathbf{H}}e^{i}$
$\displaystyle=0,$ (26)

where $\overline{\alpha}^{\top}=\alpha^{\top}S$ and
$\widehat{\mathbf{H}}=S^{\top}\mathbf{H}S.$ Note that $\widehat{\mathbf{H}}$
is symmetric. Considering the points $s^{i}-s^{n}$ for
$i\in\\{1,2,\dots,n-1\\},$ we obtain

$\overline{\alpha}^{\top}e^{i}-\overline{\alpha}^{\top}e^{n}-(e^{i})^{\top}\widehat{\mathbf{H}}e^{n}+\frac{1}{2}(e^{i})^{\top}\widehat{\mathbf{H}}e^{i}+\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}=0.$

Using (26), this simplifies to

$\displaystyle-\overline{\alpha}^{\top}e^{n}-(e^{i})^{\top}\widehat{\mathbf{H}}e^{n}+\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}=0.$
(27)

Considering the point $-s^{n},$ we find

$\displaystyle-\overline{\alpha}^{\top}e^{n}+\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}$
$\displaystyle=0,$ (28)

which reduces (27) to

$\displaystyle-(e^{i})^{\top}\widehat{\mathbf{H}}e^{n}=0\quad\mbox{for
}~{}i\in\\{1,2,\ldots,n-1\\}.$

Thus,
$(e^{i})^{\top}\widehat{\mathbf{H}}e^{n}=\widehat{\mathbf{H}}_{i,n}=\widehat{\mathbf{H}}_{n,i}=0$
for all $i\in\\{1,2,\dots,n-1\\}.$ Combining (26) at $i=n$ and (28) multiplied
by $-1$, we get

$\displaystyle\overline{\alpha}^{\top}e^{n}+\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}=0=\overline{\alpha}^{\top}e^{n}-\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n},$

which implies that
$(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}=\widehat{\mathbf{H}}_{n,n}=0$.
Considering the points $2s^{i}-s^{n}$ for $i\in\\{1,2,\dots,n-1\\},$ we get

$2\overline{\alpha}^{\top}e^{i}-\overline{\alpha}^{\top}e^{n}+2(e^{i})^{\top}\widehat{\mathbf{H}}e^{i}+\frac{1}{2}(e^{n})^{\top}\widehat{\mathbf{H}}e^{n}=0.$

Using (28), this simplifies to

$\displaystyle
2\overline{\alpha}^{\top}e^{i}+2(e^{i})^{\top}\widehat{\mathbf{H}}e^{i}=0.$

By multiplying (26) by 2 and substituting in the above equation, we get
$(e^{i})^{\top}\widehat{\mathbf{H}}e^{i}=\widehat{\mathbf{H}}_{i,i}=0$ for all
$i\in\\{1,2,\dots,n-1\\}.$ This now implies
$\overline{\alpha}_{i}=\overline{\alpha}^{\top}e^{i}=0$ for all
$i\in\\{1,2,\ldots,n\\}$, i.e., $\overline{\alpha}=\mathbf{0}$. Lastly,
consider the points $s^{i}+s^{j}-s^{n}$ for $i\neq
j,i,j\in\\{1,2,\dots,n-1\\}.$ Since $\widehat{\mathbf{H}}_{i,i}=0$ for
$i\in\\{1,2,\dots,n\\}$,
$\widehat{\mathbf{H}}_{i,n}=\widehat{\mathbf{H}}_{n,i}=0,$ for
$i\in\\{1,2,\dots,n-1\\}$ and $\overline{\alpha}=\mathbf{0}$, we obtain

$(e^{i})^{\top}\widehat{\mathbf{H}}e^{j}=0.$

Thus $\widehat{\mathbf{H}}=\mathbf{0}_{n\times n}$. Therefore, the only
solution to (1) is the trivial solution.

Case II: $\overline{T}=U_{0}.$ The proof for this case is analogous to that of
Case I. ∎

It follows from the previous proposition that when a minimal poised set for
the GSH is used, the GSH is equal to the simplex Hessian as described in [17,
Section 9.5] (see also [21, Section 3]). Therefore, the results obtained in
[21, 17] related to the simplex Hessian applies to the GSH. In this case, the
GSH (and the simplex Hessian) is equal to the Hessian of the quadratic
interpolation model passing through the sample points
$\operatorname{\mathcal{M}}(x^{0};S,U_{\ell}).$ Therefore, the results
developed in [17, 16, 21] related to the simplex Hessian and the Hessian of a
quadratic interpolation model are valid for the GSH. One of the advantage of
computing the GSH compared to the computation of the Hessian of the quadratic
interpolation model has been discussed in Remark 2.8. Another advantage of the
GSH compared to the simplex Hessian described in [17, 21] is that it provides
a simple explicit formula that is well-defined as long as the matrices of
directions employed are nonempty. Hence, the formula for the GSH provides an
approximation of the full Hessian, or a partial Hessian, for all possible
cases defined in Definition 2.9.

Next, we provide an example that serves to show that a set of $(n+1)(n+2)/2$
distinct points in $\operatorname{\mathbb{R}}^{n}$ that is poised for
quadratic interpolation is not necessarily a minimal poised set for the GSH.

######  Example 5.8.

Let $x^{0}=\begin{bmatrix}0&0\end{bmatrix}^{\top}$ be the point of interest.
Consider $\mathcal{X}=$ $\\{x^{0}=\mathbf{0},$ $e^{1},$ $e^{2},$ $-e^{1},$
$-e^{2},$ $-e^{1}-e^{2}\\}.$ Then $\mathcal{X}$ is poised for quadratic
interpolation, but cannot be expressed as a minimal poised set for the GSH at
$x^{0}$.

###### Proof.

Using a similar approach as in Proposition 5.7, one can verify that
$\operatorname{\mathcal{X}}$ is poised for quadratic interpolation. Now we
show that $\operatorname{\mathcal{X}}$ is not a minimal poised set for the GSH
at $x^{0}$ using brute force. We need to build $S=\\{s^{1},s^{2}\\}$ such that
the matrix corresponding to $S$ is full rank and $x^{0}\oplus
S\subseteq\operatorname{\mathcal{X}}$. Hence, the possible choices for $S$ are

$\displaystyle S$
$\displaystyle\in\left\\{\\{e^{1},e^{2}\\},\\{e^{1},-e^{2}\\},\\{e^{1},-e^{1}-e^{2}\\},\\{-e^{1},e^{2}\\},\right.$
$\displaystyle\quad\quad\left.\\{-e^{1}-e^{2},e^{2}\\},\\{-e^{1},-e^{2}\\},\\{-e^{1},-e^{1}-e^{2}\\},\\{-e^{1}-e^{2},-e^{2}\\}\right\\}.$

Case I: $S=\\{e^{1},e^{2}\\}.$ In this case, we need to build
$\overline{T}=\\{t^{1},t^{2}\\}$ such that the matrix corresponding to
$\overline{T}$ is full rank and

$\\{t^{1},e^{1}+t^{1},e^{2}+t^{1},t^{2},e^{1}+t^{2},e^{2}+t^{2}\\}=\operatorname{\mathcal{X}}.$

We see that the only possible choice of $t^{1}$ such that
$\\{t^{1},e^{1}+t^{1},e^{2}+t^{1}\\}\subseteq\operatorname{\mathcal{X}}$ is
$t^{1}=-e^{1}-e^{2}$ (note $t^{1}\neq\mathbf{0}$ as we require full rank).
However, the only possible choice of $t^{2}$ such that
$\\{t^{2},e^{1}+t^{2},e^{2}+t^{2}\\}\subseteq\operatorname{\mathcal{X}}$ is
$t^{2}=-e^{1}-e^{2}$. As full rank implies $t^{1}$ cannot equal $t^{2}$, we
see $S=\\{e^{1},e^{2}\\}$ cannot provide the desired properties. Cases II
through VIII: The other options for $S$ can be eliminated analogously.
Therefore, $\operatorname{\mathcal{X}}$ cannot be expressed as a minimal
poised set for the GSH at $x^{0}$. ∎

Figure 2 shows all possible directions connecting two points in the set
$\operatorname{\mathcal{X}}$ from Example 5.8. If $\operatorname{\mathcal{X}}$
were a minimal poised set for the GSH at $x^{0}=\mathbf{0}$, then it would be
possible to choose two directions (lines in blue) emerging from $x^{0}$ that
connect to other points in $\mathcal{Y}$ and these same two directions would
be emerging from two other points.

Figure 2: A set that is poised for quadratic interpolation but not a minimal
poised set for the GSH at $x^{0}$.

We conclude this section by providing formulae to obtain all the values of the
coefficients involved in the quadratic interpolation function of $f$ over a
minimal poised set for the GSH of the form
$\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})$, which we denote by
$Q_{f}(x^{0};S,U_{\ell})(x).$

######  Proposition 5.9.

Let
$f:\operatorname{dom}f\subseteq\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$.
Let $\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})\subset\operatorname{dom}f$
be a minimal poised set for the GSH at $x^{0}$ as constructed in Proposition
5.3. If $\ell\in\\{1,2,\dots,n\\},$ then the Hessian matrix $\mathbf{H}$ of
the quadratic interpolation function $Q_{f}(x^{0};S,U_{\ell})(x)$ is given by
$\mathbf{H}=S^{-\top}\widehat{\mathbf{H}}S^{-1},$ where the entries of the
symmetric matrix $\widehat{\mathbf{H}}\in\operatorname{\mathbb{R}}^{n\times
n}$ are

$\displaystyle\widehat{\mathbf{H}}_{i,\ell}$
$\displaystyle=-f(x^{0}+s^{i}-s^{\ell})+f(x^{0}+s^{i})+f(x^{0}-s^{\ell})-f(x^{0}),\quad
i\in\\{1,2,\dots,n\\}\setminus\\{\ell\\},$
$\displaystyle\widehat{\mathbf{H}}_{i,i}$
$\displaystyle=f(x^{0}+2s^{i}-s^{\ell})-2f(x^{0}+s^{i}-s^{\ell})+f(x^{0}-s^{\ell}),\quad
i\in\\{1,2,\dots,n\\}\setminus\\{\ell\\},$
$\displaystyle\widehat{\mathbf{H}}_{\ell,\ell}$
$\displaystyle=f(x^{0}+s^{\ell})+f(x^{0}-s^{\ell})-2f(x^{0}),$
$\displaystyle\widehat{\mathbf{H}}_{i,j}$
$\displaystyle=f(x^{0}+s^{i}+s^{j}-s^{\ell})-f(x^{0}+s^{i}-s^{\ell})-f(x^{0}+s^{j}-s^{\ell})+f(x^{0}-s^{\ell}),$

for all $i,j\in\\{1,2,\dots,n\\}\setminus\\{\ell\\},i\neq j.$ If $\ell=0,$
then

$\displaystyle\widehat{\mathbf{H}}_{i,i}$
$\displaystyle=f(x^{0}+2s^{i})-2f(x^{0}+s^{i})+f(x^{0}),\quad
i\in\\{1,2,\dots,n\\},$ $\displaystyle\widehat{\mathbf{H}}_{i,j}$
$\displaystyle=f(x^{0}+s^{i}+s^{j})-f(x^{0}+s^{i})-f(x^{0}+s^{j})+f(x^{0}),\quad
i,j\in\\{1,2,\dots,n\\},i\neq j.$

For all $\ell\in\\{0,1,\dots,n\\},$ the vector $\alpha$ associated to
$Q_{f}(x^{0};S,U_{\ell})(x)$ is given by $\alpha=S^{-\top}\overline{\alpha}$,
where the entries of $\overline{\alpha}\in\operatorname{\mathbb{R}}^{n}$ are

$\displaystyle\overline{\alpha}_{i}$
$\displaystyle=f(x^{0}+s^{i})-f(x^{0})-\frac{1}{2}\widehat{\mathbf{H}}_{i,i}-(x^{0})^{\top}\mathbf{H}s^{i},\quad
i\in\\{1,2,\dots,n\\}.$

The scalar $\alpha_{0}$ of $Q_{f}(x^{0};S,U_{\ell})(x)$ is

$\displaystyle\alpha_{0}$
$\displaystyle=f(x^{0})-\alpha^{\top}x^{0}-\frac{1}{2}(x^{0})^{\top}\mathbf{H}x^{0}.$

###### Proof.

The result is obtained by using Definition 2.2. Let

$\displaystyle
Q_{f}(x^{0};S,U_{\ell})(x)=\alpha_{0}+\alpha^{\top}x+\frac{1}{2}x^{\top}\mathbf{H}x,$

where
$a_{0}\in\operatorname{\mathbb{R}},\alpha\in\operatorname{\mathbb{R}}^{n}$ and
$\mathbf{H}=\mathbf{H}^{\top}\in\operatorname{\mathbb{R}}^{n\times n}.$
Suppose $\ell\in\\{1,2,\dots,n\\}.$ Evaluating $Q_{f}(x^{0};S,U_{\ell})(x)$ at
$x^{0}$, we obtain

$\displaystyle\alpha_{0}+\alpha^{\top}x^{0}+\frac{1}{2}(x^{0})^{\top}\mathbf{H}x^{0}$
$\displaystyle=f(x^{0}).$ (29)

Evaluating $Q_{f}(x^{0};S,U_{\ell})(x)$ at $x^{0}+s^{i}$ and using (29), we
obtain

$\displaystyle
f(x^{0})+\overline{\alpha}^{\top}e^{i}+(x^{0})^{\top}\overline{\mathbf{H}}e^{i}+\frac{1}{2}\widehat{\mathbf{H}}_{i,i}$
$\displaystyle=f(x^{0}+s^{i}),\quad i\in\\{1,\dots,n\\},$ (30)

where
$\overline{\alpha}^{\top}=\alpha^{\top}S,\overline{\mathbf{H}}=\mathbf{H}S$
and $\widehat{\mathbf{H}}=S^{\top}\mathbf{H}S.$ Evaluating
$Q_{f}(x^{0};S,U_{\ell})(x)$ at $x^{0}+s^{i}-s^{\ell}$ and using (29) and
(30), we obtain

$\displaystyle
f(x^{0}+s^{i})-\overline{\alpha}^{\top}s^{\ell}-(x^{0})^{\top}\overline{\mathbf{H}}e^{\ell}-\widehat{\mathbf{H}}_{i,\ell}+\frac{1}{2}\widehat{\mathbf{H}}_{\ell,\ell}=f(x^{0}+s^{i}-s^{\ell}),i\in\\{1,\dots,n\\}\setminus\\{\ell\\}.$

(31)

Evaluating $Q_{f}(x^{0};S,U_{\ell})(x)$ at $x^{0}-s^{\ell}$ and using (29), we
find

$\displaystyle-\overline{\alpha}^{\top}e^{\ell}-(x^{0})^{\top}\overline{\mathbf{H}}e^{\ell}+\frac{1}{2}\widehat{\mathbf{H}}_{\ell,\ell}$
$\displaystyle=f(x^{0}-s^{\ell})-f(x^{0}).$ (32)

Substituting (32) into (31), we obtain

$\widehat{\mathbf{H}}_{i,\ell}=\widehat{\mathbf{H}}_{\ell,i}=-f(x^{0}+s^{i}-s^{\ell})+f(x^{0}+s^{i})+f(x^{0}-s^{\ell})-f(x^{0}),i\in\\{1,\dots,n\\}\setminus\\{\ell\\}.$
(33)

Evaluating $Q_{f}(x^{0};S,U_{\ell})(x)$ at $x^{0}+2s^{i}-s^{\ell}$ and using
(29), (30) and (32), we find

(34)

Substituting (33) in (34), we find

$\displaystyle\widehat{\mathbf{H}}_{i,i}$
$\displaystyle=f(x^{0}+2s^{i}-s^{\ell})+f(x^{0}-s^{\ell})-2f(x^{0}+s^{i}-s^{\ell}),\quad
i\in\\{1,2,\dots,n\\}\setminus\\{\ell\\}.$

Evaluating $Q_{f}(x^{0};S,U_{\ell})(x)$ at $x^{0}+s^{i}+s^{j}-s^{\ell}$ and
using (29), (30) and (32), we find

$\displaystyle\widehat{\mathbf{H}}_{i,j}=\widehat{\mathbf{H}}_{j,i}$
$\displaystyle=f(x^{0}+s^{i}+s^{j}-s^{\ell})-f(x^{0}+s^{i}-s^{\ell})-f(x^{0}+s^{j}-s^{\ell})+f(x^{0}-s^{\ell}),$

for $i,j\in\\{1,2,\dots,n\\}\setminus\\{\ell\\},i\neq j.$ Rearranging (32), we
get

$\displaystyle\frac{1}{2}\widehat{\mathbf{H}}_{\ell,\ell}-f(x^{0}-s^{\ell})+f(x^{0})$
$\displaystyle=\overline{\alpha}^{\top}e^{\ell}+(x^{0})^{\top}\overline{\mathbf{H}}e^{\ell}.$
(35)

Substituting (35) into (30) for $i=\ell$, we obtain

$\displaystyle\widehat{\mathbf{H}}_{\ell,\ell}$
$\displaystyle=f(x^{0}+s^{\ell})+f(x^{0}-s^{\ell})-2f(x^{0}).$

The entries of the vector $\overline{\alpha}$ are found by isolating
$\overline{\alpha}^{\top}e^{i}$ in (30). We obtain

$\displaystyle\overline{\alpha}_{i}$
$\displaystyle=f(x^{0}+s^{i})-f(x^{0})-\frac{1}{2}\widehat{\mathbf{H}}_{i,i}-(x^{0})^{\top}\mathbf{H}s^{i},\quad
i\in\\{1,2,\dots,n\\}$

where $\alpha_{i}=S^{-\top}\overline{\alpha}_{i}.$ Lastly, the scalar
$\alpha_{0}$ is obtained from (29). We find

$\displaystyle\alpha_{0}$
$\displaystyle=f(x^{0})-\alpha^{\top}x^{0}-\frac{1}{2}(x^{0})^{\top}\mathbf{H}x^{0}.$

If $\ell=0$, a similar process can be applied to obtain $\mathbf{H},\alpha$
and $\alpha_{0}.$ ∎

Using $\overline{T}=U_{\ell}$, it is worth emphasizing that
$Q_{f}(x^{0};S,\overline{T})(x)$ can be obtained for free in terms of function
evaluations whenever $\nabla^{2}_{s}f(x^{0};S,\overline{T})$ has already been
computed. Indeed, all the coefficients of $Q_{f}(x^{0};S,\overline{T})$ are
computed using the same function evaluations used for finding
$\nabla_{s}^{2}f(x^{0};S,\overline{T}).$ Since
$\mathcal{M}(x^{0};S,\overline{T})$ is poised for quadratic interpolation,
there exists one and only one Hessian matrix $\mathbf{H}$. It follows that the
Hessian $\mathbf{H}$ of the quadratic interpolation function
$Q_{f}(x^{0};S,\overline{T})(x)$ must be equal to
$\nabla_{s}^{2}f(x^{0};S,\overline{T}).$

##  6 Conclusion

We have presented the generalized simplex Hessian (GSH), a compact and simple
formula for approximating the Hessian of a function $f$. The GSH is well-
defined as long as the matrices $S$ and $T_{i}$ are nonempty. In particular,
the matrices $S$ and $T_{i}$ do not need to be full row rank nor square to
compute the GSH. We have also presented a centered version of the GSH called
the generalized centered simplex Hessian (GCSH) and shown its relation to the
GSH (Section 3).

In Theorems 4.2 and 4.3, we developed error bounds for the GSH and GCSH. Under
some assumptions , the GSH is an order-1 accurate approximation of the partial
Hessian $\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$, and the GCSH is
an order-2 accurate approximation of the partial Hessian
$\operatorname{Proj}_{S,T_{1:m}}\nabla^{2}f(x^{0})$. Indeed, when all matrices
of directions in the set $T_{1:m}$ are equal, then we obtain error bounds that
covers all possible 16 cases defined in Definition 2.9. When all matrices
$T_{i}$ are not equal, then 12 out of 16 cases are covered. Future research
directions include an investigation of the cases not covered in Theorems 4.2
and 4.3: the four cases $S$-nondetermined \overdetermined and
$T_{1:m}$-nondetermined\underdetermined.

We note that when a matrix is nondetermined or overdetermined, it is possible
to remove columns (or rows) of the matrix to make it full rank. This may lead
to error bounds for the nondetermined cases.

As the results in this paper are reported in terms of order-$N$ accuracy, it
is possible to apply the theory of DFO calculus developed in [13] directly. In
[13], approximation based product, quotient, and chain rules are developed.
For each rule, it is shown that when an objective function can be decomposed
into a product, quotient, or composition of simpler functions, then it is
possible to construct approximate Hessians of the simpler functions and then
use those to construct approximation Hessians of the objective function.
Combining [13] with the research herein has lead to a model-based trust region
algorithm using the calculus-based approximation techniques [26].

In Section 5, we showed that if a minimal poised set for the GSH is used, the
number of distinct function evaluations required to compute the GSH is
$(n+1)(n+2)/2.$ We investigated the relationship between a minimal poised set
for the GSH and poisedness for quadratic interpolation, proving that a minimal
poised set for the GSH of the form
$\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})$ is well-poised for quadratic
interpolation, but that the converse does not necessarily hold. In this case,
we know that $\nabla^{2}_{s}f(x^{0};S,\overline{T})$ with
$\overline{T}=U_{\ell}$ is equal to the Hessian of the quadratic interpolation
function $Q_{f}(x^{0};S,U_{\ell})$. We also developed explicit formulae for
obtaining the parameters of the quadratic interpolation function of $f$ over a
minimal poised set for the GSH of the form
$\operatorname{\mathcal{M}}(x^{0};S,U_{\ell})$.

Future research should explore if the sets $U_{0},U_{1},\dots,U_{n}$ as
defined in Proposition 5.3 are the only possible choices for
$\overline{T}\in\operatorname{\mathbb{R}}^{n\times n}$ such that
$\operatorname{\mathcal{S}}(x^{0};S,\overline{T}),$ where
$S\in\operatorname{\mathbb{R}}^{n\times n}$ is full rank, is a minimal poised
set for GSH at $x^{0}$. It can be proven that it is indeed the case by using
brute force in $\operatorname{\mathbb{R}}$ and
$\operatorname{\mathbb{R}}^{2},$ but it is still unclear how to generalized
this claim in an arbitrary dimension $n$.

While it is immediately clear that a minimal poised set for GSH can be used to
construct a ‘small’ set for the GCSH (using Proposition 3.3), it is not clear
if such a set would be minimal. Future research will investigate the
properties of a minimal poised set for the GCSH and the matrices of directions
that can be used to minimize the number of sample points.

It is known that the diagonal entries of the Hessian may be obtained for free
in terms of function evaluations when a centered simplex gradient is computed.
Hence, a future research direction could be to verify if the main diagonal
entries of $\nabla^{3}f(x^{0})$ may be obtained for free when a GCSH is
computed.

On a final note, MATLAB implementations of the GSH and the GCSH are available
upon request.

## Acknowledgements

The authors express their deep gratitude to the anonymous referees for the
comprehensive reviews that improved the quality of this work.

Hare’s research is partially funded by the Natural Sciences and Engineering
Research Council (NSERC) of Canada, Discover Grant #2018-03865. Jarry-Bolduc’s
research is partially funded by the Natural Sciences and Engineering Research
Council (NSERC) of Canada, Discover Grant #2018-03865. Jarry-Bolduc would like
to acknowledge UBC for the funding received through the University Graduate
Fellowship award.

## References

  * [1] N. Amaioua, C. Audet, A. Conn, and S. Le Digabel.  Efficient solution of quadratically constrained quadratic subproblems within the mesh adaptive direct search algorithm.  European Journal of Operational Research, 268(1):13–24, 2018.
  * [2] C. Audet, V. Béchard, and S. Le Digabel.  Nonsmooth optimization through mesh adaptive direct search and variable neighborhood search.  Journal of Global Optimization, 41(2):299–318, 2008.
  * [3] C. Audet and W. Hare.  Derivative-free and Blackbox Optimization.  Springer, 2017.
  * [4] C. Audet and W. Hare.  Algorithmic construction of the subdifferential from directional derivatives.  Set-Valued and Variational Analysis, 26(3):431–447, 2018.
  * [5] C. Audet, A. Ianni, S. Le Digabel, and C. Tribes.  Reducing the number of function evaluations in mesh adaptive direct search algorithms.  SIAM Journal on Optimization, 24(2):621–642, 2014.
  * [6] C. Audet and C. Tribes.  Mesh-based Nelder–Mead algorithm for inequality constrained optimization.  Computational Optimization and Applications, 71(2):331–352, 2018\.
  * [7] A. Bagirov, B. Karasözen, and M. Sezer.  Discrete gradient method: derivative-free method for nonsmooth optimization.  Journal of Optimization Theory and Applications, 137(2):317–334, 2008.
  * [8] A. Berahas, R. Byrd, and J. Nocedal.  Derivative-free optimization of noisy functions via quasi-Newton methods.  SIAM Journal on Optimization, 2019.  To appear.
  * [9] F. Berghen and H. Bersini.  CONDOR, a new parallel, constrained extension of Powell’s UOBYQA algorithm: Experimental results and comparison with the DFO algorithm.  Journal of Computational and Applied Mathematics, 181(1):157–175, 2005.
  * [10] D. Bortz and C. Kelley.  The simplex gradient and noisy optimization problems.  In Computational Methods for Optimal Design and Control, pages 77–90. Springer, 1998.
  * [11] P. Braun, W. Hare, and G. Jarry-Bolduc.  Limiting behavior of derivative approximation techniques as the number of points tends to infinity on a fixed interval in r.  Journal of Computational and Applied Mathematics, 386:113218, 2021\.
  * [12] R. Burden, J. Faires, and A. Burden.  Numerical analysis 10/e ie.  Brooks/Cole Cengage Learning, 2016.
  * [13] Y. Chen, W. Hare, and G. Jarry-Bolduc.  Error analysis of surrogate models constructed through operations on submodels.  Mathematics of Operations Research, 2022.
  * [14] G. Cocchi, G. Liuzzi, A. Papini, and M. Sciandrone.  An implicit filtering algorithm for derivative-free multiobjective optimization with box constraints.  Computational Optimization and Applications, 69(2):267–296, 2018\.
  * [15] A. Conn, K. Scheinberg, and L. Vicente.  Geometry of interpolation sets in derivative free optimization.  Mathematical Programming, 111(1-2):141–172, 2008.
  * [16] A. Conn, K. Scheinberg, and L. Vicente.  Geometry of sample sets in derivative-free optimization: polynomial regression and underdetermined interpolation.  IMA journal of numerical analysis, 28(4):721–748, 2008.
  * [17] A. Conn, K. Scheinberg, and L. Vicente.  Introduction to derivative-free optimization.  SIAM, 2009.
  * [18] I. Coope and R. Tappenden.  Efficient calculation of regular simplex gradients.  Computational Optimization and Applications, 72(3):561–588, 2019\.
  * [19] I. Coope and R. Tappenden.  Gradient and diagonal hessian approximations using quadratic interpolation models and aligned regular bases.  Numerical Algorithms, 88:767–791, 2021.
  * [20] A. Custódio, J. Dennis, and L. Vicente.  Using simplex gradients of nonsmooth functions in direct search methods.  IMA Journal of Numerical Analysis, 28(4):770–784, 2008.
  * [21] A. Custódio and L. Vicente.  Using sampling and simplex derivatives in pattern search methods.  SIAM Journal on Optimization, 18(2):537–555, 2007.
  * [22] S. Gratton, C. Royer, and L. Vicente.  A decoupled first/second-order steps technique for nonconvex nonlinear unconstrained optimization with improved complexity bounds.  Mathematical Programming, 179:195–222, 2020.
  * [23] S. Gratton, C. Royer, L. Vicente, and Z. Zhang.  Direct search based on probabilistic feasible descent for bound and linearly constrained problems.  Computational Optimization and Applications, 72(3):525–559, 2019\.
  * [24] W. Hare.  A discussion on variational analysis in derivative-free optimization.  Set-Valued and Variational Analysis, pages 1–17, 2020.
  * [25] W. Hare and G. Jarry-Bolduc.  Calculus identities for generalized simplex gradients: rules and applications.  SIAM Journal on Optimization, 30(1):853–884, 2020.
  * [26] W. Hare and G. Jarry-Bolduc.  About the performance of a calculus-based approach to building model functions in a derivative-free trust-region algorithm.  Algorithms, 16(2), 2023.
  * [27] W. Hare, G. Jarry-Bolduc, and C. Planiden.  Limiting behaviour of the generalized simplex gradient as the number of points tends to infinity on a fixed shape in irn.  Set-Valued and Variational Analysis, 31(1):1–33, 2023.
  * [28] W. Hare, G. Jarry–Bolduc, and C. Planiden.  Error bounds for overdetermined and underdetermined generalized centred simplex gradients.  IMA Journal of Numerical Analysis, 42(1):744–770, 12 2020.
  * [29] W. Hare and J. Nutini.  A derivative-free approximate gradient sampling algorithm for finite minimax problems.  Computational Optimization and Applications, 56(1):1–38, 2013.
  * [30] W. Hare, C. Planiden, and C. Sagastizábal.  A derivative-free VU-algorithm for convex finite-max problems.  Optimization Methods and Software, 35(3):521–559, 2020.
  * [31] W. Hare, C. Sagastizábal, and M. Solodov.  A proximal bundle method for nonsmooth nonconvex functions with inexact information.  Computational Optimization and Applications, 63(1):1–28, 2016.
  * [32] G. Jarry-Bolduc.  Approximating the diagonal of a hessian: which sample set of points should be used.  Numerical Algorithms, 91(3):1349–1361, 2022.
  * [33] J. Larson, M. Menickelly, and S. Wild.  Manifold sampling for $l_{1}$ nonconvex optimization.  SIAM Journal on Optimization, 26(4):2540–2563, 2016.
  * [34] G. Liuzzi, S. Lucidi, F. Rinaldi, and L. Vicente.  Trust-region methods for the derivative-free optimization of nonsmooth black-box functions.  SIAM Journal on Optimization, 29(4):3012–3035, 2019.
  * [35] A. Maggiar, A. Wächter, I. Dolinskaya, and J. Staum.  A derivative-free trust-region algorithm for the optimization of functions smoothed via Gaussian convolution using adaptive multiple importance sampling.  SIAM Journal on Optimization, 28(2):1478–1507, 2018.
  * [36] M. Menickelly and S. Wild.  Derivative-free robust optimization by outer approximations.  Mathematical Programming, 2019.  To appear.
  * [37] J. Nocedal and S. Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [38] M. Powell.  On the lagrange functions of quadratic models that are defined by interpolation.  Optimization Methods and Software, 16(1-4):289–309, 2001.
  * [39] M. Powell.  On trust region methods for unconstrained minimization without derivatives.  Mathematical Programming, 97(3):605–623, 2003.
  * [40] M. Powell.  The BOBYQA algorithm for bound constrained optimization without derivatives.  Cambridge NA Report NA2009/06, University of Cambridge, Cambridge, pages 26–46, 2009.
  * [41] R. Regis.  The calculus of simplex gradients.  Optimization Letters, 9(5):845–865, 2015.
  * [42] R. Rockafellar and R. Wets.  Variational analysis.  Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, 1998.
  * [43] S. Shashaani, F. Hashemi, and R. Pasupathy.  ASTRO-DF: A class of adaptive sampling trust-region algorithms for derivative-free stochastic optimization.  SIAM Journal on Optimization, 28(4):3145–3176, 2018.
  * [44] A. Verdério, E. Karas, L. Pedroso, and K. Scheinberg.  On the construction of quadratic models for derivative-free trust-region algorithms.  EURO Journal on Computational Optimization, 5(4):501–527, 2017\.
  * [45] S. Wild and C. Shoemaker.  Global convergence of radial basis function trust-region algorithms for derivative-free optimization.  SIAM REVIEW, 55(2):349–371, 2013.