# Soft quasi-Newton:
Guaranteed positive definiteness by relaxing the secant constraint

Erik Berglund, Jiaojiao Zhang and Mikael Johansson CONTACT Erik Berglund.
Email<EMAIL_ADDRESS>School of Electrical Engineering and Computer Science,
KTH Royal Institute of Technology, Sweden

###### Abstract

We propose a novel algorithm, termed soft quasi-Newton (soft QN), for
optimization in the presence of bounded noise. Traditional quasi-Newton
algorithms are vulnerable to such perturbations. To develop a more robust
quasi-Newton method, we replace the secant condition in the matrix
optimization problem for the Hessian update with a penalty term in its
objective and derive a closed-form update formula. A key feature of our
approach is its ability to maintain positive definiteness of the Hessian
inverse approximation. Furthermore, we establish the following properties of
soft QN: it recovers the BFGS method under specific limits, it treats positive
and negative curvature equally, and it is scale invariant. Collectively, these
features enhance the efficacy of soft QN in noisy environments. For strongly
convex objective functions and Hessian approximations obtained using soft QN,
we develop an algorithm that exhibits linear convergence toward a neighborhood
of the optimal solution, even if gradient and function evaluations are subject
to bounded perturbations. Through numerical experiments, we demonstrate
superior performance of soft QN compared to state-of-the-art methods in
various scenarios.

###### keywords:

quasi-Newton methods; general bounded noise; secant condition; penalty

##  1 Introduction

Although the importance of uncertainty and stochasticity has been recognized
since the early days of numerical optimization (see, _e.g._ , [29, 11]), it is
fair to say that the majority of solution algorithms have been developed under
the assumption that the objective and constraint functions (and often also
their gradients and Jacobians) can be evaluated accurately. Even stochastic
programs [4], that explicitly model the uncertainty in the underlying decision
problem, are often approached via their deterministic equivalents. Techniques
for analyzing the effect of noise on iterative algorithms have a rich history
(see, _e.g._ [28]) but it is only over the last few decades that the research
devoted to algorithms that operate reliably using noisy function and gradient
evaluations has virtually exploded. This increased interest has been fueled by
an increasing number of applications in machine learning and simulation
optimization.

Machine learning, a subset of artificial intelligence, automates pattern
discovery within data sets and leverages these patterns to make inferences
based on input data [17]. Optimization plays a foundational role in various
machine learning approaches as they aim to minimize the discrepancy between
model predictions and observed data, necessitating the solution of
optimization problems [31]. A particularly challenging setting in machine
learning arises when dealing with large datasets that cannot be processed
simultaneously. In such cases, practitioners employ data sampling techniques
to construct approximations of the objective function to be minimized, rather
than evaluating it exhaustively [37, 6]. This situation leads to the
formulation of stochastic optimization problems. In addition to settings where
data is sampled, stochastic optimization problems can arise from numerical
simulations that contain stochastic components. Simulation optimization, i.e.
optimization of functions that can only be evaluated through a stochastic
simulation, is an important area with applications in diverse domains such as
operations research, manufacturing, medicine, and engineering [2]. In these
simulations, the statistical properties of the noise are not as clear, and it
is natural to consider the errors as bounded noise with otherwise unknown
properties.

First-order optimization methods are popular for solving stochastic problems.
In particular, the stochastic gradient descent method [29] and variations
thereof have become the workhorse of modern large-scale machine learning [7].
The scalability and ease of implementation of these methods contribute to
their popularity. In addition, they can achieve the optimal oracle complexity
up to a multiplicative constant in many stochastic settings [1]. However,
first-order methods struggle on ill-conditioned and highly nonlinear problems.
For example, the iteration complexity of stochastic gradient descent on
Lipschitz smooth, strongly convex problems grows proportionally to the
condition number [7]. In the non-convex setting, empirical evidence suggests
that second-order methods have advantages when it comes to the robustness of
hyper-parameter tuning and the ability to escape saddle points [35] [12] [27].

While the use of Hessian information can significantly improve the iteration
complexity of optimization algorithms, it can often be expensive to compute.
To gain some of the benefits of Newton’s method in settings where only
gradient information is available, quasi-Newton (QN) methods build a model of
the Hessian using only first-order information. In the deterministic setting,
quasi-Newton methods have proven to be much more effective than gradient
descent for small- and medium-scale problems, both in theory and practice
[25]. For large-scale problems, limited memory quasi-Newton methods, and
L-BFGS [20] in particular, have been prominent. However, traditional quasi-
Newton methods are sensitive to the presence of random noise in the gradients.
Thus, some adaptions must be made to retain their advantages in the stochastic
setting. Currently, this is a very active topic of research.

A significant effort has been devoted to adapting quasi-Newton methods to
problems where the stochastic noise comes from the sampling of gradients.
These techniques combine stochastic gradients with carefully constructed
Hessian estimations. For example, the paper [5] constructs diagonal or low-
rank Hessian approximations originating from the secant condition. The study
[8] decouples the stochastic gradient and curvature estimate calculation by
replacing the gradient difference that is typically used in quasi-Newton
methods with a Hessian-vector product. A different approach presented in [23]
introduces a regularized BFGS method that incorporates regularization
techniques to ensure positive definiteness and boundedness of Hessian
estimations for strongly convex loss functions. For non-convex problems, the
paper [32] employs damping techniques to preserve the positive definiteness
and boundedness of Hessian estimations. Another set of studies [22, 18]
investigate incremental quasi-Newton approaches, wherein a solitary Hessian
approximation associated with a chosen objective loss from the finite sum is
updated, while all other Hessian approximations remain unchanged.

Another extension of stochastic gradient descent (SGD) for constructing
descent directions is based on combining variance reduction techniques with
quasi-Newton algorithms to mitigate the noise impact introduced by stochastic
sampling [24, 36, 10].

In the setting of stochastic optimization with general bounded noise, quasi-
Newton methods have received considerably less attention. A study by Xie et
al. analyzed the behavior of BFGS in the presence of errors in the function
and gradient evaluations [34], and showed that with a lengthening technique to
modify the line-search procedure, BFGS can be made to converge linearly to a
neighborhood of the minimum of a strongly convex function. This work was
expanded upon in [30], where a modification of the line-search method was
proposed that did not require knowledge of the strong convexity constant.
Similar ideas were explored earlier in the derivative-free setting, where the
length of finite-difference interval was used to control the influence of
noise in quasi-Newton approximations [3]. Apart from adjusting the finite-
difference approximation to deal with bounded noise, another idea is to relax
the secant condition because, in the presence of noise, enforcing the secant
condition for quasi-Newton updates can be problematic. A recent work on
stochastic quasi-Newton methods takes this approach [16]. By replacing the
secant condition in the problem used to derive BFGS with a penalty term, the
authors derive a method which they call secant-penalized BFGS, or SP-BFGS for
short. They show that the resulting update formula can guarantee positive
definiteness, provided that the penalty parameter is small enough.

In this paper, we introduce a novel quasi-Newton algorithm for minimizing a
loss function in the presence of general bounded noise. Similarly to [16], we
relax the secant condition by replacing it with a penalty term, but we
consider a different matrix optimization problem and derive a new Hessian
approximation with several attractive features. We refer to the proposed
approach as soft QN, as it “softens” the constraint of having to satisfy the
secant condition. Unlike SP-BFGS, which requires its penalty parameter to be
sufficiently small to ensure positive definiteness of the Hessian
approximation, we prove that soft QN generates a positive definite matrix for
all positive penalty parameters, even when the loss functions are non-convex.
In addition, we establish several other attractive properties of soft QN,
including that it recovers BFGS under specific limits, handles positive and
negative curvature equally, and exhibits scale invariance. These
characteristics enhance the efficacy of the proposed soft QN in noisy
settings. For strongly convex and Lipschitz smooth loss functions, the
resulting quasi-Newton algorithm generates iterates that are guaranteed to
decrease the objective function at a linear rate to a neighborhood of the
optimal solution. Numerical experiments show that soft QN is competitive with
SP-BFGS and superior to other benchmark algorithms.

##  2 Preliminaries: quasi-Newton methods

In this section, we provide a brief review of the key ideas that underpin
traditional quasi-Newton methods [25]. For clarity, we focus on unconstrained
minimization problems

$\displaystyle\begin{array}[c]{ll}\underset{x\in\mathbb{R}^{n}}{\mbox{minimize}}&f(x),\end{array}$
(2)

where the objective function $f:\mathbb{R}^{n}\mapsto\mathbb{R}$ is continuous
and twice differentiable.

In essence, quasi-Newton methods are iterative optimization algorithms that
use approximations of the Hessian of the objective function to speed up their
convergence. These approximations are constructed based on the previous
iterates $(x_{k})$ generated by the method and the gradient of the objective
at these points. More specifically, given an initial point $x_{0}$, these
methods generate a sequence of iterates

$\displaystyle x_{k+1}$ $\displaystyle=x_{k}+\eta_{k}p_{k},$ (3)

where $\eta_{k}>0$ is a step-size and $p_{k}\in\mathbb{R}^{n}$ is a search
direction chosen to minimize a quadratic model $q_{k}$ of $f$ around $x_{k}$.
The model uses $f(x_{k})$, $\nabla f(x_{k})$ and an approximation $B_{k}$ of
$\nabla^{2}f(x_{k})$:

$\displaystyle q_{k}(p)$ $\displaystyle=f(x_{k})+\nabla
f(x_{k})^{\top}p+\frac{1}{2}p^{\top}B_{k}p.$

When $B_{k}$ is invertible, $q_{k}$ has the unique critical point

$p_{k}=-B_{k}^{-1}\nabla f(x_{k}).$ (4)

If $B_{k}$ is positive definite the model is minimized at this critical point
and $p_{k}$ is a descent direction for $f(x_{k})$, _i.e._ for a sufficiently
small step-size $\eta_{k}$, we have $f(x_{k}+\eta_{k}p_{k})<f(x_{k})$. Thus,
positive definiteness is a desirable property for $B_{k}.$

Different Hessian approximations $B_{k}$ lead to different quasi-Newton
methods. Most methods use Hessian approximations that satisfy the secant
condition

$y_{k}=B_{k+1}s_{k},$ (5)

where $s_{k}:=x_{k+1}-x_{k}$ and $y_{k}:=\nabla f(x_{k+1})-\nabla f(x_{k})$.
The condition is often motivated by the fact that it is satisfied by the
actual Hessian of the loss function. Specifically, since

$y_{k}=\int_{0}^{1}\nabla^{2}f(x_{k}+ts_{k})s_{k}dt,$

the mean value theorem of integrals says that there exists a $\theta\in[0,1]$
such that

$y_{k}=\nabla^{2}f(x_{k}+\theta s_{k})s_{k}.$

Requiring Equation (5) to be satisfied adds information about the true
objective function’s curvature to the model, but still leaves several degrees
of freedom for the choice of $B_{k+1}$. A key principle of quasi-Newton
methods is that they update $B_{k+1}$ from the previous Hessian approximation
$B_{k}$ by minimizing the size of the necessary modification [25]. This allows
the algorithm to accumulate approximate curvature information from previous
iterations and aids convergence when close to a local minimum, where the loss
function can be well-approximated by a quadratic one. Different metrics for
measuring closeness between $B_{k+1}$ and $B_{k}$ lead to different quasi-
Newton methods.

The most well-known quasi-Newton method, BFGS, has strong practical
performance and its limited-memory variant, L-BFGS, is used in various large-
scale optimization solvers [25]. To derive the BFGS update, one considers the
inverse of the Hessian approximation, $H_{k}$, and chooses $H_{k+1}$ to
satisfy the inverse secant condition

$H_{k+1}y_{k}=s_{k}.$ (6)

To ensure that $H_{k+1}$ is close to $H_{k}$, one considers the matrix
optimization problem

$\displaystyle\underset{H}{\text{minimize}}$
$\displaystyle\|W_{k}^{(1/2)}(H-H_{k})W_{k}^{(1/2)}\|_{F}^{2}$ (7) subject to
$\displaystyle Hy_{k}=s_{k},\;H=H^{T}.$

