# QUASI-NEWTON METHODS FOR MINIMIZING A QUADRATIC FUNCTION SUBJECT TO
UNCERTAINTY

Shen PENG11footnotemark: 1 Gianpiero CANESSA11footnotemark: 1 David
EK11footnotemark: 1 Anders FORSGREN

(Draft )

###### Abstract

We investigate quasi-Newton methods for minimizing a strictly convex quadratic
function which is subject to errors in the evaluation of the gradients. The
methods all give identical behavior in exact arithmetic, generating minimizers
of Krylov subspaces of increasing dimensions, thereby having finite
termination. A BFGS quasi-Newton method is empirically known to behave very
well on a quadratic problem subject to small errors. We also investigate
large-error scenarios, in which the expected behavior is not so clear. In
particular, we are interested in the behavior of quasi-Newton matrices that
differ from the identity by a low-rank matrix, such as a memoryless BFGS
method. Our numerical results indicate that for large errors, a memory-less
quasi-Newton method often outperforms a BFGS method. We also consider a more
advanced model for generating search directions, based on solving a chance-
constrained optimization problem. Our results indicate that such a model often
gives a slight advantage in final accuracy, although the computational cost is
significantly higher.

Keywords: Quadratic programming, quasi-Newton method, stochastic quasi-Newton
method, chance constrained model

##  1 Introduction

A strictly convex $n$-dimensional quadratic function may be written on the
form

$q(x)={\textstyle{\frac{1}{2}}}x^{T}\\!Hx+c^{T}\\!x+d,$

where $H$ is a positive definite and symmetric $n\times n$-matrix, $c$ is an
$n$-dimensional vector and $d$ is a constant. The optimization problem of
minimizing $q(x)$ is equivalent to solving $\nabla q(x)=0$, i.e., solving the
linear equation $Hx+c=0$.

One way to do so by a direct method is to find an initial point $x_{0}$ and
associated gradient $g_{0}=Hx_{0}+c$. Then generate $x_{k}$ and $g_{k}$, with
$g_{k}=Hx_{k}+c$, such that $x_{k}$ is the minimizer of $q(x)$ on
$x_{0}+{\mathcal{K}}_{k}(g_{0},H)$, where

${\mathcal{K}}_{0}(g_{0},H)=\\{0\\},\quad{\mathcal{K}}_{k}(g_{0},H)=\mathop{\operator@font
span}\\{g_{0},Hg_{0},H^{2}g_{0},\dots,H^{k-1}g_{0}\\},\quad k=1,2,\dots.$

This is equivalent to $g_{k}$ being orthogonal to
${\mathcal{K}}_{k}(g_{0},H)$, so that $g_{0}$, $g_{1}$, …, $g_{k}$ form an
orthogonal basis for ${\mathcal{K}}_{k+1}(g_{0},H)$. Since there can be at
most $n$ nonzero orthogonal vectors, there is an $r$, with $r\leq n$, such
that $g_{r}=0$. Consequently, $x_{r}$ is the minimizer of $q(x)$.

A method for computing $x_{1}$, $x_{2}$, …, $x_{r}$ this way is characterized
by the search direction $p_{k}$ leading from $x_{k}$ to $x_{k+1}$. Given
$p_{k}$, the step length $\alpha_{k}$ is given by minimizing $q(x_{k}+\alpha
p_{k})$, i.e.,

$\alpha_{k}=-\frac{g_{k}^{T}p_{k}}{p_{k}^{T}Hp_{k}}.$ (1.1)

Therefore, it suffices to characterize $p_{k}$. One characterization is given
by the conditions that

$\displaystyle(i)\quad$ $p_{k}$ is a linear combination of $g_{0}$, $g_{1}$,
$\dots$, $g_{k}$, in addition to (1.2a) $\displaystyle(ii)\quad$ satisfying
$g_{i}^{T}p_{k}=-g_{k}^{T}g_{k}$, $i=0,\dots,k,$ (1.2b)

see, e.g., [EF21, Lemma 1]. There could be an arbitrary scaling
$g_{i}^{T}p_{k}=c_{k}$ for a nonzero scalar $c_{k}$. Throughout, we will use
$c_{k}=-g_{k}^{T}g_{k}$.

The method of conjugate gradients gives a short recursion for the search
direction $p_{k}$ satisfying (1.2). It may be written on the form

$p_{k}=\begin{cases}-g_{0},&k=0,\\\
-g_{k}+\frac{g_{k}^{T}g_{k}}{g_{k-1}^{T}g_{k-1}}p_{k-1},&k=1,2,\dots,r-1.\end{cases}$
(1.3)

See, e.g., [Saa03, Chapter 6] for an introduction to the method of conjugate
gradients.

An alternative way of computing the search direction $p_{k}$ satisfying (1.2)
is through a quasi-Newton method, in which $p_{k}$ is defined by a linear
equation $B_{k}p_{k}=-g_{k}$, for $B_{k}$ a symmetric positive-definite
matrix. A well-known method for which $B_{k}$ gives $p_{k}$ satisfying (1.2)
is the BFGS quasi-Newton method, in which $B_{0}=I$, and $B_{k}$ is formed by
adding a symmetric rank-2 matrix to $B_{k-1}$. In our setting, the BFGS method
may be viewed as the “ideal” update that dynamically transforms $B_{k}$ from
$I$ to $H$ in $n$ steps. Identity curvature is transformed to $H$ curvature in
one dimension at each step, and this information is maintained throughout. The
precise formulas will be given in Section 2.

In exact arithmetic, the method of conjugate gradients and a BFGS quasi-Newton
method compute identical iterates in the setting we consider, minimizing a
strictly convex quadratic problem using exact linesearch. In this situation,
the recursion formula for the method of conjugate gradients is to prefer,
since the computational cost for solving with $B_{k}$ for the BFGS method
increases with $k$.

In finite precision arithmetic, the BFGS quasi-Newton method may still be
expected to compute search directions of very high quality as the Hessian is
approximated on subspaces of higher dimension. This is not to be expected for
the method of conjugate gradients. Our interest is to study the situation
where noise is added to the gradients, thereby considering a situation with
significantly higher level of error than finite precision arithmetic. In this
situation, it is not clear that the BFGS quasi-Newton method is superior, in
the sense that the Hessian approximation may become inaccurate. It is here
also of interest to consider the method of steepest descent, where the search
direction is the negative of the gradient. For exact arithmetic, the search
directions satisfying (1.2) will outperform steepest descent due to the
property of minimizing the quadratic objective function over expanding
subspaces. In the case of large noise, this is not at all clear.

The reason for noise in the gradients can be seen in different perspectives.
Firstly, as mentioned above, the finite precision arithmetic gives a residual
between the evaluated gradients and the true gradients. Secondly, in many
practical problem, such as PDE-constrained optimization, the objective
function often contains computational noise created by an inexact linear
system solver, adaptive grids, or other internal computations. In other cases,
the noise in the gradients can be due to stochastic errors. For example, when
minimizing $q(x)=\mathbb{E}[f(x;\xi)]$, where $\xi$ is a random variable. With
given sample set $\Xi=\\{\xi^{i},i=1,\cdots,N\\}$, instead of $q(x)$, the
following empirical expectation will be considered:

$\tilde{q}(x)=\frac{1}{N}\sum_{\xi^{i}\in\Xi}f(x;\xi^{i}).$

Due to the randomness of samples, $\tilde{q}(x)=q(x)+\varepsilon$, and
$\nabla\tilde{q}(x)=\nabla q(x)+\epsilon$, where $\varepsilon,\epsilon$ are
random noise.

The basis for the methods we consider is that they compute search directions
identical to those of the BFGS method and the method of conjugate gradients,
i.e., satisfying (1.2), in exact arithmetic. In particular, we are interested
in a setting where $B_{k}$ is the identity matrix plus a symmetric matrix of
rank two. We will refer to such a quasi-Newton matrix as a _low-rank_ quasi-
Newton matrix and name the corresponding method a low-rank quasi-Newton
method. Our intention is to investigate the behavior of a low-rank quasi-
Newton method compared to the method of steepest descent, thereby mimicking a
situation where two more vectors in addition to $g_{k}$ are used at iteration
$k$. The corresponding search direction can then be computed from a two-by-two
system. For comparison, we also compare to the BFGS quasi-Newton method. Our
choice of quadratic problem allows an environment where the behavior of the
methods in infinite precision is known, and we can study the effect of noise.