Here, $\|\cdot\|_{F}$ denotes the Frobenius norm, $W_{k}$ is any positive
definite matrix that satisfies $W_{k}s_{k}=y_{k}$, and $W_{k}^{(1/2)}$ is its
principal square root. Setting $H_{k+1}=H^{\star}$, where $H^{\star}$ is the
solution to (7), leads to the update formula

$\displaystyle H_{k+1}$
$\displaystyle=\left(I-\dfrac{s_{k}y_{k}^{T}}{y_{k}^{T}s_{k}}\right)H_{k}\left(I-\dfrac{y_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}\right)+\frac{s_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}.$
(8) For a derivation of this, see e.g. [15]. By the Sherman-Morrison-Woodbury
formula, the corresponding update for $B_{k}$ is $\displaystyle B_{k+1}$
$\displaystyle=B_{k}-\dfrac{B_{k}s_{k}s_{k}^{T}B_{k}}{s_{k}^{T}B_{k}s_{k}}+\dfrac{y_{k}y_{k}^{T}}{y_{k}^{T}s_{k}}.$
(9)

An important property of the BFGS update is that if $H_{k}$ is positive
definite, then $H_{k+1}$ will be positive definite if and only if
$s_{k}^{\top}y_{k}>0$ [25, p. 140]. In fact, a positive definite weight matrix
$W_{k}$ satisfying $W_{k}s_{k}=y_{k}$ exists if and only if
$s_{k}^{\top}y_{k}>0$, so the matrix optimization problem (7) is only defined
when this curvature condition holds.

The BFGS update also solves the matrix optimization problem

$\displaystyle\underset{B\succ 0}{\text{minimize}}$
$\displaystyle\operatorname{tr}(BB_{k}^{-1})-\log(\det(BB_{k}^{-1}))$ (10)
subject to $\displaystyle Bs_{k}=y_{k}.$

In this problem, the weighted Frobenius norm has been replaced by the log-det
divergence to measure the distance between $B$ and $B_{k}$. Since the
constraints on $B$ require that $B\succ 0$ and that $Bs_{k}=y_{k}$, it only
admits a solution if $s_{k}^{T}y_{k}>0$.

In the deterministic setting, BFGS is typically combined with a line-search
method to determine the step-size $\eta_{k}$ in Equation (3). If the line-
search technique enforces the Wolfe conditions [25, Equations (3.6a) and
(3.6b)], the algorithm ensures not only descent in each iteration but also
that $s_{k}^{T}y_{k}>0$ holds. In the stochastic setting, where $y_{k}$ and
$s_{k}$ are perturbed by noise, it may be impossible to satisfy this
condition. This contributes to making regular quasi-Newton methods very
sensitive to noise, to the extent that even a slight perturbation in the
gradients might cause them to break down.

As outlined in Section 1, multiple attempts have been made at adapting quasi-
Newton methods for stochastic problems. In the next section, we introduce a
method that modifies the matrix optimization problem that defines the BFGS
update to be more robust to noise, yet results in a simple update formula for
the Hessian inverse approximation.

##  3 The soft QN method

In this section, we introduce soft QN, a quasi-Newton method derived by
solving a matrix optimization problem with a “softened” secant condition. The
resulting algorithm is similar to the classical BFGS method but it is more
robust to noise and guaranteed to always yield a positive definite Hessian
inverse approximation. Below, we pose the matrix optimization problem, derive
the corresponding update formula, and discuss key properties of the method.

###  3.1 A relaxed matrix optimization problem and its solution

As discussed above, if $s_{k}^{T}y_{k}\leq 0$ there is no matrix that is both
positive definite and satisfies the secant equation (5). Therefore, the matrix
optimization problems (7) and (10) that are traditionally used for deriving
quasi-Newton updates are only justified when $s_{k}^{T}y_{k}>0$. However,
suppose that one replaces the secant constraint in those problems with a
penalty term that punishes deviations from the secant condition. It may then
be possible to obtain a positive definite $H_{k+1}$ that is close to
satisfying the secant equation. To this end, we consider the following
relaxation of the matrix optimization problem (10)

$\displaystyle\underset{B}{\text{minimize}}$ $\displaystyle\Upsilon(B):=$
$\displaystyle\
\operatorname{tr}(BH_{k})-\log(\det(BH_{k}))+\alpha_{k}\upsilon(s_{k}-B^{-1}y_{k})$
(11) subject to $\displaystyle B\succ 0.$

The penalty function $\upsilon:\mathbb{R}^{n}\mapsto\mathbb{R}_{\geq 0}$
replaces the hard secant condition and is assumed to be convex with
$\upsilon(0)=0$. The parameter $\alpha_{k}\geq 0$ is used to determine how
large a penalty is imposed for violating the secant equation. There are
various ways that one could choose the form of the penalty term. As we will
show later in this section, letting

$\upsilon(s_{k}-B^{-1}y_{k})=\|s_{k}-B^{-1}y_{k}\|_{B}^{2},$ (12)

where $\|z\|_{B}^{2}=z^{T}Bz$, leads to an explicit update formula that is
invariant under linear changes of variables. Furthermore, since this
particular penalty function can be written as both
$(s_{k}-B^{-1}y_{k})^{T}B(s_{k}-B^{-1}y_{k})$ and
$(Bs_{k}-y_{k})^{T}B^{-1}(Bs_{k}-y_{k})$, it can be seen as a penalty on both
the residual of the secant condition (5) and of the inverse secant condition
(6). The main result of this subsection is given by the following theorem.

######  Theorem 3.1.

For every $\alpha_{k}>0$ and every $H_{k}\succ 0$, there exists a unique
positive definite solution $B^{\star}$ to (11) with the function $\upsilon$
defined in (12). Letting $H_{k+1}=(B^{\star})^{-1}$ leads to the recursive
update

$H_{k+1}=H_{k}+\alpha_{k}s_{k}s_{k}^{T}-\frac{\alpha_{k}}{\gamma_{k}^{2}}(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})^{T},$
(13)

where

$\gamma_{k}=0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}.$

###### Proof.

We begin by showing that the objective function $\Upsilon(B)$ of (11) is
strictly convex. Since $\operatorname{tr}(BH_{k})-\log(\det(BH_{k}))$ has been
proven to be strictly convex on the cone of positive definite matrices in
[13], we only need to show that the penalty term

$\alpha_{k}\|s_{k}-B^{-1}y_{k}\|_{B}^{2}=\alpha_{k}(s_{k}^{T}Bs_{k}-2s_{k}^{T}y_{k}+y_{k}^{T}B^{-1}y_{k}),$

is convex. We see that the first term is linear in $B$ and therefore convex,
while the second is constant. Regarding the third term, taking the inverse of
a positive definite matrix is a convex operation [26] in the sense that for
any positive definite matrices $X$ and $Y$, and scalar $\theta\in[0,1]$, it
holds that

$(\theta X+(1-\theta)Y)^{-1}\preceq\theta X^{-1}+(1-\theta)Y^{-1}.$

From this, it follows that

$y_{k}^{T}(\theta X+(1-\theta)Y)^{-1}y_{k}\leq\theta
y_{k}^{T}X^{-1}y_{k}+(1-\theta)y_{k}^{T}Y^{-1}y_{k},$

showing that the third term is convex too. As a result, $\Upsilon(B)$ is
strictly convex on the cone of positive definite matrices, and the optimal
solution to (11) is unique and must satisfy the first-order optimality
condition

$\nabla\Upsilon(B)=H_{k}-B^{-1}+\alpha_{k}(s_{k}s_{k}^{T}-B^{-1}y_{k}y_{k}^{T}B^{-1})=0.$
(14)

Since the difference between $B^{-1}$ and $H_{k}$ is a symmetric
rank-2-matrix, it is natural to make the following ansatz:

$B^{-1}=H_{k}+\alpha_{k}(s_{k}s_{k}^{T}-u_{k}u_{k}^{T}).$

From Equation (14), we can see that the above equality holds if
$u_{k}=B^{-1}y_{k}$ or $u_{k}=-B^{-1}y_{k}$. Without loss of generality, we
choose

$u_{k}=B^{-1}y_{k}=H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k}-\alpha_{k}(u_{k}^{T}y_{k})u_{k}.$
(15)

To determine $u_{k}$, we first determine the unknown quantity

$u_{k}^{T}y_{k}=y_{k}^{T}u_{k}=y_{k}^{T}H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})^{2}-\alpha_{k}(u_{k}^{T}y_{k})^{2}.$

Solving for $u_{k}^{T}y_{k}$ yields

$u_{k}^{T}y_{k}=-\frac{1}{2\alpha_{k}}\pm\sqrt{\frac{1}{4\alpha_{k}^{2}}+\frac{y_{k}^{T}H_{k}y_{k}}{\alpha_{k}}+(s_{k}^{T}y_{k})^{2}}.$

There are two possible solutions, but for $B^{-1}$ to be positive definite, we
need that $y_{k}^{T}B^{-1}y_{k}=u_{k}^{T}y_{k}>0$, so we use the positive
solution. By Equation (15),

$u_{k}=\frac{1}{1+\alpha_{k}u_{k}^{T}y_{k}}(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k}),$

so plugging in the obtained value of $u_{k}^{T}y_{k}$ gives

$u_{k}=\frac{H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k}}{0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}.$

Thus, setting $H_{k+1}$ to the value of $B^{-1}$ at the critical point leads
to the update formula

$H_{k+1}=H_{k}+\alpha_{k}s_{k}s_{k}^{T}-\frac{\alpha_{k}}{\gamma_{k}^{2}}(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})^{T}$
(16)

with

$\gamma_{k}=0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}.$

It only remains to validate that $H_{k+1}$ is positive definite if $H_{k}$ is.
To this end, we first observe that $H_{k+1}$ is the Schur complement of
$\frac{\gamma_{k}^{2}}{\alpha_{k}}$ in the matrix

$S:=\begin{bmatrix}H_{k}+\alpha_{k}s_{k}s_{k}^{T}&(H_{k}+\alpha_{k}s_{k}s_{k}^{T})y_{k}\\\
y_{k}^{T}(H_{k}+\alpha_{k}s_{k}s_{k}^{T})&\frac{\gamma_{k}^{2}}{\alpha_{k}}\end{bmatrix}.$

Since $\frac{\gamma_{k}^{2}}{\alpha_{k}}$ is positive, $H_{k+1}$ is positive
definite if and only if $S$ is. Furthermore, $S$ is positive definite if and
only if $H_{k}+\alpha_{k}s_{k}s_{k}^{T}$ and its Schur complement in $S$ is.
But $H_{k}+\alpha_{k}s_{k}s_{k}^{T}$ must be positive definite since $H_{k}$
is, so it remains to check whether its Schur complement is positive. The Schur
complement is

$\displaystyle\frac{\gamma_{k}^{2}}{\alpha_{k}}-y_{k}^{T}(H_{k}+\alpha_{k}s_{k}s_{k}^{T})(H_{k}+\alpha_{k}s_{k}s_{k}^{T})^{-1}(H_{k}+\alpha_{k}s_{k}s_{k}^{T})y_{k}$
$\displaystyle=$
$\displaystyle\frac{0.5+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}{\alpha_{k}}$
$\displaystyle-y_{k}^{T}(H_{k}+\alpha_{k}s_{k}s_{k}^{T})y_{k}$
$\displaystyle=$
$\displaystyle\frac{1}{2\alpha_{k}}+\sqrt{\frac{1}{4\alpha_{k}^{2}}+\frac{y_{k}^{T}H_{k}y_{k}}{\alpha_{k}}+(s_{k}^{T}y_{k})^{2}},$

which is indeed positive. The proof is complete. ∎

###  3.2 Properties of the soft QN update