In addition, we investigate the potential for improving performance of the
quasi-Newton method by formulating robust optimization problems of chance-
constraint type for computing the search directions. These methods become of
higher interest in the case of large noise and multiple copies of the
gradients. Our interest is to capture the essence of the behavior, and try to
understand the interplay between quality in computed direction compared to
robustness given by the chance constraints. The computational cost will always
be significantly higher, but our interest is to see if we can gain in terms of
robustness and accuracy of the computed solution.

###  1.1 Background and related work

The paper builds on previous work in the setting of exact arithmetic and
finite precision arithmetic. Forsgren and Odland [FO18] have studied exact
linesearch quasi-Newton methods for minimizing a strictly convex quadratic
function, and given necessary and sufficient conditions for a quasi-Newton
matrix $B_{k}$ to generate a search direction which is parallel to that of
(1.2) in exact arithmetic. With exact linesearch methods, Ek and Forsgren
[EF21] have studied certain limited-memory quasi-Newton Hessian approximations
for minimizing a convex quadratic function in the setting of finite precision
arithmetic. Dennis and Walker [DW84] have considered the use of bounded-
deterioration quasi-Newton methods implemented in floating-point arithmetic
where only inaccurate values are available. In contrast, our work allows for
large noise and we study performance on a set of test problems.

In the present manuscript, we consider a situation where the function values
and gradients cannot be easily obtained and only noisy information about the
gradient is available. To handle this situation, some stochastic methods are
proposed to minimize the objective function with inaccurate information. Our
setting is minimizing a strictly convex quadratic function.

For strongly convex problems, Mokhtari and Ribeiro [MR13] have proposed a
regularized stochastic BFGS method and analyzed its convergence, and Mokhtari
and Ribeiro [MR15] have further studied an online L-BFGS method. Berahas,
Nocedal and Takac [BNT16] have considered the stable performance of quasi-
Newton updating in the multi-batch setting, illustrated the behavior of the
algorithm and studied its convergence properties for both the convex and
nonconvex cases. Byrd et al.  [BHNS16] have proposed a stochastic quasi-Newton
method in limited memory form through subsampled Hessian-vector products. Shi
et al.  [SXBN20] have proposed practical extensions of the BFGS and L-BFGS
methods for nonlinear optimization that are capable of dealing with noise by
employing a new linesearch technique. More recently, Xie, Byrd and Nocedal
[XBN20] have considered the convergence analysis of quasi-Newton methods when
there are (bounded) errors in both function and gradient evaluations, and
established conditions under which an Armijo-Wolfe linesearch on the noisy
function yields sufficient decrease in the true objective function.

Unlike the stochastic quasi-Newton methods, which are based on the subsampled
gradients or Hessians, there are also other stochastic tools to reduce the
effect of noise when generating the search direction. Lucchi et al.  [LMH15]
have studied quasi-Newton method by incorporating variance reduction technique
to reduce the effect of noise in Hessian matrices by proposing a variance-
reduced stochastic Newton method. This method keeps the variance under control
in the use of a multi-stage scheme. Moritz et al.  [MNJ16] have proposed a
linearly convergent method that integrates the L-BFGS method to alleviate the
effect of noisy gradients with the variance reduction technique by adding the
residual between subsample gradient and full gradient to the noisy gradient.

In addition, chance constraint is a natural approach to handle the effect of
random noise [AX18]. Therefore, chance constraint has the potential to reduce
the effect of noise when generating the search direction. By integrating
chance constraints in the design of quasi-Newton methods, we investigate the
ability to improve robustness into the computation of the search direction in
the presence of noise.

##  2 Suggestions on quasi-Newton matrices

As mentioned in Section 1, a well-known method for which $B_{k}$ gives $p_{k}$
satisfying (1.2) is the BFGS quasi-Newton method. In the BFGS method,
$B_{0}=I$ and

$\displaystyle B_{k}$
$\displaystyle=B_{k-1}+\frac{1}{g_{k-1}^{T}p_{k-1}}g_{k-1}g_{k-1}^{T}$
$\displaystyle+\frac{1}{\alpha_{k-1}(g_{k}-g_{k-1})^{T}p_{k-1}}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T},\quad
k=1,\dots,r.$ (2.1)

Expansion gives

$\displaystyle B_{k}$
$\displaystyle=B_{k-1}+\frac{1}{g_{k-1}^{T}p_{k-1}}g_{k-1}g_{k-1}^{T}+\frac{1}{\alpha_{k-1}(g_{k}-g_{k-1})^{T}p_{k-1}}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T}$
$\displaystyle=I+\sum_{i=0}^{k-1}\frac{1}{g_{i}^{T}\\!p_{i}}g_{i}g_{i}^{T}+\sum_{i=0}^{k-1}\frac{1}{\alpha_{l}(g_{i+1}-g_{i})^{T}p_{i}}(g_{i+1}-g_{i})(g_{i+1}-g_{i})^{T}.$
(2.2)

For the quadratic case with exact linesearch, the BFGS matrix of (2.2) takes
the form

$\displaystyle B_{k}$
$\displaystyle=B_{k-1}-\frac{1}{g_{k-1}^{T}g_{k-1}}g_{k-1}g_{k-1}^{T}+\frac{1}{p_{k-1}^{T}Hp_{k-1}}Hp_{k-1}p_{k-1}^{T}H$
$\displaystyle=I-\sum_{i=0}^{k-1}\frac{1}{g_{i}^{T}\\!g_{i}}g_{i}g_{i}^{T}+\sum_{i=0}^{k-1}\frac{1}{p_{i}^{T}Hp_{i}}Hp_{i}p_{i}^{T}H,$
(2.3)

see, e.g., [EF21]. If $n$ steps are taken, then with $P_{n}=(p_{0}\ p_{1}\
\cdots\ p_{n-1})$, it follows that $P_{n}$ is square and nonsingular, so that

$\displaystyle H$
$\displaystyle=HH^{-1}H=HP_{n}P_{n}^{-1}H^{-1}P_{n}^{-T}\\!P_{n}^{T}H$
$\displaystyle=HP_{n}(P_{n}^{T}HP_{n})^{-1}P_{n}^{T}H=\sum_{i=0}^{n-1}\frac{1}{p_{i}^{T}Hp_{i}}Hp_{i}p_{i}^{T}H,$
(2.4)

where the conjugacy of the $p_{i}$s is a consequence of (1.2). In this
situation, (2.3) may therefore be seen as a dynamic way of generating the true
Hessian in $n$ steps, if the method does not converge early, as $B_{n}=H$ by

$I-\sum_{i=0}^{n-1}\frac{1}{g_{i}^{T}\\!g_{i}}g_{i}g_{i}^{T}=0\text{and}\sum_{i=0}^{n-1}\frac{1}{p_{i}^{T}Hp_{i}}Hp_{i}p_{i}^{T}H=H.$

This is a consequence of the orthogonal gradients then spanning the whole
space in combination with (2.4). Consequently, the BFGS method may be viewed
as the “ideal” update that dynamically transforms $B_{k}$ from $I$ to $H$ in
$n$ steps. Identity curvature is transformed to $H$ curvature in one dimension
at each step, and this curvature information is maintained throughout.

The discussion above may be generalized to a general class of quasi-Newton
matrices $B_{k}$ of the form

$B_{k}=V_{k}+\sum_{i=0}^{k-1}\rho_{i}(g_{i+1}-g_{i})(g_{i+1}-g_{i})^{T},$
(2.5)

where $V_{k}p_{k}=-g_{k}$ for $p_{k}$ satisfying (1.2), and $\rho_{i}$,
$i=0,\dots,k-1$, are nonnegative scalars. In exact arithmetic and under exact
linesearch, the specific values of $\rho_{i}$, $i=0,\dots,k-1$, have no impact
on the search direction, due to (1.2). In a noisy framework, they will make a
difference, and we will pay attention to how they are selected.

As discussed in Section 1, we are particularly interested in low-rank quasi-
Newton matrices. We will therefore consider a _memoryless BFGS_ quasi-Newton
method, in which

$V_{k}=I-\frac{1}{p_{k-1}^{T}p_{k-1}}p_{k-1}p_{k-1}^{T},$ (2.6)

in addition to $\rho_{i}=0$, $i=0,\dots,k-2$, so that

$B_{k}=I-\frac{1}{p_{k-1}^{T}p_{k-1}}p_{k-1}p_{k-1}^{T}+\rho_{k-1}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T},$
(2.7)

where the value of $\rho_{k-1}$ is given by the _secant condition_
$\alpha_{k-1}B_{k}p_{k-1}=g_{k}-g_{k-1}$. We denote this particular value of
$\rho_{k-1}$ by $\hat{\rho}_{k-1}$. For exact linesearch, $B_{k}$ of (2.7)
gives

$\hat{\rho}_{k-1}=-\frac{1}{\alpha_{k-1}g_{k-1}^{T}p_{k-1}}.$ (2.8)

Then $V_{k}p_{k}=-g_{k}$ for $p_{k}$ satisfying (1.2) in the case of exact
arithmetic, see, e.g., [FO18, Proposition 1]. The $V_{k}$ of the memoryless
BFGS matrix given by (2.6) is analogous to the first two terms of the BFGS
matrix of (2.3), but the matrix is singular with its nullspace restricted to
the one-dimensional span of $p_{k-1}$, as opposed to the span of all previous
gradients. In addition, the memoryless BFGS matrix $B_{k}$ of (2.7) is
nonsingular as $\rho_{k-1}>0$ and $(g_{k}-g_{k-1})^{T}p_{k-1}\neq 0$.

We will also interpret the method of conjugate gradients in a quasi-Newton
framework, by forming the _symmetric CG_ quasi-Newton matrix $B_{k}$ given by

$B_{k}=\left(I-\frac{1}{g_{k-1}^{T}p_{k-1}}g_{k}p_{k-1}^{T}\right)\left(I-\frac{1}{g_{k-1}^{T}p_{k-1}}p_{k-1}g_{k}^{T}\right).$
(2.9)

This matrix is formed by rewriting the recursion (1.3) and making an
additional symmetrization, see [FO18]. It can be put in the matrix family
given by (2.5) by setting $V_{k}=B_{k}$ and $\rho_{i}=0$, $i=0,\dots,k-1$.

In summary, we will consider quasi-Newton matrices of the form (2.5). In
particular, we will consider two specific low-rank matrices. The quasi-Newton
matrices $B_{k}$ given by _memoryless BFGS_ in (2.7) and _symmetric CG_ in
(2.9) are both low-rank quasi-Newton matrices that differ from the identity by
a symmetric rank two matrix and fulfill the conditions we require, (1.2). They
will be used in our computational study.

##  3 A chance-constrained model for finding the search direction

In addition to investigating the behavior of the quasi-Newton methods
discussed so far, we are also interested in investigating the potential of
increasing the performance of the quasi-Newton methods in the presence of
noise by selecting parameters from a chance-constrained optimization problem.

The aim is to design a quasi-Newton matrix $B_{k}$, with $B_{k}\succ 0$, so as
to compute a search direction of high quality. For the case of exact
arithmetic, i.e., no noise, our model direction is the direction $p_{k}$ that
satisfies (1.2). The interest is now to investigate and design quasi-Newton
matrices in the presence of noise. In particular, we are interested in
studying the performance of different methods for different noise levels. For
a given quasi-Newton matrix $B_{k}$, the search direction $p_{k}$ is computed
from $B_{k}p_{k}=-g_{k}$.

As there exists noise in each iteration, it means that the obtained gradient
is not accurate, it is the combination of the true gradient and some white
noise. The update of search direction may result in a non-descent direction
because of the influence by the noise. Then, we have the following proposition
to show that the search direction satisfying $B_{k}p_{k}=-g_{k}$ is a descent
direction under certain conditions even with noise in each iteration.

###### Proposition 3.1

Consider iteration $k$ of a quasi-Newton method for minimizing $q(x)$. Let
$g_{k}=\bar{g}_{k}+\epsilon$, where $\bar{g}_{k}$ is the true gradient and
$\epsilon$ is the noise generated with mean equal to $0$. If $B_{k}\succ 0$
and $\|\epsilon\|<\frac{1}{\|B_{k}^{-1}g_{k}\|}g_{k}^{T}B_{k}^{-1}g_{k}$, the
direction $p_{k}$, satisfying $B_{k}p_{k}=-g_{k}$, is a descent direction at
point $x_{k}$.

Proof. The direction $p_{k}$ is a descent direction if
$\bar{g}_{k}^{T}p_{k}<0$. Since $B_{k}p_{k}=-g_{k}$, we have

$\bar{g}_{k}^{T}p_{k}=(g_{k}-\epsilon)^{T}p_{k}=-(g_{k}-\epsilon)^{T}B_{k}^{-1}g_{k}=-g_{k}^{T}B_{k}^{-1}g_{k}+\epsilon^{T}B_{k}^{-1}g_{k}.$

Hence, $p_{k}$ is a descent direction if
$g_{k}^{T}B_{k}^{-1}g_{k}-\epsilon^{T}B_{k}^{-1}g_{k}>0$.

As $g_{k}^{T}B_{k}^{-1}g_{k}-\epsilon^{T}B_{k}^{-1}g_{k}\geq
g_{k}^{T}B_{k}^{-1}g_{k}-\|\epsilon\|\|B_{k}^{-1}g_{k}\|$,
$\|\epsilon\|<\frac{1}{\|B_{k}^{-1}g_{k}\|}g_{k}^{T}B_{k}^{-1}g_{k}$ implies
$p_{k}$ is a descent direction. This concludes the proof.

A consequence is that the property of $p_{k}$ being a descent direction may be
lost if the termination criteria on $\|\tilde{g}_{k}\|$ is set smaller than
the noise level, as observed in the following remark.

###### Remark 1

Proposition 3.1 shows that when the noise $\epsilon$ is small enough, the
direction satisfying $B_{k}p_{k}=-g_{k}$ is a descent direction. However, if
the noise $\epsilon$ is large, the direction obtained from $B_{k}p_{k}=-g_{k}$
may not be a descent direction. In particular, let $\tau>0$ denote the
tolerance level of the stopping criterion based on $\|g_{k}\|$. Consider the
situation when $\|\epsilon\|>\tau$. If $x_{k}$ is close to termination
solution, the value of $\|g_{k}\|$ is close to $\tau$. In this case, we could
have $\tau<\|g_{k}\|\leq\|\epsilon\|$. Since
$\frac{1}{\|B_{k}^{-1}g_{k}\|}g_{k}^{T}B_{k}^{-1}g_{k}\leq\frac{1}{\|B_{k}^{-1}g_{k}\|}\|g_{k}\|\|B_{k}^{-1}g_{k}\|=\|g_{k}\|\leq\|\epsilon\|$,
the conditions in Proposition 3.1 can not hold. Therefore, the direction
satisfying $B_{k}p_{k}=-g_{k}$ may be not a descent direction anymore.

In our quasi-Newton setting, the aim is not only to generate descent
directions, but also to generate search directions of high quality. Suppose
$r_{k}(p_{k})$ is a quality measure for the search direction $p_{k}$ at
iteration $k$. The aim is to minimize the quality measure $r_{k}(p_{k})$ such
that $B_{k}p_{k}=-g_{k},B_{k}\succ 0$. Therefore, for a given quality measure
$r_{k}(p_{k})$, we can generate a search direction with highest quality by
solving the following optimization problem:

$\begin{array}[]{ll}{\displaystyle\mathop{\operator@font
minimize}_{p_{k},B_{k}}}&{r_{k}(p_{k})}\\\ \operator@font subject\
to&B_{k}p_{k}=-g_{k},\\\ &B_{k}\succ 0.\end{array}$ (3.1)