We will now establish four attractive properties of the soft QN update: it
always generates descent directions, it recovers BFGS in the limit, it treats
positive and negative curvature equally, and it is scale invariant. Some of
these properties are immediate, while others require more detailed proofs and
arguments. Nevertheless, we believe that they are all essential features that
contribute to the efficacy of the method.

####  3.2.1 Soft QN always generates descent directions

A remarkable feature of the soft QN update is that it guarantees that $H_{k}$
remains positive definite for any $\alpha_{k}>0$, even if $s_{k}^{T}y_{k}\leq
0$. As mentioned in Section 2, this guarantees that $-H_{k}\nabla f(x_{k})$ is
a descent direction, and that it is the unique minimizer of the quadratic
model function used to approximate the actual loss function.

####  3.2.2 Soft QN recovers BFGS in the limit

The proposed soft QN method is derived as the solution to a matrix
optimization problem that is a relaxed version of the one solved by BFGS. It
is therefore natural to expect that the soft QN update converges to the BFGS
update when $\alpha_{k}\to\infty$, provided that the BFGS update preserves
positive definiteness. The next theorem confirms that this is indeed the case.

######  Theorem 3.2.

If $s_{k}^{T}y_{k}>0$, the right hand side of Equation (13) converges to

$\displaystyle\left(I-\frac{s_{k}y_{k}^{T}}{s_{k}^{T}y_{k}}\right)H_{k}\left(I-\frac{y_{k}s_{k}^{T}}{s_{k}^{T}y_{k}}\right)+\frac{s_{k}s_{k}^{T}}{s_{k}^{T}y_{k}}$

as $\alpha_{k}\to\infty$.

###### Proof.

Equation (13) can be rewritten as

$\displaystyle H_{k+1}=$ $\displaystyle H_{k}+\alpha_{k}s_{k}s_{k}^{T}$
$\displaystyle-\frac{\alpha_{k}}{\gamma_{k}^{2}}(H_{k}y_{k}y_{k}^{T}H_{k}+\alpha_{k}(s_{k}^{T}y_{k})(H_{k}y_{k}s_{k}^{T}+s_{k}y_{k}^{T}H_{k})+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}s_{k}s_{k}^{T})$
$\displaystyle=$
$\displaystyle\left(I-\frac{\alpha_{k}^{2}s_{k}^{T}y_{k}}{\gamma_{k}^{2}}s_{k}y_{k}^{T}\right)H_{k}\left(I-\frac{\alpha_{k}^{2}s_{k}^{T}y_{k}}{\gamma_{k}^{2}}y_{k}s_{k}^{T}\right)$
$\displaystyle+\alpha_{k}\left(1-\frac{\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}{\gamma_{k}^{2}}-\frac{\alpha_{k}^{3}(s_{k}^{T}y_{k})^{2}y_{k}^{T}H_{k}y_{k}}{\gamma_{k}^{4}}\right)s_{k}s_{k}^{T}-\frac{\alpha_{k}}{\gamma_{k}^{2}}H_{k}y_{k}y_{k}^{T}H_{k}$

and should be compared to the BFGS update:

$\displaystyle
H_{k+1}=\left(I-\frac{1}{s_{k}^{T}y_{k}}s_{k}y_{k}^{T}\right)H_{k}\left(I-\frac{1}{s_{k}^{T}y_{k}}y_{k}s_{k}^{T}\right)+\frac{1}{s_{k}^{T}y_{k}}s_{k}s_{k}^{T}.$

To show that the former update converges to the latter, we will show that the
first two terms in our update converge to the corresponding terms in the BFGS
update, while the third term vanishes. We first show that
$\alpha_{k}^{2}s_{k}^{T}y_{k}/\gamma_{k}^{2}$ converges to $1/s_{k}^{T}y_{k}$.
Since

$\gamma_{k}^{2}=0.5+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}},$

the dominant term in the denominator becomes
$\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}$ as $\alpha_{k}\to\infty$, and therefore,

$\lim_{\alpha_{k}\to\infty}\frac{\alpha_{k}^{2}s_{k}^{T}y_{k}}{\gamma_{k}^{2}}=\lim_{\alpha_{k}\to\infty}\frac{\alpha_{k}^{2}s_{k}^{T}y_{k}}{\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}=\frac{1}{s_{k}^{T}y_{k}}.$

Hence, the first term in our update converges to the first term in the BFGS
update. We also see that the third term in our update,
$-\frac{\alpha_{k}}{\gamma_{k}^{2}}H_{k}y_{k}y_{k}^{T}H_{k}$, vanishes as
$\alpha_{k}\to\infty$ since the numerator grows proportionally to $\alpha_{k}$
while the denominator grows as $\alpha_{k}^{2}$. It remains to show that the
second term in our update converges to the second term in the BFGS update,
_i.e._ that

$\lim_{\alpha_{k}\to\infty}\alpha_{k}\left(1-\frac{\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}{\gamma_{k}^{2}}-\frac{\alpha_{k}^{3}(s_{k}^{T}y_{k})^{2}y_{k}^{T}H_{k}y_{k}}{\gamma_{k}^{4}}\right)=\frac{1}{s_{k}^{T}y_{k}}.$

We have that

$\displaystyle\alpha_{k}\left(1-\frac{\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}{\gamma_{k}^{2}}-\frac{\alpha_{k}^{3}(s_{k}^{T}y_{k})^{2}y_{k}^{T}H_{k}y_{k}}{\gamma_{k}^{4}}\right)$
$\displaystyle=$
$\displaystyle\alpha_{k}\Bigg{(}\frac{0.5+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}{\gamma_{k}^{2}}-\frac{\alpha_{k}^{3}(s_{k}^{T}y_{k})^{2}y_{k}^{T}H_{k}y_{k}}{\gamma_{k}^{4}}\Bigg{)}$
$\displaystyle=$
$\displaystyle\frac{\alpha_{k}}{\gamma_{k}^{2}}\left(0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}\right)+\frac{\alpha_{k}^{2}}{\gamma_{k}^{2}}\left(y_{k}^{T}H_{k}y_{k}\left(1-\frac{(s_{k}^{T}y_{k})^{2}\alpha_{k}^{2}}{\gamma_{k}^{2}}\right)\right).$

When $\alpha_{k}\to\infty$, the term $\alpha_{k}^{2}(s_{k}^{T}y_{k})$ will be
the dominant term under the square root, and the first term will approach
$\alpha_{k}^{2}|s_{k}^{T}y_{k}|/\gamma_{k}^{2}$. When $s_{k}^{T}y_{k}>0$, this
is $\alpha_{k}^{2}s_{k}^{T}y_{k}/\gamma_{k}^{2}$ which, as we have already
shown, converges to $1/s_{k}^{T}y_{k}$ as $\alpha_{k}\to\infty$. The second
term will converge to 0, as $\alpha_{k}^{2}/\gamma_{k}^{2}$ converges to
$1/(s_{k}^{T}y_{k})^{2}$ and
$1-(s_{k}^{T}y_{k})^{2}\alpha_{k}^{2}/\gamma_{k}^{2}$ goes to 0. Thus, we
recover the BFGS update formula as $\alpha_{k}\rightarrow\infty$. ∎

####  3.2.3 Soft QN treats positive and negative curvature equally

If we repeat the argument in Theorem 3.2 for $s_{k}^{T}y_{k}<0$, we find that
soft QN tends to a BFGS update where $y_{k}$ is replaced by $-y_{k}$. Thus,
the limiting behavior of the soft QN update as $\alpha_{k}\rightarrow\infty$
is a BFGS update with a reversed sign for $y_{k}$ whenever $s_{k}^{T}y_{k}<0$.
One can see that soft QN treats the curvature as if it was positive instead of
negative in this limiting case. This is no coincidence: the soft QN update
formula (13) does not depend on the sign of $y_{k}$ or $s_{k}$, _i.e._
replacing either $y_{k}$ by $-y_{k}$ or $s_{k}$ by $-s_{k}$ would result in
the same update. This means that any quasi-Newton algorithm based on this
update is only concerned with the absolute value of the curvature in the
direction $s_{k}$ and not the sign, and will re-scale the gradient to take
longer steps in directions with lower absolute curvature.

To see why this could be advantageous, consider the argument laid out in [12]
to motivate a modification of Newton’s method from trust-region principles.
The authors consider an iterative algorithm $x_{k+1}=x_{k}+p_{k}$ where
$p_{k}$ is the solution to the trust-region problem

$\displaystyle\underset{p}{\text{minimize}}$ $\displaystyle
f(x_{k})+p^{T}\nabla f(x_{k})$ (17) subject to $\displaystyle
p^{T}|\nabla^{2}f(x_{k})|p\leq c.$

Here, $|\nabla^{2}f(x_{k})|$ denotes the matrix obtained by taking the Hessian
evaluated at $x_{k}$ and setting all its eigenvalues to their absolute values
while retaining its eigenvectors, while $c>0$ is the trust-region radius. As
usual for trust-region problems, the objective function is replaced by a model
based on its Taylor expansion, although trust-region problems typically used
both first- and second-order terms. An upper bound of the second-order term is
instead used to determine the shape of the trust-region. In contrast to trust-
region methods that just use a 2-norm constraint, this formulation takes into
account that the first-order model becomes inaccurate quicker in directions of
high curvature, whether positive or negative. Solving this trust-region update
gives that $p_{k}$ is proportional to $-|\nabla^{2}f(x_{k})|^{-1}\nabla
f(x_{k})$. The method proposed in [12], called saddle-free Newton, is based on
this insight. A similar idea is the basis for the non-convex Newton method
proposed in [27]. While the practical methods proposed in these works make
additional modifications to the iteration
$x_{k+1}=x_{k}-|\nabla^{2}f(x_{k})|^{-1}\nabla f(x_{k})$ to be scalable in
higher dimensions, both works argue that it is advantageous to move further in
directions of low absolute curvature in order to escape from saddle points.

While both saddle-free Newton and non-convex Newton, in contrast to the
regular Newton method, escape saddle points by moving far in directions
corresponding to eigenvalues with a low absolute value, they are unable to
detect directions of low curvature that are linear combinations of
eigenvectors with eigenvalues of opposite sign. Following the trust-region
principle, it would be desirable for a method to identify such directions and
re-scale the gradient accordingly to obtain a suitable step. Here, a quasi-
Newton method that updates its Hessian approximation based on the absolute
value of the curvature can have an advantage. To illustrate this, we consider
the 2-dimensional toy-example of optimizing

$\displaystyle f(x,y)$
$\displaystyle=(x-0.7)^{2}((x+0.7)^{2}+0.1)+(y+0.7)^{2}((y-0.7)^{2}+0.1).$

Figure 1:  The soft QN method detects a diagonal direction of low curvature
and takes a large step along it. The saddle-free Newton method is unable to
perceive the low curvature in that direction and moves further toward the
saddle point before changing its course.

This function has one local maximum, four saddle points and four local minima,
with one of them also being a global minimum. The function is separable in the
two variables, so the Hessian is diagonal. In Figure 1, we show plots of 500
iterations of $x_{k+1}=x_{k}-\eta|\nabla^{2}f(x_{k})|^{-1}\nabla f(x_{k})$
(denoted as saddle-free Newton), against iterations of $x_{k+1}=x_{k}-\eta
H_{k}\nabla f(x_{k})$ where $H_{k}$ is updated according to Equation (13) and
$H_{0}$ is $|\nabla^{2}f(x_{0})|^{-1}$. We set $x_{0}=(-0.05,0.08)$, the step-
size $\eta=0.01$ and the penalty parameter $\alpha_{k}=8\cdot 10^{5}$. The
parameters have been chosen to highlight the behavior that we want to
illustrate. While the step-size of Newton-type methods is typically set to 1
or determined by line-search, the short step-size and high penalty parameter
ensure that the methods will take long steps if and only if they discover a
direction of low curvature. This phenomenon is illustrated in Figure 1: at the
point $(0.543,0.0574)$, the Hessian has the values $1.78$ and $-1.72$ on its
diagonal. The eigenvalues are of opposite signs but very close in magnitude,
meaning that saddle-free Newton treats the Hessian almost as a scaled identity
matrix. Soft QN, on the contrary, detects that the curvature is close to 0 for
a diagonal direction, takes a long step in that direction to end up at
$(0.827,-0.230)$ in the next iteration. This sets it far ahead of saddle-free
Newton and allows it to find the local minimum within 500 iterations.