In addition, we typically require $B_{k}$ to have some additional properties,
as problem (3.1) becomes highly complex otherwise.

The model (3.1) is actually a deterministic model, where the noisy gradients
are deterministic. However, as mentioned in Section 1, in some practical
problems, the gradients themselves are random because of the randomness in the
original objective quadratic function. Therefore, it is more natural to view
the gradients as random vectors in model (3.1). At iteration $k$, it is
assumed that $\tilde{g}_{k}=\bar{g}_{k}+\tilde{\epsilon}$, where
$\tilde{\epsilon}$ is a random noise with mean equal to $0$. Hence,
considering the randomness and to overcome the shortcoming of model (3.1) as
mentioned in Remark 1, the model (3.1) can be formulated as the following
chance constrained model:

$\begin{array}[]{lll}{\displaystyle\mathop{\operator@font
minimize}_{t,p_{k},B_{k}}}&t\\\ \operator@font subject\
to&\mathbb{P}\left\\{r_{k}(p_{k})\leq t,B_{k}p_{k}=-\tilde{g}_{k},B_{k}\succ
0\right\\}\geq 1-\beta,\end{array}$ (3.2)

where $\beta\in(0,1)$ is a given probability level. The value of $\beta$
indicates the risk-aversion of the decision maker, where 0 is the most
conservative approach as we need to comply with the supremum value of the
underlying random vector, while higher values would make our solution averse
to risk. Even small values of $\beta$ can have significant impact to the
results [BHdMM+16], therefore studying the behaviour of the solution while
$\beta$ is 0 and close to 0 (typically 0.01 or 0.05) is the usual approach.
The chance constraint in problem (3.2) not only guarantees the validation of
quasi-Newton setting with probability at least $1-\beta$, but also controls
the quality measure of the search direction.

To show that the search direction obtained by solving problem (3.2) is a
descent direction with a high probability, we have the following proposition.

###### Proposition 3.2

Denote $\tilde{g}_{k}=\bar{g}_{k}+\tilde{\epsilon}$, where $\bar{g}_{k}$ is
the true gradient and $\tilde{\epsilon}$ is a random noise with mean equal to
$0$. If $B_{k}$ is invertible and
$\mathbb{P}\left\\{\|\tilde{\epsilon}\|<\frac{1}{\|B_{k}^{-1}\tilde{g}_{k}\|}\tilde{g}_{k}^{T}B_{k}^{-1}\tilde{g}_{k}\right\\}\geq
1-\alpha$, the direction $p_{k}$ satisfying the chance constraint in problem
(3.2) is a descent direction at point $x_{k}$ with probability at least
$\max\\{1-\alpha-\beta,0\\}$.

Proof. From Proposition 3.1, we have that the direction $p_{k}$ is a descent
direction if the following constraints are satisfied:

$B_{k}p_{k}=-\tilde{g}_{k},B_{k}\succ 0.$

Then, we have

$\displaystyle\mathbb{P}\left\\{\begin{array}[]{l}B_{k}p_{k}=-\tilde{g}_{k},B_{k}\succ
0,\\\
\|\epsilon\|<\frac{1}{\|B_{k}^{-1}\tilde{g}_{k}\|}\tilde{g}_{k}^{T}B_{k}^{-1}\tilde{g}_{k}\end{array}\right\\}$
$\displaystyle\geq$
$\displaystyle\mathbb{P}\left\\{\begin{array}[]{l}r_{k}(p_{k})\leq
t,B_{k}p_{k}=-\tilde{g}_{k},B_{k}\succ 0,\\\
\|\epsilon\|<\frac{1}{\|B_{k}^{-1}\tilde{g}_{k}\|}\tilde{g}_{k}^{T}B_{k}^{-1}\tilde{g}_{k}\end{array}\right\\}$
$\displaystyle\geq$
$\displaystyle\max\Bigg{\\{}\mathbb{P}\left\\{r_{k}(p_{k})\leq
t,B_{k}p_{k}=-\tilde{g}_{k},B_{k}\succ 0\right\\}$
$\displaystyle+\mathbb{P}\left\\{\|\epsilon\|<\frac{1}{\|B_{k}^{-1}\tilde{g}_{k}\|}\tilde{g}_{k}^{T}B_{k}^{-1}\tilde{g}_{k}\right\\}-1,0\Bigg{\\}}$
$\displaystyle\geq$ $\displaystyle\max\\{1-\alpha-\beta,0\\}.$

The second inequality comes from Fréchet inequality. Therefore, the conclusion
can be obtained.

In contrast to the deterministic case, we can still maintain a descent
direction with positive probability even for $\|\tilde{g}_{k}\|<\epsilon$, as
observed in the following remark.

###### Remark 2

From Proposition 3.2, we can observe that even $x_{k}$ is close to the optimal
solution and the value of $\|\tilde{g}_{k}\|$ is small, there still exists a
constant $0<\bar{\alpha}<1$ such that
$\mathbb{P}\left\\{\|\epsilon\|<\frac{1}{\|B_{k}^{-1}\tilde{g}_{k}\|}\tilde{g}_{k}^{T}B_{k}^{-1}\tilde{g}_{k}\right\\}\geq
1-\bar{\alpha}$. This implies that problem (3.2) can always provide a descent
direction $p_{k}$ with some probability no matter how close $x_{k}$ is to the
optimal solution.

A chance-constrained model is expected to be computationally expensive. Our
interest is to study its behavior, in particular to see how such an approach
might be able to achieve improved accuracy, as indicated by Remark 2.

###  3.1 Suggestion on quality measure

With these propositions concerning the search direction, an interesting issue
is how to determine the quality measure of the search direction $p_{k}$ when
the method is applied to an unconstrained quadratic optimization problem. We
suggest an approach based on the characterization given in (1.2). In the exact
arithmetic case, (1.2) gives $(g_{i+1}-g_{i})^{T}p_{k}=0$, $i=0,\dots,k-1$.
Therefore, in this situation, the desired $p_{k}$ would give a global minimum
zero with respect to the measure

$r_{k}(p_{k})=\sum_{i=1}^{k-1}\left((g_{i+1}-g_{i})^{T}p_{k}\right)^{2}.$
(3.5)

In the exact arithmetic case, $p_{k}$ is orthogonal to the affine span of the
generated gradients. In the noisy setting, we will use this measure to show
how close the direction is to the characterization in (1.2).

Motivated by the discussion above, at iteration $k$, for a given positive
semidefinite matrix $V_{k}$, we will consider quasi-Newton matrices in the
family given by (2.5) and wish to find $\rho$ as the solution of

$\begin{array}[]{ll}{\displaystyle\mathop{\operator@font
minimize}_{\rho,p_{k},B_{k}\succ
0}}&\displaystyle{\sum_{i=0}^{k-1}((g_{i+1}-g_{i})^{T}p_{k})^{2}}\\\
\operator@font subject\ to&B_{k}p_{k}=-g_{k},\\\
&B_{k}=V_{k}+\sum_{i=0}^{k-1}\rho_{i}(g_{i+1}-g_{i})(g_{i+1}-g_{i})^{T}.\end{array}$
(D)

For the given $V_{k}$, in exact arithmetic and under exact linesearch, the
specific values of $\rho$ have no impact on the search direction. However, the
specific values may ensure nonsingularity of $B_{k}$ and numerical stability.
Note also that the sum of rank-one matrices in (2.5) is similar to terms
present in the BFGS Hessian approximation of (2.2).

Considering the possible randomness in the gradient $\tilde{g}_{k}$, the
chance constrained model (3.2) should take equations (2.5) and (3.5) into
consideration. In addition, from (2.9), we can notice that the matrix $V_{k}$
can be dependent on the noisy gradient $\tilde{g}_{k}$ in some cases, which
implies that $V_{k}$ is random. Therefore, in the random case, we denote

$B_{k}=\tilde{V}_{k}+\sum_{i=0}^{k-2}\rho_{i}(g_{i+1}-g_{i})(g_{i+1}-g_{i})^{T}+\rho_{k-1}(\tilde{g}_{k}-g_{k-1})(\tilde{g}_{k}-g_{k-1})^{T},$
(3.6)

where $\tilde{V}_{k}$ is a random matrix. Assume that the gradients
$g_{i},i=0,\dots,k-2$, have been evaluated, which can be the realized values
of noisy gradients or average values of noisy gradient samples. Then, based on
model (D), the chance constrained model for finding the search direction,
associated with model (3.2), can be formulated as

$\begin{array}[]{lll}{\displaystyle\mathop{\operator@font
minimize}_{t,\rho,p_{k},B_{k}\succ
0}}&\displaystyle\sum_{i=0}^{k-1}t_{i}^{2}\\\ \operator@font subject\
to&\mathbb{P}\left\\{\begin{array}[]{l}-t_{i}\leq(g_{i+1}-g_{i})^{T}p_{k}\leq
t_{i},\ i=0,\dots,k-2,\\\ -t_{k-1}\leq(\tilde{g}_{k}-g_{k-1})^{T}p_{k}\leq
t_{k-1},\\\ B_{k}p_{k}=-\tilde{g}_{k},\\\
B_{k}=\tilde{V}_{k}+\sum_{i=0}^{k-2}\rho_{i}(g_{i+1}-g_{i})(g_{i+1}-g_{i})^{T}\\\
\phantom{B_{k}=}+\rho_{k-1}(\tilde{g}_{k}-g_{k-1})(\tilde{g}_{k}-g_{k-1})^{T}\end{array}\right\\}\geq
1-\beta.\end{array}$ (C)

In model (C), only $\tilde{g}_{k}$ is the random vector, while
$g_{i},i=0,\dots,k-2$, are constant values.

Propositions 3.1 and 3.2 guarantee that the search directions obtained by
solving model (D) and model (C) are descent directions under certain
conditions, respectively.

##  4 Reformulation for the low-rank quasi-Newton setting

Models (D) and (C) are hard to solve given the non-convex nature of the
equality-constraints created by the condition $B_{k}p_{k}=-g_{k}$ and the
associated condition on $B_{k}$, and the nature of the chance constraints. In
the low-rank setting, however, the constraints related to the quasi-Newton
matrix can be simplified significantly.

In our low-rank setting, the only unknown parameter is for the memoryless BFGS
matrix, where we may write

$B_{k}=V_{k}+\rho(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T},$

and treat $\rho$ as a variable. For this case, we may allow further
simplification by circumventing the possible singularity of $V_{k}$ by letting
$\hat{\rho}_{k-1}$ be the value given by the secant condition (2.8), and
writing

$B_{k}=\hat{B}_{k}+(\rho-\hat{\rho}_{k-1})(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T},$
(4.1)

for

$\hat{B}_{k}=V_{k}+\hat{\rho}_{k-1}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T}.$ (4.2)

For the memoryless BFGS matrix, it holds that $V_{k}$ is positive semidefinite
with at most one zero eigenvalue in addition to $V_{k}(g_{k}-g_{k-1})\neq 0$,
so that $\hat{B}_{k}\succ 0$ as $\hat{\rho}_{k-1}>0$. The point of introducing
$\hat{\rho}_{k-1}$ and $\hat{B}_{k}$ is to give a nonsingular and positive
definite $\hat{B}_{k}$, which may be used as a foundation for optimizing over
$\rho$. We may therefore view the optimization over $\rho$ as the potential
for improving over utilizing the secant condition.

For a $\hat{\rho}_{k-1}$ such that $\hat{B}_{k}\succ 0$ in (4.1), the Sherman-
Morrison formula gives

$B_{k}^{-1}=\hat{B}_{k}^{-1}+\gamma\hat{B}_{k}^{-1}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T}\hat{B}_{k}^{-1}$
(4.3)

for

$\gamma=-\frac{(\rho-\hat{\rho}_{k-1})}{1+(\rho-\hat{\rho}_{k-1})\left(g_{k}-g_{k-1}\right)^{T}\hat{B}_{k}^{-1}\left(g_{k}-g_{k-1}\right)},$
(4.4)

so that an explicit expression for $p_{k}$ may be given as

$p_{k}=-B_{k}^{-1}g_{k}=-\hat{B}_{k}^{-1}g_{k}-\gamma(g_{k}-g_{k-1})^{T}\hat{B}_{k}^{-1}g_{k}\hat{B}_{k}^{-1}(g_{k}-g_{k-1}).$

Note that there is a one-to-one correspondence between $\gamma$ and $\rho$ as
(4.4) gives

$\rho-\hat{\rho}_{k-1}=-\frac{\gamma}{1+\gamma\left(g_{k}-g_{k-1}\right)^{T}\hat{B}_{k}^{-1}\left(g_{k}-g_{k-1}\right)}.$
(4.5)

In addition, if $\hat{B}_{k}\succ 0$, then (4.4) and (4.5) show that
$B_{k}\succ 0$ if and only if the equivalent conditions

$\gamma>-\frac{1}{(g_{k}-g_{k-1})^{T}\hat{B}_{k}^{-1}(g_{k}-g_{k-1})}\text{and}\rho-\bar{\rho}>-\frac{1}{(g_{k}-g_{k-1})^{T}\hat{B}_{k}^{-1}(g_{k}-g_{k-1})}$

hold. This is a consequence of these lower bounds defining an interval around
$\hat{B}_{k}$ and $\hat{B}_{k}^{-1}$ respectively, where $B_{k}$ and
$B_{k}^{-1}$ are well defined.

Summarizing, we may formulate the simplified problem as

$\begin{array}[]{ll}{\displaystyle\mathop{\operator@font
minimize}_{\gamma,p_{k}}}&\displaystyle\sum_{i=0}^{k-1}((g_{i+1}-g_{i})^{T}p_{k})^{2}\\\
\operator@font subject\
to&p_{k}=-\hat{B}_{k}^{-1}g_{k}-\gamma(g_{k}-g_{k-1})^{T}\hat{B}_{k}^{-1}g_{k}\hat{B}_{k}^{-1}(g_{k}-g_{k-1}),\\\\[2.84526pt]
&\gamma>-\frac{1}{\left(g_{k}-g_{k-1}\right)^{T}\hat{B}_{k}^{-1}\left(g_{k}-g_{k-1}\right)},\end{array}$
(DS)

which is a convex constrained quadratic program if a tolerance is introduced
for the strict lower bound on $\gamma$. For this problem, we may eliminate
$p_{k}$ to get one variable only, $\gamma$. Note the one-to-one correspondence
given by (4.5) which allows us to recover $\rho$ from $\gamma$.

For the chance-constrained model (C), analogous to (3.6), we denote

$B_{k}=\tilde{B}_{k}+(\rho-\hat{\rho}_{k-1})(\tilde{g}_{k}-g_{k-1})(\tilde{g}_{k}-g_{k-1})^{T},$
(4.6)

for

$\tilde{B}_{k}=V_{k}+\hat{\rho}_{k-1}(\tilde{g}_{k}-g_{k-1})(\tilde{g}_{k}-g_{k-1})^{T},$
(4.7)

which is random due to the randomness of $\tilde{g}_{k}$. $V_{k}$ is
deterministic as defined in (2.6). Then, analogous simplification and
reformulation of (C) gives

$\begin{array}[]{llll}{\displaystyle\mathop{\operator@font
minimize}_{t,\gamma,p_{k}}}&\displaystyle\sum_{i=0}^{k-1}t_{i}^{2}\\\
\operator@font subject\
to&\displaystyle\mathbb{P}\left\\{\begin{array}[]{l}-t_{i}\leq(g_{i+1}-g_{i})^{T}p_{k}\leq
t_{i},\ i=0,\dots,k-2,\\\\[2.84526pt]
-t_{k-1}\leq(\tilde{g}_{k}-g_{k-1})^{T}p_{k}\leq t_{k-1},\\\
p_{k}=-\tilde{B}_{k}^{-1}\tilde{g}_{k}\\\
\phantom{p_{k}=}-\gamma(\tilde{g}_{k}-g_{k-1})^{T}\tilde{B}_{k}^{-1}\tilde{g}_{k}\tilde{B}_{k}^{-1}(\tilde{g}_{k}-g_{k-1}),\\\\[2.84526pt]
\displaystyle\gamma>-\frac{1}{\left(\tilde{g}_{k}-g_{k-1}\right)^{T}\tilde{B}_{k}^{-1}\left(\tilde{g}_{k}-g_{k-1}\right)}\end{array}\right\\}\geq
1-\beta,\\\ \end{array}$ (CS)