####  3.2.4 Soft QN is scale invariant

Scale invariance, _i.e._ invariance under linear transformations of the
variables, is a desirable property of optimization algorithms since it
guarantees that the algorithm is insensitive to ill-conditioning of the
Hessian. To understand why, consider a strongly convex problem with an ill-
conditioned Hessian. By a linear transformation of variables, we can transform
such a problem into an equivalent one with an arbitrary positive definite
Hessian. If the algorithm is scale invariant, then the iterations will not be
affected by the variable transformation, and the method will perform equally
well on the original (ill-conditioned) problem formulation as the transformed
(well-conditioned) one. The next result shows that a scale-invariance property
holds for the soft QN update.

######  Proposition 3.3.

The soft QN update is scale-invariant, i.e., under a change of variables
$\tilde{x}_{k}=Ax_{k}$, where $A$ is an invertible square matrix,
$\tilde{H}_{k}=AH_{k}A^{T}$ implies that $\tilde{H}_{k+1}=AH_{k+1}A^{T}$,
where $\tilde{H}_{k}$ is the Hessian inverse approximation associated with the
variables $\tilde{x}_{k}.$

###### Proof.

Defining $\tilde{y}_{k}$ and $\tilde{s}_{k}$ from $\tilde{x}_{k}$ analogously
to how $s_{k}$ and $y_{k}$ are defined from $x_{k}$, we have that
$\tilde{s}_{k}=As_{k}$ and $\tilde{y}_{k}=A^{-T}y_{k}$. Since
$\tilde{y_{k}}^{T}\tilde{s}_{k}=y_{k}A^{-1}As_{k}=y_{k}^{T}s_{k}$ and
$\tilde{y}_{k}^{T}\tilde{H}_{k}\tilde{y}_{k}=y_{k}A^{-1}AH_{k}A^{T}A^{-T}y_{k}=y_{k}^{T}H_{k}y_{k}$,
so the quantity $\gamma_{k}$ from Equation (13) is scale invariant. Therefore,

$\displaystyle\tilde{H}_{k+1}=$
$\displaystyle\tilde{H}_{k}+\alpha_{k}\tilde{s}_{k}\tilde{s}_{k}^{T}-\frac{\alpha_{k}}{\gamma_{k}^{2}}(\tilde{H}_{k}\tilde{y}_{k}+\alpha_{k}(\tilde{s}_{k}^{T}\tilde{y}_{k})\tilde{s}_{k})(\tilde{H}_{k}\tilde{y}_{k}+\alpha_{k}(\tilde{s}_{k}^{T}\tilde{y}_{k})\tilde{s}_{k})^{T}$
$\displaystyle=$ $\displaystyle AH_{k}A^{T}+\alpha_{k}As_{k}s_{k}^{T}A^{T}-$
$\displaystyle\frac{\alpha_{k}}{\gamma_{k}^{2}}(AH_{k}A^{T}A^{-T}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})As_{k})(AH_{k}A^{T}A^{-T}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})As_{k})^{T}$
$\displaystyle=$ $\displaystyle
A(H_{k}+\alpha_{k}s_{k}s_{k}^{T}-\alpha_{k}(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})^{T})A^{T}$
$\displaystyle=$ $\displaystyle AH_{k+1}A^{T}.$

∎

The preservation of scale-invariance is a property that soft QN shares with
all quasi-Newton methods in the Broyden class [25]. However, if $H_{0}$ is not
scale-invariant then the subsequent Hessian inverse approximations will not be
scale-invariant either. While it is possible to choose a scale-invariant
$H_{0}$, e.g. setting it equal to the true Hessian at the initial point, this
is typically computationally expensive and it is common to simply set $H_{0}$
to the identity matrix, which corresponds to $\tilde{H}_{0}=AA^{T}$ under the
change of variables $\tilde{x}_{k}=Ax_{k}$. However, as the iteration number
increases, the influence of the initial matrix diminishes, and one can argue
that the Hessian inverse approximation becomes “more scale invariant” and
thereby more resistant to ill-conditioning. This is confirmed empirically with
traditional quasi-Newton methods outperforming gradient descent on ill-
conditioned problems [25].

###  3.3 Comparison with SP-BFGS

The SP-BFGS method [16] is similar to soft QN, but it is based on solving a
relaxed version of problem (7) instead of problem (10). Specifically, it is
derived as the solution to

$\displaystyle\underset{H}{\text{minimize}}$
$\displaystyle\|W_{k}^{(1/2)}(H-H_{k})W_{k}^{(1/2)}\|_{F}^{2}+\frac{\beta_{k}}{2}\|W_{k}^{(1/2)}(Hy_{k}-s_{k})\|_{2}^{2}$
(18) subject to $\displaystyle H=H^{T}$

for some strictly positive penalty parameter $\beta_{k}>0$. In this
formulation, the inverse secant condition is weighted by the (principal square
root of the) weight matrix $W_{k}$, which can be any positive definite matrix
that satisfies the secant condition. This weighting ensures scale-invariance.
For any such $W_{k}$, the solution of this problem leads to the update formula

$H_{k+1}=\left(I-\omega_{k}s_{k}y_{k}^{T}\right)H_{k}\left(I-\omega_{k}y_{k}s_{k}^{T}\right)+\omega_{k}\left(\dfrac{\pi_{k}}{\omega_{k}}+(\pi_{k}-\omega_{k})y_{k}^{T}H_{k}y_{k}\right)s_{k}s_{k}^{T},$
(19)

where

$\pi_{k}=\dfrac{1}{s_{k}^{T}y_{k}+\frac{1}{\beta_{k}}}\ ,\
\omega_{k}=\dfrac{1}{s_{k}^{T}y_{k}+\frac{2}{\beta_{k}}}.$

By Lemma 1 in [16], the update (19) preserves positive definiteness when
$s_{k}^{T}y_{k}>-\beta_{k}^{-1}$. This means that it can still yield positive
definite matrices when $s_{k}^{T}y_{k}<0$, as long as the penalty parameter
$\beta_{k}$ is not too large, and the update formula can be used as long as
care is taken to choose $\beta_{k}$ appropriately. However, as discussed in
Section 2, there is no matrix $W_{k}$ that satisfies the secant condition when
$s_{k}^{T}y_{k}<0$, so the matrix optimization problem (18) is not justified
and it is not clear in what sense $H_{k+1}$ is an optimal update of $H_{k}$.
In contrast, $H_{k+1}$ given by the proposed update formula (13) is positive
definite for any value of $\alpha_{k}>0$, and it always exists as an optimal
solution to the problem (11).

When the penalty parameters of soft QN and SP-BFGS approach infinity, the
behavior of the update formulas is similar, but with one important difference:
While the SP-BFGS update always converges to the BFGS update, soft QN
converges to a modified update where the sign of $y_{k}$ is changed if
$s_{k}^{T}y_{k}$ is negative (cf. the discussion in Section 3.2.2).

##  4 Convergence

In this section, we demonstrate how the soft QN update can be integrated into
a quasi-Newton method with guaranteed convergence rate for strongly convex
problems also in the presence of bounded errors in gradient and function
evaluations. The integration relies on a uniform bound on the Hessian inverse
approximation generated by soft QN.

###  4.1 Bounding the Hessian inverse approximation

In this subsection, we derive a method to bound $H_{k}$ in a way that takes
into account its size in the current iteration, allowing larger changes to the
matrix if it is far from a desired upper or lower bound. We do so by having a
different penalty parameter $\alpha_{k}$ in each iteration $k$ and choosing it
based on upper and lower bounds that depend on a size measure of $H_{k}$.
Several convergence proofs, e.g. [32], [24], and [16], rely on upper and lower
bounds on the eigenvalues of $H_{k}$, and we will do that here too. We denote
the maximum eigenvalue of a matrix $A$ with real-valued eigenvalues by
$\lambda_{\max}(A)$ and the minimum eigenvalue by $\lambda_{\min}(A)$.

######  Lemma 4.1.

Given $\psi I\preceq H_{0}\preceq\Psi I$ with some constants
$0<\psi\leq\Psi<\infty$, if the penalty parameter $\alpha_{k}$ satisfies

$\alpha_{k}\leq\min\left\\{\frac{\lambda_{\min}(H_{k})-\psi}{(\|s_{k}\|_{2}+\|H_{k}y_{k}\|_{2})^{2}},~{}\dfrac{\Psi-\lambda_{\max}(H_{k})}{\|s_{k}\|_{2}^{2}}\right\\},$
(20)

then the sequence $\\{H_{k}\\}_{k}$ generated by the proposed soft QN
satisfies

$\psi I\preceq H_{k}\preceq\Psi I,~{}\forall k.$ (21)

###### Proof.

From Equation (13), we have

$\lambda_{\max}(H_{k+1})\leq\lambda_{\max}(H_{k})+\alpha_{k}\|s_{k}\|_{2}^{2}.$

Hence, we can ensure that the eigenvalues are smaller than some constant
$\Psi$ by letting

$\alpha_{k}\leq\dfrac{\Psi-\lambda_{\max}(H_{k})}{\|s_{k}\|_{2}^{2}}.$

To lower bound $H_{k}$, we observe from Equation (13) that

$\lambda_{\min}(H_{k+1})\geq\lambda_{\min}(H_{k})-\alpha_{k}\|u_{k}\|_{2}^{2}$

with

$u_{k}=\frac{H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k}}{0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}.$

Since $u_{k}$ depends on $\alpha_{k}$, we want to eliminate $\|u_{k}\|_{2}$
from the above expression to derive a practical bound. We can do so by
observing that

$\displaystyle\|u_{k}\|_{2}$
$\displaystyle=\frac{\|H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k}\|_{2}}{0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}$
$\displaystyle\leq\frac{\|H_{k}y_{k}\|_{2}+\alpha_{k}\|(s_{k}^{T}y_{k})s_{k}\|_{2}}{0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}$
$\displaystyle\leq\frac{\|H_{k}y_{k}\|_{2}}{1}+\frac{\alpha_{k}|s_{k}^{T}y_{k}|\|s_{k}||_{2}}{\sqrt{\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}}=\|H_{k}y_{k}\|_{2}+\|s_{k}\|_{2}.$

With this bound, $\lambda_{\min}(H_{k+1})$ is guaranteed to stay above $\psi$
if

$\alpha_{k}\leq\frac{\lambda_{\min}(H_{k})-\psi}{(\|s_{k}\|_{2}+\|H_{k}y_{k}\|_{2})^{2}}.$

This completes the proof. ∎

While Lemma 4.1 shows that $H_{k}$ will remain bounded if $\alpha_{k}$ is
small enough, the computational demands are something to consider if we want
to use it in practice. While $\lambda_{\min}(H_{k})$ and
$\lambda_{\max}(H_{k})$ can be computed with power iteration, doing so might
be slow if there are several eigenvalues that are close in magnitude. However,
note that one can replace the maximum eigenvalue of $H_{k}$ by an upper bound
and the minimum eigenvalue by a lower bound and still ensure boundedness of
$H_{k}$. There are several ways that such bounds can be calculated [21] [33],
but to ensure that computing it does not overshadow other operations performed
in the algorithm, e.g. the multiplication of $H_{k}$ with the gradient, we
want the evaluation of the bound to scale at most as $O(n^{2})$, where $n$ is
the number of variables in the problem. An upper bound that satisfies this and
works for any $n\times n$ complex matrix with real eigenvalues $A$ is

$\lambda_{\max}(A)\leq\frac{\operatorname{tr}(A)}{n}+\sqrt{\dfrac{\frac{\operatorname{tr}(A^{2})}{n}-\big{(}\frac{\operatorname{tr}(A)}{n}\big{)}^{2}}{n-1}}.$