Our proposed model (CS) can be read as follows: the obtained solution
$\gamma$, which will be transformed to $\rho$ by (4.5), will have a
probability of at least $1-\beta$ to obtain a direction $p_{k}$ from
$\tilde{g}_{k}$ that will be a descent direction. This means our approach
obtains a value of $\gamma$ robust enough to point us in a descent direction
$p_{k}$ with a probability bounded given the probabilistic nature of the
gradient.

Chance constrained models are often non-convex in general and hard to solve
[SDR14]. However, different equivalent formulations can be applied to obtain
an analytical solution or approximate the chance constraints. In general, it
is often difficult to get an analytical solution, since it always requires
strict assumptions on the probability distribution of random variables and the
structure of chance constraints, which making the situation specific and not
general enough. In constrast to the analytical solution, sample average
approximation and scenario approximation are two general approaches without
much assumptions on the random variables, which will be applied to solve the
chance constrained model in the following sections.

###  4.1 Deterministic equivalent formulation

Model (CS) can not be solved directly in its current state, as it is not a
deterministic problem. Therefore, a sample average approximation (SAA)
approach is proposed, given its flexibility to work under any type of
stochastic variables. The first step is to formulate it as a deterministic
problem that approximates the solution of (CS) and that can be solved by some
solvers. Let $\Omega$ be the set of sample, $|\Omega|=S$,
$\tilde{g}_{k\omega},\omega\in\Omega$ be the i.i.d. noisy gradient samples,
$\delta>0$ a sufficiently small real number and $0\leq K<k$ be an integer
number that represent the time window to be considered in the model. All
samples have a probability of $1/S$. Then, a deterministic equivalent
formulation of (CS) using sample average approximation (SAA) is as follows:

$\begin{array}[]{lll}{\displaystyle\mathop{\operator@font
minimize}}&\sum_{i=\max(0,k-K)}^{k-1}t_{i}^{2}\\\ \operator@font subject\
to&(g_{i+1}-g_{i})^{T}p_{k\omega}\geq-t_{i}-Mz_{\omega},\ i\in I_{K},\
\omega\in\Omega,\\\ &(g_{i+1}-g_{i})^{T}p_{k\omega}\leq t_{i}+Mz_{\omega},\
i\in I_{K},\ \omega\in\Omega,\\\
&-t_{k-1}-Mz_{\omega}\leq(\tilde{g}_{k\omega}-g_{k-1})^{T}p_{k\omega}\leq
t_{k-1}+Mz_{\omega},\ \omega\in\Omega,\\\
&p_{k\omega}=-\tilde{B}_{k\omega}^{-1}\tilde{g}_{k\omega}\\\
&\phantom{p_{k\omega}=}-\gamma(\tilde{g}_{k\omega}-g_{k-1})^{T}\tilde{B}_{k\omega}^{-1}\tilde{g}_{k\omega}\tilde{B}_{k\omega}^{-1}(\tilde{g}_{k\omega}-g_{k-1}),\
\omega\in\Omega,\\\\[2.84526pt]
&\displaystyle\gamma+Mz_{\omega}\geq-\frac{1}{\left(\tilde{g}_{k\omega}-g_{k-1}\right)^{T}\tilde{B}_{k\omega}^{-1}\left(\tilde{g}_{k\omega}-g_{k-1}\right)}+\delta,\
\omega\in\Omega,\\\ &\sum_{\omega\in\Omega}z_{\omega}\leq\lfloor
S\beta\rfloor,\\\ &z_{\omega}\in\\{0,1\\},\ \omega\in\Omega,\end{array}$ (CSA)

where $I_{K}=\\{\max(0,k-K),\dots,k-2\\}$ and $K\leq k$ indicates the number
of gradients to be considered in the model. Once we obtain the optimal value
of $\gamma$, then $\rho_{k}$ is approximated by

$\rho_{k-1}=\frac{\gamma}{1+\gamma\left(\bar{g}_{k}-g_{k-1}\right)^{T}\bar{B}_{k}^{-1}\left(\bar{g}_{k}-g_{k-1}\right)},$
(4.8)

where $\bar{g}_{k}=\frac{1}{S}\sum_{\omega\in\Omega}g_{k\omega}$ and
$\bar{B}_{k}$ is obtained by replacing $\tilde{g}_{k}$ in (4.7) with
$\bar{g}_{k}$. Finally, $p_{k}$ is obtained using the equation

$p_{k}=-B_{k}^{-1}\bar{g}_{k}=-\bar{B}_{k}^{-1}\bar{g}_{k}-\gamma(\bar{g}_{k}-g_{k-1})^{T}\bar{B}_{k}^{-1}\bar{g}_{k}\bar{B}_{k}^{-1}(\bar{g}_{k}-g_{k-1}),$
(4.9)

i.e., using the average of the gradients at iteration $k$.

It can be noticed that if $\beta>0$ then (CSA) is a mixed integer program
whose complexity will be tied to the number of dimensions and samples used to
solve the problem. Since at every iteration new gradients are added, the
dimensionality of the problem will grow at each step regardless. And to
guarantee the quality of solution by SAA, the sample size should not be too
small if the dimension is large.

Therefore the complexity of this approach grows at each step, leading to
increasing solving times which will be an issue on long runs with low
convergence speed. However, by using the value of parameter $K$ with values
greater than 0, we can use a limited memory or memoryless, implementation to
avoid this.

If $\beta=0$, then all binary variables must be set to zero and can be
eliminated from the problem, creating a continuous linear program which is
much simpler to solve. This is commonly referred to the scenario approach,
where all possible sampled scenarios of the random variables are being
considered. This also implies that the solution will be closely tied to the
most conservative of the sampled scenarios. There are other approaches to
simplify and approximate this model present in the literature [AX18].

##  5 Computational results for the quasi-Newton methods

Two sets of results are presented: First, we consider a set of randomly
generated problems intended to illustrate the properties and methodologies
proposed in this paper. The second set of problems are real-life instances
from the CUTEst test set [GOT15], to test the applicability of these methods
in a more realistic environment.

A comparison of the results is provided using different models and/or
approximation formulations, and we discuss the practical implications obtained
with each method. All models are implemented in Python 3.7.10, using Gurobi
9.1 as a solver for the resulting optimization problems and all computation
were done on an Intel(R) i7 @ 2.7 GHz and 16 GB of memory over macOS 10.

The methods chosen in our experiments are as follows:

  * •

Steepest Descent (SD):

$B_{k}=I.$

  * •

Conjugate Gradients (CG): The symmetric CG method as presented in (2.9),

$B_{k}=\left(I-\frac{1}{g_{k-1}^{T}p_{k-1}}g_{k}p_{k-1}^{T}\right)\left(I-\frac{1}{g_{k-1}^{T}p_{k-1}}p_{k-1}g_{k}^{T}\right).$

  * •

BFGS, as presented in (2.2) and where $\rho_{k-1}$ is given by the secant
condition in (2.8),

$B_{k}=B_{k-1}+\frac{1}{g_{k-1}^{T}p_{k-1}}g_{k-1}g_{k-1}^{T}+\frac{1}{\alpha_{k-1}(g_{k}-g_{k-1})^{T}p_{k-1}}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T}.$

  * •

Memoryless-BFGS (ml-BFGS), as presented in (2.7) and where $\rho_{k-1}$ is
given by the secant condition (2.8),

$B_{k}=I-\frac{1}{p_{k-1}^{T}p_{k-1}}p_{k-1}p_{k-1}^{T}+\rho_{k-1}(g_{k}-g_{k-1})(g_{k}-g_{k-1})^{T}.$

  * •

Chance-Constrained Quasi Newton (CCQN $\beta$). The search direction is
obtained by first solving (CSA) for $K=0$ and $\beta$, then obtaining
$\rho_{k-1}$ from (4.8) and finally $p_{k}$ from (4.9).

  * •

Limited Memory Chance-Constrained Quasi Newton (lm-CCQN $\beta$). Same as
CCQN, but $0<K<k$.

In order to standardize our results across different experiments, _performance
profiles_ are used in two different analyses: First, for the number of steps
required by the method to break a certain gradient norm value which will
define as the tolerance level (denoted as $tol$), then by taking the method
with the lowest possible value of steps as the comparison point, we show how
much larger the values of the other methods are compared to this minimum.
Next, a performance profile is created to detect the minimum value of the
gradient norm. Finally, the total amount of times the different experiments
using a set method were able to reach a set thresholds is summed and
presented: Starting from 1 (being the minimum) up to 20 (i.e. 20 times the
value of the minimum). This methodology is fully expanded in [EF21].

Since noise will severely distort the gradient norms once it reaches a certain
point in the run, a set of performance profiles are created for tolerance
values close to the noise variance, i.e. if the variance $\sigma^{2}=10^{-2}$
then these performance profiles will be set to tolerances such as
$10^{-1},10^{-2}$ and $10^{-3}$. Higher does not show significant differences
to the deterministic case, and lower can cause the method to not reach any
threshold.

For every experiment, we applied the following algorithm: At iteration $k$,
$tol$ denotes our tolerance threshold of precision for the norm of the
gradient of the solution obtained, $\bar{g}_{k}$ is the average value of the
gradient $\tilde{g}_{k}$, $K\leq k$ denotes the number of gradients to be used
in the calculations and $MaxK$ is the maximum number of steps. When
calculating the step length $\alpha$, an exact line search approach is used
with the value of the gradient without noise (the deterministic value of
$g_{k}$), as the objective is to isolate the effects of each method in finding
a descent direction. The algorithm goes as follows:

Algorithm 1 General solving algorithm.

$k\leftarrow 0$;

$x_{k}\leftarrow$ initial point;

$\tilde{g}_{k\omega}\leftarrow
Hx_{k}+c+\epsilon_{\omega},\quad\omega\in\Omega$;

$\bar{g}_{k}\leftarrow\frac{1}{S}\sum_{\omega\in\Omega}\tilde{g}_{k\omega}$;

while $\|g_{k}\|_{2}>tol\ \textbf{and}\ k\leq MaxK$ do

if method = CCQN then

$\gamma\leftarrow$ solution to (CSA) using $K$ gradients;

$\rho\leftarrow$ solution to (4.8);

end if

$B_{k}\leftarrow$ from method;

$p_{k}\leftarrow$ solution to $B_{k}p_{k}=-\bar{g}_{k}$;

$\alpha_{k}\leftarrow-\frac{g_{k}^{T}p_{k}}{p_{k}^{T}Hp_{k}}$;

$x_{k+1}\leftarrow x_{k}+\alpha_{k}p_{k}$;

$k\leftarrow k+1$;

$\tilde{g}_{k\omega}\leftarrow
Hx_{k}+c+\epsilon_{\omega},\quad\omega\in\Omega$;

$\bar{g}_{k}\leftarrow\frac{1}{S}\sum_{\omega\in\Omega}\tilde{g}_{k\omega}$;

end while

All experiments are repeated 30 times using different random number generator
seeds and using 20 samples of noisy gradients at each step, therefore the
performance profiles also separate each method and experiment by seed. We used
a maximum amount of steps $MaxK=500$ and $K=10$ for the lm-CCQN method.

###  5.1 Results for randomly created problems

The first experiment is a set of randomly generated unconstrained quadratic
problems. For each problem, the Hessian matrix $H$ is defined as
$H=Q^{T}Q+\epsilon\mathop{\operator@font diag}(U_{1,n})$, where
$Q=aJ_{1,n}+(b-a)U_{n,n}$, $a,b\in\mathbb{R}$, $J_{n,m}$ is the unit $n\times
m$ matrix, $U_{n,m}$ is a $n\times m$ matrix, each component of $U_{n,m}$ is
randomly generated following a Uniform(0,1) distribution and $\epsilon>0$ is a
sufficiently small number. The vector $c$ is randomly generated as
$c=U_{n,1}$. In our experiments, we defined $a=-1,b=1,n=100$ and
$\epsilon=0.3$.

(a) $\sigma^{2}=10^{-6}$

(b) $\sigma^{2}=10^{-2}$

Figure 1: Average log norm of the gradient at step $k$ for each tested method
with different noise variances.

Figure 1 shows the behaviour of each method in the two selected noise levels.
The traditional approaches, such as CG or BFGS, appear to converge faster but
the average log norm of the gradient can not surpass the tolerance threshold
$tol$ smaller than the noise variance $\sigma^{2}$, while for CCQN and lm-CCQN
(regardless of the value $\beta$) the average log norm of the gradient can
surpass this barrier, albeit slower compared to the results presented by ml-
BFGS.

In this experiment, the chosen value of $MaxK$ was not large enough for SD to
show convergence as seen in figure 1, however we ran the same experiment for
this method using a larger value, showing that average log norm of the
gradient found by SD can suprass the $tol\leq\sigma^{2}$ barrier, similar to
CCQN, lm-CCQN and ml-BFGS.

(a) $tol=10^{-5},\sigma^{2}=10^{-6}$

(b) $tol=10^{-1},\sigma^{2}=10^{-2}$

(c) $tol=10^{-6},\sigma^{2}=10^{-6}$

(d) $tol=10^{-2},\sigma^{2}=10^{-2}$

Figure 2: Performance profiles for different tolerances and noise variance
levels.

Figure 2 shows the performance profiles of each method for two noise variance
levels under two tolerance thresholds. In the traditional approaches, we can
observe that ml-BFGS performs better overall, while CG and BFGS perform well
under low noise variance but poorly under larger values. On the other hand,
the performance of CCQN and lm-CCQN methods do not show significant
differences for different tolerance and noise variance levels. Furthermore,
the difference between using $\beta>0$ or 0 is not significant, which means we
implement a convex approximation using the scenario approach ($\beta=0$),
avoiding the usage of a mixed integer linear program which becomes difficult
to solve for larger problems.

(a) $\sigma^{2}=10^{-6}$

(b) $\sigma^{2}=10^{-2}$

Figure 3: Performance profile of the minimum gradient norm for different noise
variance levels.

Figure 3 shows the performance profiles of the minimum gradient norm found in
the set of problems with different seeds for two noise variance levels. When
$\sigma^{2}=10^{-6}$, we can observe that ml-BFGS is able to obtain the
minimum value consistently, i.e. for every problem and every seed, which is
followed by CCQN and lm-CCQN. When $\sigma^{2}=10^{-2}$, lm-CCQN performs
better than ml-BFGS, and followed by CCQN. The classical methods BFGS, CG and
SD perform worse than the previously discussed methods independent of the
noise variance, and specifically the minimum gradient norms found by SD are
larger than any of minimum norms found by all the other methods.

###  5.2 Results for CUTEst problems

In our experiments, we will compare the performance of the same approaches
presented in the last section applied to different problems from the CUTEst
test set [GOT15], specifically quadratic, unconstrained and number of
variables chosen by the user (QUV using CUTEst classification system).
However, only 6 problems families fall in this category, therefore we
implemented a second batch of problems by adding those unconstrained sum of
squares problems (SUV using CUTEst classification system) which had a positive
definite Hessian at the starting point and left it constant throughout the
solving scheme. This brought the total amount of problems to 22\. In our
numerical experiments, results for noise variance $10^{-2}$ and $10^{-6}$ did
not show significant differences. Therefore, in this section we only focus on
results with a noise variance of $10^{-2}$.

(a) $tol=10^{0}$

(b) $tol=10^{-1}$

(c) $tol=10^{-2}$

(d) $tol=10^{-3}$

Figure 4: Performance profiles for different tolerance levels of the CUTEst
problems.