An inexpensive way of enforcing a lower bound of the lowest eigenvalue is to
add a small bias term $\lambda I$ to the nominal $H_{k}$ before multiplying it
with the gradient. Since $H_{k}$ is guaranteed to be positive definite no
matter the values of $s_{k}$ and $y_{k}$, $\lambda_{\min}(H_{k}+\lambda
I)\geq\lambda$.

###  4.2 Convergence guarantees

With the uniform bound on the Hessian inverse approximation established in
Lemma 4.1, it is safe to use the Hessian inverse approximations generated by
soft QN to construct approximate Newton directions for solving the
minimization problem (2) in the presence of general bounded noise. More
specifically, assume that we want to optimize a function $\phi(x)$ but only
have access to corrupted measurements of the function and the gradient

$f(x):=\phi(x)+n_{f}(x),~{}g(x):=\nabla\phi(x)+n_{g}(x),$ (22)

where $n_{f}(x)$ and $n_{g}(x)$ are noises in the function and gradient,
respectively. We make the following assumption on the noises.

######  Assumption 1.

There exist positive constants $e_{f}$ and $e_{g}$ such that

$|n_{f}(x)|\leq e_{f},~{}\|n_{g}(x)\|_{2}\leq e_{g},\quad\forall x.$

The soft QN algorithm for solving problem (2) is summarized in Algorithm 1.

Algorithm 1 Soft QN Algorithm

1:Set $\eta_{k}$, $H_{0}\succ 0$, $x_{0}$, and $g_{0}$

2:for $k$=$0$, $1$, …, $K-1$ do

3: Update $p_{k}=-H_{k}g_{k}$

4: Update $x_{k+1}=x_{k}+\eta_{k}p_{k}$

5: Update $g_{k+1}=g(x_{k+1})$

6: Update $s_{k}=x_{k+1}-x_{k}$

7: Update $y_{k}=g_{k+1}-g_{k}$

8: Choose penalty parameter $\alpha_{k}$ satisfying (20)

9: Update
$\gamma_{k}=0.5+\sqrt{0.25+\alpha_{k}y_{k}^{T}H_{k}y_{k}+\alpha_{k}^{2}(s_{k}^{T}y_{k})^{2}}$

10: Update
$H_{k+1}=H_{k}+\alpha_{k}s_{k}s_{k}^{T}-\frac{\alpha_{k}}{\gamma_{k}^{2}}(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})(H_{k}y_{k}+\alpha_{k}(s_{k}^{T}y_{k})s_{k})^{T}$

11:end for

12:Output $x_{K}$

For twice continuously differentiable, $L$-smooth and $\mu$-strongly convex
functions, with bounded noise on gradient and function evaluations, the
authors of [16] show that general quasi-Newton methods with bounded Hessian
inverse approximations converge linearly to a neighborhood of the optimal
solution when the step-size is fixed or is generated by line-search. As these
theorems are also applicable here, we state them and discuss how they apply to
our quasi-Newton method.

######  Assumption 2.

The gradient $\nabla\phi$ is Lipschitz continuous with Lipschitz constant $L$;
namely for any $x,y\in\mathbb{R}^{n}$,

$\|\nabla\phi(x)-\nabla\phi(y)\|\leq L\|x-y\|.$

######  Assumption 3.

The function $\phi$ is $\mu$-strongly convex, i.e., there exists a positive
constant $\mu$ such that for all $x,y\in\mathbb{R}^{n}$:

$\phi(y)\geq\phi(x)+\nabla\phi(x)^{T}(y-x)+\frac{\mu}{2}\|x-y\|^{2}.$

The following theorem provides convergence guaranteed under a fixed step-size
$\eta_{k}=\eta$.

######  Theorem 4.2.

(Fixed step-size $\eta_{k}=\eta$) Consider Algorithm 1. Under Assumptions 1–3,
if the step-size $\eta$ satisfies

$0<\eta\leq\frac{\varepsilon}{L\Psi^{2}}<\frac{\psi}{L\Psi}\left[\frac{A}{((1+A)\Psi+\psi)}\right],$

where $0<\varepsilon<\frac{A\psi\Psi}{((1+A)\Psi+\psi)}$ and $A>0$ is some
constant, then for any $x_{k}\notin\mathcal{N}_{1}\left(\psi,\Psi,A\Psi
e_{g}\right)$, where

$\displaystyle\mathcal{N}_{1}\left(\psi,\Psi,A\Psi
e_{g}\right)=\left\\{x\mid\phi(x)\leq\phi^{\star}+\frac{1}{2\mu}\left(\frac{(1+A)\Psi
e_{g}}{\psi}\right)^{2}\right\\},$

we have

$\displaystyle\phi(x_{k+1})-\left(\phi(x^{\star})+\frac{1}{2\mu}\left(\frac{(1+A)\Psi
e_{g}}{\psi}\right)^{2}\right)$ $\displaystyle\leq$
$\displaystyle(1-\eta\varepsilon\mu)\left(\phi(x_{k})-\left(\phi(x^{\star})+\frac{1}{2\mu}\left(\frac{(1+A)\Psi
e_{g}}{\psi}\right)^{2}\right)\right).$

###### Proof.

See [16, Theorem 4]. ∎

Besides a fixed step-size, we can also use line-search to generate an
$\eta_{k}$ for iteration $k$. To account for the noise in the function
evaluations, we use the relaxed Armijo condition introduced in [16]

$f(x_{k}+\eta_{k}p_{k})\leq
f(x_{k})+\eta_{k}cp_{k}^{T}g(x_{k})+2\varepsilon_{\rm tol}.$ (23)

Here, $c>0$, $\varepsilon_{\rm tol}>0$, and the step-size $\eta_{k}$ is the
maximum value in $\left\\{\tau^{j}\eta_{0}\mid j=0,1,2,\ldots\right\\}$ that
satisfies the relaxed Armijo condition (23) for some $\eta_{0}>0$ and
$\tau\in(0,1)$. When we use such step-sizes $\eta_{k}$ in Algorithm 1 we get a
soft QN Algorithm with line-search that has the following guarantees.

######  Theorem 4.3.

(Line-search step-size $\eta_{k}$ satisfying (23)) Consider Algorithm 1. Under
Assumptions 1–3, let $\varepsilon_{\rm tol}>\epsilon_{f}$ and

$0<c\leq\frac{\varepsilon}{2\Psi}\frac{\left(1-\frac{\psi}{(1+A)\Psi}\right)}{\left(1+\frac{\psi}{(1+A)\Psi}\right)}<\frac{\psi}{2}\frac{A((1+A)\Psi-\psi)}{((1+A)\Psi+\psi)^{2}},$

where $0<\varepsilon<\frac{A\psi\Psi}{((1+A)\Psi+\psi)}$ and $A>0$ is some
constant. Then, for any
$x_{k}\notin\mathcal{N}_{1}\left(\psi,\Psi,A\Psi{e}_{g}\right)\cup\mathcal{N}_{2}(\bar{\eta})$,
where

$\mathcal{N}_{2}(\bar{\eta})=\left\\{x\mid\phi(x)\leq\phi^{\star}+\bar{\eta}\right\\},~{}\bar{\eta}=\frac{\Psi
L}{c\tau\left(1-\frac{\psi}{(1+A)\Psi}\right)^{2}\varepsilon\mu}\left(\varepsilon_{\rm
tol}+\epsilon_{f}\right),$

we have

$\phi(x_{k+1})-\left(\phi(x^{\star})+\bar{\eta}\right)\leq\left(1-\frac{2c\tau\left(1-\frac{\psi}{(1+A)\Psi}\right)^{2}}{\Psi
L}\varepsilon\mu\right)\left(\phi(x_{k})-\left(\phi(x^{\star})+\bar{\eta}\right)\right).$

###### Proof.

See [16, Theorem 5]. ∎

Theorem 4.2 establishes that the soft QN iterates are guaranteed to decrease
the objective value at a linear rate to a neighborhood around the optimal
value. However, the presence of noise in the gradient estimates eventually
impedes further progress. From the definition of ${\mathcal{N}}_{1}$, we see
that the upper bound on the residual error is proportional to $e_{g}^{2}$.
Since $\varepsilon<\psi$, it holds that
$\eta\varepsilon\mu\leq\frac{\psi}{L\Psi^{2}}\psi\mu=\frac{\mu}{L}\frac{\psi^{2}}{\Psi^{2}}$
and we can infer a convergence factor of at least
$1-O(\frac{\mu}{L}\frac{\psi^{2}}{\Psi^{2}})$. Moreover, there exists a trade-
off between the convergence rate and the proximity to the optimal value that
can be influenced by adjusting the free parameter $A$. Specifically, reducing
$A$ results in a smaller neighborhood, but it also reduces the admissible
step-size $\eta$ and consequently yields slower convergence. In particular, in
the absence of gradient computation noise, i.e., when $e_{g}=0$, linear
convergence to the optimal value can be attained.

The bounds presented in Theorem 4.3 are more complicated to interpret, due to
the influence of noise on both gradient and function evaluations. When
$\mathcal{N}_{2}$ dominates (i.e., $x_{k}\in\mathcal{N}_{2}$ but
$x_{k}\notin\mathcal{N}_{1}$), the iterates produced by the algorithm generate
objective function values that converge linearly to a neighborhood surrounding
the optimal value that is inescapable due to the presence of $\epsilon_{f}$.
However, in scenarios where $\mathcal{N}_{1}$ dominates, the algorithm lacks
theoretical guarantees to further converge towards $\mathcal{N}_{2}$.
Additionally, the magnitude of the function value decreases ensured by the
line-search impacts the convergence factors and the guaranteed precision.
Decreasing the parameter $c$ and increasing $\varepsilon_{\rm tol}$ relaxes
the Armijo condition more but results in a slower rate of convergence and a
larger neighborhood.

######  Remark 1.

Unlike most existing stochastic QN methods that consider noise from random
sampling [23, 32, 24, 36, 10], this paper addresses scenarios involving more
general bounded noise. Due to the unclear characteristics of general noise, we
assume only boundedness in Assumption 1, resulting in convergence results, as
shown in Theorems 4.2 and 4.3, which are no better than first-order methods.
While some studies assume unbiased gradient estimates and employ variance
reduction techniques to mitigate the impact of noise on convergence in the
case of random sampling [24, 36, 10], these techniques are not applicable for
general bounded noise. Moreover, the use of stochastic gradients to
simultaneously update a Hessian inverse approximation and define the search
direction of the quasi-Newton algorithm results in a situation where $H_{k}$
and $g_{k}$ are correlated, which poses significant challenges for analysis.
One way to ensure independence of $g_{k}$ and $H_{k}$ is to calculate
additional gradients for generating $H_{k}$, instead of using the same
gradients $g_{k}$ that are used to define the search directions and impose
additional assumptions such as independence of gradient noises at different
evaluations. If $H_{k}$ and $g_{k}$ are uncorrelated and upper and lower
bounds on $H_{k}$ are enforced, a soft QN algorithm can be shown to fit the
framework for stochastic quasi-Newton methods of [32]. Hence, it can be shown
to be convergent in a setting with a non-convex loss function and gradient
noise bounded only in variance. However, using different gradient evaluations
for updating $H_{k}$ and computing $g_{k}$ increases the number of gradient
evaluations and limits the applicable noise scenarios. As we will show in the
numerical experiments, our soft QN matches the performance of stochastic
gradient descent with random sampling noise (see Fig. 2) and significantly
outperforms it with additive noise (see Fig. 3).

##  5 Numerical results