Figure 4 shows the performance profiles of different methods solving CUTEst
problems under different tolerance levels. We can observe that the performance
of BFGS and CG can be good and SD performs worst among all the tested methods,
when the tolerance level is larger than the noise variance. However, as shown
in Figure 4 (d), SD can perform better than ml-BFGS, BFGS and CG. ml-BFGS can
show its efficiency and effectiveness in most cases. CCQN and lm-CCQN present
their robustness under different tolerance levels, and perform better when the
tolerance level is close or smaller than the noise variance. When
$tol=10^{-3}$, CCQN and lm-CCQN perform best, with the latter being slightly
ahead.

When the tolerance level becomes much smaller, the performance profile of each
method will depend on the problem, except CCQN and lm-CCQN. The ml-BFGS method
will not always be the best, while the performances of CCQN and lm-CCQN are
among the top three methods, which shows their robustness if the noise
variance is larger the tolerance level. This behaviour can be observed in
Appendix A for each individual problem.

Figure 5: Performance profile of the minimum gradient norm found of the CUTEst
problems.

Figure 5 the performance profiles of the minimum gradient norm found in the
CUTEst set of problems with different seeds. CCQN is able to obtain the
minimum value consistently, followed by the lm-CCQN, SD and ml-BFGS (in that
order), while BFGS and CG perform worst than all the other studied methods.
The chance-constrained methods were able to obtain the minimum value faster
and more frequently, even if compared it to ml-BFGS which had a good
performance across our previously shown profiles.

##  6 Conclusion

In this paper, we have studied low-rank quasi-Newton methods for minimizing a
strictly convex quadratic function in a noisy framework. We have considered a
memoryless BFGS method and compared to a BFGS method, the method of conjugate
graddients and steepest descent. In order to potentially improve the
performance of the low-rank quasi-Newton method, a chance constrained
stochastic optimization model has also been formulated. The secant condition
is here replaced by solving a one-dimensional convex quadratic programming
problem. The proposed chance constrained model, which can be solved
effectively by sample average approximation method or scenario approach, has
been proven to provide a descent search direction with a high probability in
the random noisy framework, while the deterministic model may fail to provide
a descent direction, if the noise level is large.

In the numerical experiments, we compare classical methods and the proposed
chance constrained model in a noisy setting. Results of ml-BFGS and CCQN show
promise when solving problems with uncertainty in the gradient, however the
latter is more consistent and its performance appear to be independent of the
problem, while the former does not. The performance of chance-constrained
model (and its different iterations) appears to be in the top three in terms
of convergence speed under different tolerance thresholds. Furthermore, while
studying the behaviour of all the models, the minimal value of gradient norm
was consistently found by the approach based on chance constrained model.
Therefore, we believe that the usage of more advanced solving algorithms than
the one presented (i.e. stochastic inexact linesearch) could further improve
the results presented in this paper.

Finally, our intention is to investigate the behavior and the interplay
between quality and robustness of the low-rank quasi-Newton method, especially
in the case of large noise and multiple copies of gradients. Both the
theoretical and numerical results show that we can gain the robustness and
accuracy of the computed solution with the chance constrained model, although
the computational cost can be high. This shows the potential to be further
considered and explored in convex optimization problems.

## References

  * [AX18] S. Ahmed and W. Xie.  Relaxations and approximations of chance constraints under finite distributions.  Mathematical Programming, 170(1), 43–65, 2018.
  * [BHdMM+16] J. Barrera, T. Homem-de Mello, E. Moreno, B. K. Pagnoncelli, and G. Canessa.  Chance-constrained problems and rare events: an importance sampling approach.  Mathematical Programming, 157(1), 153–189, 2016.
  * [BHNS16] R. H. Byrd, S. L. Hansen, J. Nocedal, and Y. Singer.  A stochastic quasi-Newton method for large-scale optimization.  SIAM Journal on Optimization, 26(2), 1008–1031, 2016\.
  * [BNT16] A. S. Berahas, J. Nocedal, and M. Takac.  A multi-batch L-BFGS method for machine learning.  In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems 29 (NIPS 2016), volume 29, 2016.
  * [DW84] J. E. Dennis and H. F. Walker.  Inaccuracy in quasi-Newton methods - local improvement theorems.  Mathematical Programming Study, 22(DEC), 70–85, 1984\.
  * [EF21] D. Ek and A. Forsgren.  Exact linesearch limited-memory quasi-Newton methods for minimizing a quadratic function.  Computational Optimization and Applications, 79(3), 789–816, 2021.
  * [FO18] A. Forsgren and T. Odland.  On exact linesearch quasi-Newton methods for minimizing a quadratic function.  Computational Optimization and Applications, 69(1), 225–241, 2018\.
  * [GOT15] N. I. M. Gould, D. Orban, and P. L. Toint.  CUTEst: a constrained and unconstrained testing environment with safe threads for mathematical optimization.  Computational Optimization and Applications, 60(3), 545–557, 2015\.
  * [LMH15] A. Lucchi, B. McWilliams, and T. Hofmann.  A variance reduced stochastic Newton method.  Preprint arXiv:1503.08316 [cs.LG], 2015.
  * [MNJ16] P. Moritz, R. Nishihara, and M. Jordan.  A linearly-convergent stochastic L-BFGS algorithm.  In Artificial Intelligence and Statistics, pages 249–258. PMLR, 2016.
  * [MR13] A. Mokhtari and A. Ribeiro.  Regularized stochastic BFGS algorithm.  In 2013 IEEE Global Conference on Signal and Information Processing, pages 1109–1112. IEEE, 2013.
  * [MR15] A. Mokhtari and A. Ribeiro.  Global convergence of online limited memory BFGS.  The Journal of Machine Learning Research, 16(1), 3151–3181, 2015\.
  * [Saa03] Y. Saad.  Iterative methods for sparse linear systems.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2003\.
  * [SDR14] A. Shapiro, D. Dentcheva, and A. Ruszczyński.  Lectures on stochastic programming: modeling and theory.  SIAM, Philadelphia, 2014.
  * [SXBN20] H.-J. M. Shi, Y. Xie, R. Byrd, and J. Nocedal.  A noise-tolerant quasi-Newton algorithm for unconstrained optimization.  Preprint arXiv:2010.04352 [math.OC], 2020.
  * [XBN20] Y. Xie, R. H. Byrd, and J. Nocedal.  Analysis of the BFGS method with errors.  SIAM Journal on Optimization, 30(1), 182–209, 2020.

##  Appendix A Average log norms gradients of the CUTEst problems

In this section, the average log norm of each CUTEst problem is presented. The
objective is to make evident the difference between problems as stated in the
experiment section 5.2. All of these results are presented for
$\sigma^{2}=10^{-2}$.

(a) ARGLINA

(b) BDQRTIC

Figure 6: Average log norm of the gradient at step $k$ for each tested method.

(a) CHAINWOO

(b) EXTROSNB

Figure 7: Average log norm of the gradient at step $k$ for each tested method.

(a) LIARWHD

(b) PENALTY1

Figure 8: Average log norm of the gradient at step $k$ for each tested method.

(a) PENALTY2

(b) SROSENBR

Figure 9: Average log norm of the gradient at step $k$ for each tested method.

(a) TQUARTIC

(b) WOODS

Figure 10: Average log norm of the gradient at step $k$ for each tested
method.

(a) INTEQNELS

(b) CHNROSNB

Figure 11: Average log norm of the gradient at step $k$ for each tested
method.

(a) ERRINROS

(b) ERRINRSM

Figure 12: Average log norm of the gradient at step $k$ for each tested
method.

(a) DIXON3DQ_n1000

(b) DIXON3DQ_n100

Figure 13: Average log norm of the gradient at step $k$ for each tested
method.

(a) DIXON3DQ

(b) DQDRTIC_n100

Figure 14: Average log norm of the gradient at step $k$ for each tested
method.

(a) DIXON3DQ

(b) DQDRTIC_n100

Figure 15: Average log norm of the gradient at step $k$ for each tested
method.

(a) TRIDIA_n100

(b) TESTQUAD_n1000

Figure 16: Average log norm of the gradient at step $k$ for each tested
method.

(a) HILBERTB_n50

(b) TOINTQOR_n50

Figure 17: Average log norm of the gradient at step $k$ for each tested
method.