In this section, we assess the effectiveness of the proposed soft QN method on
logistic regression problems, randomly generated quadratic problems, and
problems from the CUTEst collection [14]. The type of noise is different in
all these tests: In the case of logistic regression, the noise arises due to
random sampling of data. For both the quadratic problems and the CUTEst
problems, artificial random noise is generated and added to the gradient. For
the CUTEst problems, noise is added to the objective function as well. For the
quadratic problems, the noise is sampled from the standard normal
distribution, while for the CUTEst problems, it is sampled uniformly from a
hypersphere. For all the tests, we used a constant $\alpha_{k}$ for soft QN to
keep it simple. For SP-BFGS, we used methods suggested in [16] to determine
$\beta_{k}$. Different methods for step-size selection are also used in the
three different tests. For all the quasi-Newton methods used, the initial
Hessian inverse approximation is an identity matrix. Since all the problems
are stochastic, the algorithm does not have access to the true objective
function and its gradient, but approximations of them. When we describe the
stochastic problems below, we use the same notation as in (22): $\phi(x)$ is
the objective function that we want to optimize, $\nabla\phi(x)$ is its
gradient, while $f(x)$ and $g(x)$ are the stochastic approximations of
$\phi(x)$ and $\nabla\phi(x)$ available to the solver.

###  5.1 Logistic regression

In the first setting, we consider noise that originates from random sampling.
To do so, we let the algorithms solve a logistic regression problem

$\displaystyle\underset{x\in\mathbb{R}^{n+1}}{\mbox{minimize}}~{}\phi(x)=\frac{1}{N}\sum_{i=1}^{N}\phi(x;\\{z_{i},y_{i}\\})+\rho\|x_{1:n}\|^{2},$
$\displaystyle\text{with}~{}\phi(x;\\{z_{i},y_{i}\\}):=\log(1+\exp(-y_{i}(x_{0}+x_{1:n}^{T}z_{i}))).$

Here, $x\in\mathbb{R}^{n+1}$ is the decision variable whose first component,
$x_{0}\in\mathbb{R}$, is the model weight for the intercept and the remaining
components, $x_{1:n}\in\mathbb{R}^{n}$, are feature weights. Further,
$(z_{i},y_{i})$ are the feature-label pairs in the training data, where
$z_{i}\in\mathbb{R}^{n}$ are features and $y_{i}\in\\{-1,1\\}$ are the
associated labels. The scalar $\rho>0$ is a regularization parameter.

To compute gradient approximations, we randomly sample a subset
$\mathcal{S}\subset\\{1,...,N\\}$ with cardinality $|\mathcal{S}|$ and
approximate the full gradient by the stochastic gradient $g(x)$ given by
$g(x):=\frac{1}{|\mathcal{S}|}\sum_{\mathcal{S}_{i}\in\mathcal{S}}\nabla\phi(x;(y_{\mathcal{S}_{i}},z_{\mathcal{S}_{i}}))$.
We use the data set ijcnn1 from libsvm [9], where $n=22$ and $N=49990$. The
data was pre-processed using feature-wise 2-norm normalization.

We compare our soft QN with SP-BFGS, SGD, and BFGS with gradients computed by
our stochastic approximation, which we refer to as stochastic BFGS. We set
$|\mathcal{S}|=1000$ and $\rho=0.1$. For soft QN, we set the penalty parameter
$\alpha_{k}=0.5$. For SP-BFGS, we set the penalty parameter
$\beta_{k}=0.1\|s_{k}\|_{2}+10^{-10}$. As suggested in [16], such $\beta_{k}$
can adaptively relax the secant equation when $s_{k}^{T}y_{k}$ is dominated by
noise (corresponding to small $\|s_{k}\|_{2}$ ) and tighten the secant
equation when $s_{k}^{T}y_{k}$ is dominated by curvature information
(corresponding to large $\|s_{k}\|_{2}$). Both the constant value of
$\alpha_{k}$ and the coefficient of $\|s_{k}\|_{2}^{2}$ in $\beta_{k}$ were
determined by tuning experiments conducted for one particular realization of
the random sampling. For all the algorithms, we set the step-size
$\eta_{k}=0.1$. We conduct a Monte Carlo simulation with 100 trials, in which
we provide the algorithms with only sampling-based gradient approximations,
while also recording the actual gradient norm. For each trial and each method,
the iterate vector is initialized to be a vector of all zeros.

Figure 2:  Logistic regression experiments. The 10-logarithm of the gradient
norm is plotted against iteration number. Solid lines indicate the mean over
all instances in the Monte Carlo simulation, while shaded areas represent
three standard deviation confidence intervals.

Fig. 2 depicts the average 10-logarithm of the gradient norm for a certain
iteration number, with shaded areas around the graphs representing values
within three sample standard deviations of the estimator. It shows that all
the other algorithms outperform stochastic BFGS. Although stochastic BFGS
manages to decrease the gradient, the adaptions made to derive SP-BFGS and
soft QN led to an improvement. We did not find a statistically significant
difference between the performance of SP-BFGS, soft QN, and SGD in this test.
This stands in contrast to the results of the subsequent experiments, which
show a clear advantage for soft QN. We believe that this is due to the
regularization and the pre-processing of the data, which makes the logistic
regression problem well-conditioned and limits the advantage of methods that
take the curvature of the objective function into account with Hessian inverse
approximations.

###  5.2 Quadratic problems

In the second setting, we test soft QN on problems with noise that has bounded
variance, added in the gradient evaluations. We generate quadratic problems on
the form

$\underset{x\in\mathbb{R}^{n}}{\mbox{minimize}}~{}\phi(x):=\frac{1}{2}x^{T}Hx+b^{T}x.$
(24)

The procedure for generating these problems is based on [19]. We construct the
matrix $H\in{\mathbb{R}}^{n\times n}$ by first defining its eigenvectors from
the columns of the factor Q in the QR decomposition of a random matrix where
each element follows a standard normal distribution. We then fix the minimum
eigenvalue of $H$ to be 0.01 and the maximum eigenvalue to be 1, and draw the
remaining $n-2$ eigenvalues uniformly from the interval $[0.01,1]$. Finally,
we choose the vector $b\in\mathbb{R}^{n}$ in such a way that the optimal
solution to the problem (24) is ${\bf 1}$, representing a vector with all
elements equal to 1. In our tests, we use $n=100$ variables.

We consider a scenario where the algorithm has access to a gradient
approximation $g(x)$ that is perturbed by noise $n_{g}\sim{\mathcal{N}}(0,I)$,

$g(x)=\nabla\phi(x)+n_{g}.$

As shown in [28], for unconstrained quadratic problems with additive normal
noise on the gradients, Newton’s method with the exact Hessian, using a step-
size $\eta_{k}=1/k$ where $k$ is the iteration number, is an asymptotically
optimal algorithm. Therefore, we use Newton’s method as a baseline, utilizing
a precise Hessian matrix but with noisy gradient $g(x)$. In addition, we
compare the proposed soft QN with SGD, stochastic BFGS, and SP-BFGS. Note that
the choice of additive normal noise is made to ensure the optimality of
Newton’s method. While the noise components are theoretically unlimited in
this case, they are limited in any practical implementation. Tests with
bounded noise will be presented in the next section.

For soft QN, we set the penalty parameter $\alpha_{k}=10^{-4}$. For SP-BFGS,
we set the penalty parameter $\beta_{k}=10^{-2}$ if $s_{k}^{T}y_{k}\geq 0$ and
$\beta_{k}=\frac{-0.9}{s_{k}^{T}y_{k}}$ if $s_{k}^{T}y_{k}<0$. Here, we apply
the SP-BFGS relaxation of the positive curvature condition [16, Equation (51)]
when $s_{k}^{T}y_{k}<0$. The constant value of $\alpha_{k}$ and the value of
$\beta_{k}$ used when $s_{k}^{T}y_{k}\geq 0$ were determined by tuning
experiments. For stochastic BFGS, we implement the fail-safe measure that if
$s_{k}^{T}y_{k}<0$ due to random noise, we do not update the Hessian inverse
approximation. For all algorithms, we use a diminishing step-size
$\eta_{k}=\frac{1}{k}$. We run a Monte Carlo simulation with 100 trials, with
the initial variable vector set to ${\bf 0}$ in each trial. Note that from the
way that the problems are generated, the squared norm of the initial true
gradient is upper bounded by $n$, and the expected squared norm of the noise
vector is $n$ throughout the iterations. Thus, noise has a significant impact
on the gradient estimate even in the initial iteration. This explains why the
value of $\alpha_{k}$ obtained in the initial tuning experiments was
significantly lower than the one obtained for logistic regression: a small
penalty parameter is needed in order to avoid overfitting to the noise. The
results of our Monte Carlo simulations are shown in Fig. 3. Each graph shows
the average 10-logarithm of the normalized suboptimality, i.e.
$\log_{10}((\phi(x_{k})-\phi({\bf 1}))/(\phi({\bf 0})-\phi({\bf 1})))$, with
shaded areas representing three standard deviation confidence intervals for
the mean estimates.

Newton’s method performs the best due to its use of accurate Hessians;
stochastic BFGS actually increases the suboptimality by several orders of
magnitude, as its use of noisy gradients fails to effectively capture
curvature information. SGD converges slowly due to its reliance on only first-
order information, and the problem is generated to have a high condition
number. SP-BFGS and the proposed soft QN are faster than SGD since the second-
order information accelerates convergence. As the confidence intervals do not
overlap, we see that soft QN outperforms SP-BFGS with statistical
significance.

Figure 3:  Experiments on quadratic problems. The 10-logarithm of the
normalized suboptimality is plotted against iteration number. Solid lines show
the mean over all instances in the Monte Carlo simulation, shaded areas
represent three standard deviation confidence intervals.

###  5.3 CUTEst problems

In the third setting, we consider bounded noise in the function and gradient
evaluations. Specifically, we let the solver have access to $f(x)$ and $g(x)$,
where

$f(x):=\phi(x)+n_{f}(x),~{}g(x):=\nabla\phi(x)+n_{g}(x).$ (25)

Here, $n_{f}$ is sampled randomly from the uniform distribution over the
interval $[-e_{f},e_{f}]$, and $n_{g}$ is sampled uniformly from the
$n$-dimensional sphere with radius $e_{g}$. We compare the performance of soft
QN and SP-BFGS for solving a set of problems from the CUTEst suite [14] The
names and the dimensions of the problems are listed in Table 1 in Appendix A.
The updates of soft QN and SP-BFGS were combined with a backtracking line-
search procedure similar to the one described in [16], summarized in Algorithm
2.

Algorithm 2 Line-search procedure

1:procedure line-search under noise
($p_{k},x_{k},g(x_{k}),\varepsilon_{tol},c,\tau,T,f(\cdot)$)

2: $\eta_{k}\leftarrow 1$

3: $t\leftarrow 0$

4: while
$f(x_{k}+\eta_{k}p_{k})>f(x_{k})+\eta_{k}cp_{k}^{T}g(x_{k})+2\varepsilon_{tol}\And
t<T$ do

5: $\eta_{k}\leftarrow\tau\eta_{k}$

6: $t\leftarrow t+1$

7: end while

8: if $f(x_{k}+\eta_{k}p_{k})<f(x_{k})+2\varepsilon_{tol}$ then

9: return $\eta_{k}$

10: else

11: return 0

12: end if

13:end procedure

To scale the noise according to the problem, we set
$e_{f}=10^{-4}|\phi(x_{0})|$ and $e_{g}=10^{-4}\|\nabla\phi(x_{0})\|_{2}$,
where $x_{0}$ is the initial point. We set $\tau=0.5$, $c=10^{-4}$, $T=45$,
and $\varepsilon_{tol}=e_{f}$. As for the values of $\alpha_{k}$ for soft QN
and $\beta_{k}$ for SP-BFGS, we utilize different strategies. After some test
runs on the problem ARWHEAD, we determined that $\alpha_{k}=10^{6}$ for all
the test problems. In contrast to the experiments with quadratic problems, the
noise is low compared to the norm of the initial gradient, allowing a
significantly larger $\alpha_{k}$ without underfitting curvature information.
For $\beta_{k}$, we set
$\beta_{k}=\frac{10^{8}}{\epsilon_{g}}\|s_{k}\|_{2}+10^{-10}$, which has been
heuristically discovered by [16] to work well in practice for various CUTEst
problems. However, if this choice would result in a non-positive definite
matrix, i.e., $s_{k}^{T}y_{k}\leq-\frac{1}{\beta_{k}}$, we instead refrain
from updating the Hessian inverse approximation for that iteration.

For each test problem, we run the algorithms 30 times and track the measure

$\Delta^{i}_{j}:=\phi(x^{i}_{j})-\phi^{\star},$

where $x^{i}_{j}$ is the variable vector after $j$ function evaluations for
test-run $i$, and $\phi^{\star}$ is the optimal value that is available in the
file defining the CUTEst problem. The initial point from which the iterations
are started is also found in the same file. Note that while we use the exact
function values for the purpose of tracking suboptimality to evaluate
performance, the algorithms in the test are only given access to noisy
function values. The line-search that we apply involves a tolerance term for
the noise and does not guarantee descent in each step, so the suboptimality
may actually increase from one instance to another. Nevertheless, we observe a
clear downward trend in suboptimality in the test runs; see Figs 4(a) and
4(b), which show the performance for soft QN and SP-BFGS on the test problem
DIXMAANA.

In Figures 4(a) and 4(b), we show the median of all $\Delta^{i}_{j},\
i\in\\{1,\dots,30\\}$ against the number of function evaluations $j$, and we
show the range of all suboptimalities as well as the middle quartiles. The
reason for using this way of visualizing the results, in contrast to previous
plots, is the varying number of function evaluations in each iteration of the
algorithm. For a specific number $j$ of function evaluations, the algorithm in
test run $i$ might still be in the process of performing a line-search. The
value of $x^{i}_{j}$ is then the last iterate before $j$ function evaluations
has been performed, at which the algorithm had finished a line-search. Since
the number of function evaluations performed before computing $x^{i}_{j}$ can
vary considerably between test runs, one can expect the distribution of
$\Delta^{i}_{j}$ to be skewed, making the median a better representation of
the typical suboptimality than the mean in this case. Showing the quartiles as
well as the full range of the distribution of $\Delta^{i}_{j}$ gives an idea
of the skewness of the distribution.

For the test perfomed on the DIXMAANA problem, we can see that soft QN has a
consistently better performance than SP-BFGS; Its median suboptimality
decreases faster, the variability of its suboptimality is lower, and at the
final iteration, the highest suboptimality for all of the test runs for soft
QN is lower than the lowest suboptimality for SP-BFGS.

(a) soft QN

(b) SP-BFGS

Figure 4:  Comparision of soft QN and SP-BFGS on the CUTEst-problem DIXMAANA.
Interval of suboptimality for different test runs, plotted on a logarithmic
scale against number of function evaluations.

For all the tests with the CUTEst test problems, we terminate the algorithms
as soon as 2000 function evaluations have been exceeded, interrupting any
line-searches and determining whether the previous iterate or the current
trial iterate in the line-search should be considered the final one, in the
same way as if the usual maximum number of function evaluations in the line
search $T$ has been exceeded. We then compute summary statistics for of
$\Delta^{i}_{2000}$ for all the test runs. The results are given in Tables 2
and 3 in Appendix A.

##  6 Conclusion

We have proposed a novel quasi-Newton method to address the challenges of
minimizing loss functions in the presence of general bounded noise. By
replacing the traditionally hard secant condition with a penalty on its
violation, we developed a recursive Hessian inverse approximation that is
guaranteed to be positive definite for all penalty parameters, even for non-
convex loss functions. This soft QN update exhibits several attractive
properties: it recovers BFGS under specific limits, treats positive and
negative curvature equally, and is scale invariant. These features
collectively enhance the efficacy of soft QN in noisy environments. For
strongly convex loss functions, we have detailed an algorithm that has linear
convergence to a neighborhood of the optimal solution. Numerical experiments
demonstrate the superiority of soft QN compared with existing algorithms,
validating its practical effectiveness and highlighting its potential for
various optimization tasks in noisy scenarios.

Although this paper has provided a relatively complete treatment of the soft
QN update, several problems warrant further attention. These include exploring
ways to adapt the penalty parameter dynamically during iterations to help
escape saddle points, as well as developing statistically efficient techniques
in the presence of random function and gradient noise. Furthermore, designing
a limited memory version of soft QN would be crucial for extending its
applicability to large-scale settings. We hope to return to these issues in
our future work.

## Disclosure statement

No potential conflict of interest was reported by the authors.

## Funding

This work was supported in part by Swedish Research Council (VR) under grant
2019-05319, by the Knut and Alice Wallenberg Foundation, and Digital Futures.

## Notes on contributor(s)

_Erik Berglund_ received the B.Sc. degree in Engineering Physics and the M.Sc.
degree in Applied and Computational Mathematics from KTH Royal Institute of
Technology, Stockholm, Sweden, in 2015 and 2017 respectively. He is currently
a PhD student in the Division of Decision and Control Systems at KTH Royal
Institute of Technology. His research is focused on the use of Hessians and
Hessian approximations in optimization algorithms.

_Jiaojiao Zhang_ received the B.E. degree in automation from School of
Automation, Harbin Engineering University, Harbin, China, in 2015. She
received the master degree in control theory and control engineering from
University of Science and Technology of China, Hefei, China, in 2018. She
received the Hong Kong PhD Fellowship Scheme (HKPFS) in August 2018. She
received Ph.D. degree in the Department of Systems Engineering and Engineering
Management, Chinese University of Hong Kong, in 2022. She is currently a
postdoctoral researcher in the Division of Decision and Control Systems at KTH
Royal Institute of Technology, Stockholm, Sweden. Her current research
interests include distributed optimization and algorithm design.

_Mikael Johansson_ received the M.Sc. and Ph.D. degrees in electrical
engineering from Lund University, Lund, Sweden, in 1994 and 1999,
respectively. He held postdoctoral positions with Stanford University,
Stanford, CA, USA, and the University of California at Berkeley, Berkeley, CA,
USA, before joining the KTH Royal Institute of Technology, Stockholm, Sweden,
in 2002, where he is currently a Full Professor. He has played a leading role
in several national and international research projects on optimization,
control, and communications. He has coauthored two books and over 200 papers,
several of which are highly cited and have received recognition in terms of
best paper awards. His research interests revolve around large-scale and
distributed optimization, autonomous decision-making, control, and machine
learning. Dr. Johansson has served on the Editorial Boards of Automatica and
the IEEE Transactions on Control of Network Systems and on the program
committee for several top conferences organized by IEEE and ACM.

## References

  * [1] Agarwal, A., Bartlett, P. L., Ravikumar, P., and Wainwright, M. J. Information-theoretic lower bounds on the oracle complexity of stochastic convex optimization.  IEEE Transactions on Information Theory 58, 5 (2012), 3235–3249.
  * [2] Amaran, S., Sahinidis, N. V., Sharda, B., and Bury, S. J. Simulation optimization: a review of algorithms and applications.  Annals of Operations Research 240 (2016), 351–380.
  * [3] Berahas, A. S., Byrd, R. H., and Nocedal, J. Derivative-free optimization of noisy functions via quasi-Newton methods.  SIAM Journal on Optimization 29, 2 (2019), 965–993.
  * [4] Birge, J. R., and Louveaux, F. Introduction to stochastic programming.  Springer Science & Business Media, 2011.
  * [5] Bordes, A., Bottou, L., and Gallinari, P. SGD-QN: Careful quasi-Newton stochastic gradient descent.  Journal of Machine Learning Research 10 (2009), 1737–1754.
  * [6] Bottou, L. Large-scale machine learning with stochastic gradient descent.  177–186.
  * [7] Bottou, L., Curtis, F. E., and Nocedal, J. Optimization methods for large-scale machine learning.  SIAM review 60, 2 (2018), 223–311.
  * [8] Byrd, R. H., Hansen, S. L., Nocedal, J., and Singer, Y. A stochastic quasi-Newton method for large-scale optimization.  SIAM Journal on Optimization 26, 2 (2016), 1008–1031.
  * [9] Chang, C.-C., and Lin, C.-J. Libsvm: a library for support vector machines.  ACM transactions on intelligent systems and technology (TIST) 2, 3 (2011), 1–27.
  * [10] Chen, X., and Feng, H. A modified stochastic quasi-Newton algorithm for summing functions problem in machine learning.  Journal of Applied Mathematics and Computing 69, 2 (2023), 1491–1506.
  * [11] Dantzig, G. B. Linear programming under uncertainty.  Management science 1, 3-4 (1955), 197–206.
  * [12] Dauphin, Y. N., Pascanu, R., Gulcehre, C., Cho, K., Ganguli, S., and Bengio, Y. Identifying and attacking the saddle point problem in high-dimensional non-convex optimization.  Advances in neural information processing systems 27 (2014).
  * [13] Fletcher, R. A new variational result for quasi-Newton formulae.  SIAM Journal on Optimization 1, 1 (1991), 18–21.
  * [14] Gould, N. I., Orban, D., and Toint, P. L. CUTEst: a constrained and unconstrained testing environment with safe threads for mathematical optimization.  Computational optimization and applications 60 (2015), 545–557.
  * [15] Guler, O., Gürtuna, F., and Shevchenko, O. Variational problems in quasi–Newton methods.
  * [16] Irwin, B., and Haber, E. Secant penalized BFGS: a noise robust quasi-Newton method via penalizing the secant condition.  Computational Optimization and Applications 84, 3 (2023), 651–702.
  * [17] Jordan, M. I., and Mitchell, T. M. Machine learning: Trends, perspectives, and prospects.  Science 349, 6245 (2015), 255–260.
  * [18] Lahoti, A., Senapati, S., Rajawat, K., and Koppel, A. Sharpened lazy incremental quasi-Newton method.  arXiv preprint arXiv:2305.17283 (2023).
  * [19] Lenard, M. L., and Minkoff, M. Randomly generated test problems for positive definite quadratic programming.  ACM Trans. Math. Softw. 10, 1 (jan 1984), 86–96.
  * [20] Liu, D. C., and Nocedal, J. On the limited memory BFGS method for large scale optimization.  Mathematical programming 45, 1-3 (1989), 503–528.
  * [21] Merikoski, J. K., and Virtanen, A. Bounds for eigenvalues using the trace and determinant.  Linear algebra and its applications 264 (1997), 101–108.
  * [22] Mokhtari, A., Eisen, M., and Ribeiro, A. IQN: An incremental quasi-Newton method with local superlinear convergence rate.  SIAM Journal on Optimization 28, 2 (2018), 1670–1698.
  * [23] Mokhtari, A., and Ribeiro, A. RES: Regularized stochastic BFGS algorithm.  IEEE Transactions on Signal Processing 62, 23 (2014), 6089–6104.
  * [24] Moritz, P., Nishihara, R., and Jordan, M. A linearly-convergent stochastic L-BFGS algorithm.  249–258.
  * [25] Nocedal, J., and Wright, S. J. Numerical optimization, 2e ed.  Springer, New York, NY, USA, 2006.
  * [26] Nordström, K. Convexity of the inverse and Moore–Penrose inverse.  Linear algebra and its applications 434, 6 (2011), 1489–1512.
  * [27] Paternain, S., Mokhtari, A., and Ribeiro, A. A Newton-based method for nonconvex optimization with fast evasion of saddle points.  SIAM Journal on Optimization 29, 1 (2019), 343–368.
  * [28] Polyak, B. T. Introduction to optimization.  New York, Optimization Software, 1987.
  * [29] Robbins, H., and Monro, S. A stochastic approximation method.  The annals of mathematical statistics (1951), 400–407.
  * [30] Shi, H.-J. M., Xie, Y., Byrd, R., and Nocedal, J. A noise-tolerant quasi-Newton algorithm for unconstrained optimization.  SIAM Journal on Optimization 32, 1 (2022), 29–55.
  * [31] Sra, S., Nowozin, S., and Wright, S. J. Optimization for machine learning.  Mit Press, 2012.
  * [32] Wang, X., Ma, S., Goldfarb, D., and Liu, W. Stochastic quasi-Newton methods for nonconvex stochastic optimization.  SIAM Journal on Optimization 27, 2 (2017), 927–956.
  * [33] Wolkowicz, H., and Styan, G. P. Bounds for eigenvalues using traces.  Linear algebra and its applications 29 (1980), 471–506.
  * [34] Xie, Y., Byrd, R. H., and Nocedal, J. Analysis of the BFGS method with errors.  SIAM Journal on Optimization 30, 1 (2020), 182–209.
  * [35] Xu, P., Roosta, F., and Mahoney, M. W. Second-order optimization for non-convex machine learning: An empirical study.  199–207.
  * [36] Yasuda, S., Indrapriyadarsini, S., Ninomiya, H., Kamio, T., and Asai, H. addHessian: Combining quasi-Newton method with first-order method for neural network training.  Nonlinear Theory and Its Applications, IEICE 13, 2 (2022), 361–366.
  * [37] Zhang, T. Solving large scale linear prediction problems using stochastic gradient descent algorithms.  In Proceedings of the twenty-first international conference on Machine learning (2004), p. 116.

## Appendix A

The problem names and dimensions of the considered test problems in the CUTEst
suite are given in Table 1. Results for all the CUTEst test problems are
summarized in Tables 2 and 3. In these tables, we compare soft QN and SP-BFGS
by various metrics: The minimum, maximum, mean, median, and sample variance of
all $\Delta^{i}_{2000}$ for each problem. In each table, we have marked each
cell with a lower value than the corresponding cell in the table for the other
algorithm in green. From these tables, one can see that out of the 32 test
problems, soft QN had a lower min $\Delta^{i}_{2000}$ on 20, a lower max
$\Delta^{i}_{2000}$ on 20, a lower mean $\Delta^{i}_{2000}$ on 19, a lower
median $\Delta^{i}_{2000}$ on 22 and a lower variance of $\Delta^{i}_{2000}$
on 20 of the problems. This indicates that soft QN has a competitive
performance compared to SP-BFGS in this setting.

Problem name | $n$ | Problem name | $n$ | Problem name | $n$
---|---|---|---|---|---
ARWHEAD | 100 | DIXMAANI | 90 | GENROSE | 100
BDQRTIC | 100 | DIXMAANJ | 90 | MOREBV | 100
CRAGGLVY | 100 | DIXMAANK | 90 | NONDIA | 100
DIXMAANA | 90 | DIXMAANL | 90 | NONDQUAR | 100
DIXMAANB | 90 | DIXMAANM | 90 | QUARTC | 100
DIXMAANC | 90 | DIXMAANN | 90 | SPARSQUR | 100
DIXMAAND | 90 | DIXMAANO | 90 | TQUARTIC | 100
DIXMAANE | 90 | DIXMAANP | 90 | TRIDIA | 100
DIXMAANF | 90 | EIGENALS | 110 | WATSON | 31
DIXMAANG | 90 | EIGENBLS | 110 | WOODS | 100
DIXMAANH | 90 | EIGENCLS | 30 |  |
Table 1: The name and dimension of test problems from the CUTEst suite. Problem name | Min $\Delta_{2000}^{i}$ | Max $\Delta_{2000}^{i}$ | Mean $\Delta_{2000}^{i}$ | Median $\Delta_{2000}^{i}$ | Var $\Delta_{2000}^{i}$
---|---|---|---|---|---
ARWHEAD | 9.22E-07 | 1.52E-04 | 7.45E-06 | 2.17E-06 | 7.19E-10
BDQRTIC | 5.76E-04 | 2.07 | 1.07E-01 | 1.36E-03 | 1.15E-01
CRAGGLVY | 2.67E-03 | 8.40E+01 | 3.12 | 1.16E-02 | 2.27E+02
DIXMAANA | 1.47E-07 | 2.84E-07 | 2.23E-07 | 2.28E-07 | 1.62E-15
DIXMAANB | 5.81E-07 | 9.25E-06 | 1.57E-06 | 1.11E-06 | 2.44E-12
DIXMAANC | 3.42E-06 | 2.05E-04 | 1.83E-05 | 6.35E-06 | 1.40E-09
DIXMAAND | 1.45E-05 | 1.15E-02 | 4.98E-04 | 3.11E-05 | 4.25E-06
DIXMAANE | 3.37E-07 | 8.90E-07 | 5.48E-07 | 5.36E-07 | 2.05E-14
DIXMAANF | 1.46E-06 | 1.20E-05 | 3.34E-06 | 2.39E-06 | 5.59E-12
DIXMAANG | 2.25E-04 | 1.34E-03 | 8.15E-04 | 8.30E-04 | 8.37E-08
DIXMAANH | 3.91E-05 | 1.41E-03 | 1.69E-04 | 7.50E-05 | 6.95E-08
DIXMAANI | 1.37E-04 | 2.54E-04 | 1.95E-04 | 1.97E-04 | 9.75E-10
DIXMAANJ | 2.36E-03 | 9.03E-03 | 5.53E-03 | 5.46E-03 | 3.96E-06
DIXMAANK | 5.61E-04 | 2.96E-03 | 1.73E-03 | 1.76E-03 | 4.47E-07
DIXMAANL | 2.03E-03 | 2.21E-02 | 6.45E-03 | 5.67E-03 | 1.26E-05
DIXMAANM | 1.42E-04 | 2.20E-04 | 1.85E-04 | 1.89E-04 | 4.76E-10
DIXMAANN | 1.84E-04 | 2.65E-04 | 2.20E-04 | 2.19E-04 | 3.07E-10
DIXMAANO | 1.64E-04 | 2.16E-04 | 1.87E-04 | 1.87E-04 | 1.71E-10
DIXMAANP | 2.44E-04 | 2.40E-03 | 5.38E-04 | 4.81E-04 | 1.32E-07
EIGENALS | 3.69E-06 | 5.88E-05 | 2.25E-05 | 1.93E-05 | 2.77E-10
EIGENBLS | 2.81E-08 | 8.24E-07 | 1.58E-07 | 9.19E-08 | 3.02E-14
EIGENCLS | 1.47E-09 | 5.62E-09 | 3.83E-09 | 3.93E-09 | 9.65E-19
GENROSE | 6.80E-09 | 2.71E-08 | 1.61E-08 | 1.52E-08 | 2.58E-17
MOREBV | 8.55E-07 | 8.56E-07 | 8.55E-07 | 8.55E-07 | 1.08E-20
NONDIA | 1.78E-04 | 2.57E-02 | 1.28E-03 | 3.85E-04 | 2.06E-05
NONDQUAR | 4.42E-04 | 6.83E-04 | 5.76E-04 | 5.78E-04 | 4.02E-09
QUARTC | 1.27E+02 | 2.87E+05 | 1.64E+04 | 2.16E+02 | 3.12E+09
SPARSQUR | 5.00E-05 | 3.05E-02 | 1.59E-03 | 1.10E-03 | 3.18E-05
TQUARTIC | 1.74E-10 | 1.98E-09 | 1.06E-09 | 1.06E-09 | 2.11E-19
TRIDIA | 2.53E-07 | 1.49E-06 | 4.66E-07 | 4.60E-07 | 5.46E-14
WATSON | 2.39E-05 | 3.17E-04 | 1.29E-04 | 1.16E-04 | 5.51E-09
WOODS | 4.23E-02 | 6.30E-01 | 1.17E-01 | 1.00E-01 | 1.02E-02
Table 2:  Results from the test runs with soft QN on a set of CUTEst problems. Instances where soft QN performed better than SP-BFGS are marked in green. Problem name | Min $\Delta_{2000}^{i}$ | Max $\Delta_{2000}^{i}$ | Mean $\Delta_{2000}^{i}$ | Median $\Delta_{2000}^{i}$ | Var$\Delta_{2000}^{i}$
---|---|---|---|---|---
ARWHEAD | 2.05E-06 | 5.59E-05 | 9.43E-06 | 5.11E-06 | 1.24E-10
BDQRTIC | 5.53E-04 | 5.86E-02 | 4.61E-03 | 1.48E-03 | 1.11E-04
CRAGGLVY | 3.35E-03 | 8.40E+01 | 3.12 | 1.29E-02 | 2.27E+02
DIXMAANA | 3.86E-07 | 3.81E-04 | 1.47E-05 | 1.20E-06 | 4.65E-09
DIXMAANB | 9.57E-07 | 3.12E-05 | 4.54E-06 | 2.50E-06 | 3.11E-11
DIXMAANC | 2.97E-06 | 6.92E-04 | 4.27E-05 | 9.50E-06 | 1.58E-08
DIXMAAND | 1.55E-05 | 3.00E-02 | 1.13E-03 | 3.42E-05 | 2.89E-05
DIXMAANE | 7.61E-07 | 5.26E-05 | 4.61E-06 | 1.96E-06 | 9.02E-11
DIXMAANF | 1.92E-06 | 9.11E-05 | 1.06E-05 | 4.84E-06 | 3.06E-10
DIXMAANG | 2.90E-04 | 1.11E-01 | 6.03E-03 | 1.09E-03 | 4.11E-04
DIXMAANH | 3.60E-05 | 1.19E-01 | 4.39E-03 | 7.41E-05 | 4.55E-04
DIXMAANI | 1.38E-04 | 2.79E-04 | 2.05E-04 | 2.08E-04 | 1.38E-09
DIXMAANJ | 2.35E-03 | 8.98E-03 | 5.27E-03 | 4.88E-03 | 3.78E-06
DIXMAANK | 6.27E-04 | 3.14E-03 | 1.73E-03 | 1.80E-03 | 5.10E-07
DIXMAANL | 2.14E-03 | 1.06E-02 | 5.72E-03 | 5.57E-03 | 2.67E-06
DIXMAANM | 1.23E-04 | 2.36E-04 | 1.90E-04 | 1.93E-04 | 7.49E-10
DIXMAANN | 1.88E-04 | 6.02E-03 | 5.11E-04 | 2.23E-04 | 1.28E-06
DIXMAANO | 1.58E-04 | 2.67E-04 | 1.97E-04 | 1.94E-04 | 3.85E-10
DIXMAANP | 2.29E-04 | 6.82E-04 | 4.74E-04 | 4.92E-04 | 9.86E-09
EIGENALS | 2.33E-06 | 7.63E-05 | 2.46E-05 | 1.86E-05 | 3.46E-10
EIGENBLS | 4.65E-08 | 2.77E-07 | 1.31E-07 | 1.05E-07 | 3.18E-15
EIGENCLS | 4.22E-09 | 8.63E-08 | 1.83E-08 | 1.48E-08 | 2.25E-16
GENROSE | 8.39E-09 | 1.15E-07 | 5.42E-08 | 5.30E-08 | 5.19E-16
MOREBV | 4.78E-12 | 3.86E-09 | 2.88E-10 | 9.31E-11 | 4.70E-19
NONDIA | 2.05E-04 | 6.26E-03 | 6.99E-04 | 4.38E-04 | 1.12E-06
NONDQUAR | 4.40E-04 | 2.37E-03 | 6.40E-04 | 5.68E-04 | 1.15E-07
QUARTC | 1.31E+02 | 4.96E+04 | 2.39E+03 | 2.43E+02 | 8.24E+07
SPARSQUR | 5.76E-05 | 2.05E-03 | 2.02E-04 | 9.80E-05 | 1.28E-07
TQUARTIC | 1.06E-10 | 1.99E-09 | 7.38E-10 | 7.00E-10 | 1.42E-19
TRIDIA | 1.23E-06 | 6.28E-06 | 3.05E-06 | 2.76E-06 | 1.55E-12
WATSON | 2.42E-05 | 3.46E-04 | 1.32E-04 | 1.12E-04 | 7.26E-09
WOODS | 4.08E-02 | 1.58E-01 | 9.19E-02 | 8.98E-02 | 8.00E-04
Table 3:  Results from the test runs with SP-BFGS on a set of CUTEst problems.
Instances where SP-BFGS performed better than soft QN are marked in green.