Approximating Hessian matrices using Bayesian inference: a new approach for quasi-Newton methods in stochastic optimization

Abstract

Using quasi-Newton methods in stochastic optimization is not a trivial task given the difficulty of extracting curvature information from the noisy gradients. Moreover, pre-conditioning noisy gradient observations tend to amplify the noise. We propose a Bayesian approach to obtain a Hessian matrix approximation for stochastic optimization that minimizes the secant equations residue while retaining the extreme eigenvalues between a specified range. Thus, the proposed approach assists stochastic gradient descent to converge to local minima without augmenting gradient noise. We propose maximizing the log posterior using the Newton-CG method. Numerical results on a stochastic quadratic function and an ${\ell }_2$-regularized logistic regression problem are presented. In all the cases tested, our approach improves the convergence of stochastic gradient descent, compensating for the overhead of solving the log posterior maximization. In particular, pre-conditioning the stochastic gradient with the inverse of our Hessian approximation becomes more advantageous the larger the condition number of the problem is.

Keywords:

• Stochastic optimization
• quasi-Newton
• Monte Carlo
• variance reduction
• control variates
• machine learning

AMS subject classifications:

• 65K10
• 90C15
• 62F15
• 90C53

1. Introduction

First-order optimization methods are a common choice for performing local search when gradients are available. In the strongly convex case, the gradient descent method converges linearly, requiring a number of iterations proportional to the condition number [21, Theorem 2.1.15], that is, the ratio between the upper and lower bounds of the eigenvalues of the Hessian matrix in the search domain. However, first-order methods may not be a practical choice when the condition number for a problem is large, advancing little towards the optimum each iteration. One way to mitigate this dependency is to use quasi-Newton methods, which use the secant equation to obtain an approximation of the Hessian matrix of the objective function from subsequent gradient evaluations. Then, one may precondition the gradient with the inverse of the Hessian approximation, furnishing a Newton-like method without the cost of computing the true Hessian or using a finite differences type of approximation [10].

In the case of stochastic optimization, the stochastic gradient descent method (SGD) is the canonical first-order method. The asymptotic convergence of SGD is dominated by noise in gradient evaluations (cf. [24, Theorem 3]), thus controlling the condition number for the problem does not improve convergence. However, controlling the relative statistical error of gradient estimates using a uniform bound recovers linear convergence per iteration for the minimization of strongly convex functions [8,24]. In this case, the convergence rate depends on the condition number of the problem; hence preconditioning the gradient with the Hessian inverse might be desirable. The natural approach is to use quasi-Newton methods to approximate the Hessian matrix, but devising quasi-Newton methods for stochastic optimization is challenging due to its noisy gradient observations. Naively fitting the secant equations on noisy gradient observations often leads to catastrophic errors in Hessian approximations [24, Section 4.3.3].

There are many works addressing the development of quasi-Newton methods in stochastic optimization. One of the first attempts is that of Schraudolph, Yu, and Günther [26] where they adapt the L-BFGS method to the stochastic setting by adding damping parameters that control the effect of noise, furnishing the oL-BFGS method. Their approach is empirical, and not much can be said about the quality of the Hessian approximation. Bordes and Pierre [3] propose a diagonal rescaling matrix for linear support vector machines based on oL-BFGS, providing an analysis of their approach. In their numerical tests, their method matched oL-BFGS in convergence per epoch but with a significant improvement in terms of runtimes. Hennig [15] uses a Bayesian approach to learn about an objective function from noisy gradient evaluations. Byrd et al. [6] use a sample average approximation to compute the search direction in Newton-CG and L-BFGS. As a consequence, they do not have a statistical error on their estimates but a bias. Hennig [15] models the prior distribution of the objective function as a Gaussian process with a squared exponential kernel and, using the linearity of Gaussian processes, derives prior distributions for the gradient and the Hessian of the objective function. Then, assuming Gaussian likelihoods for both the gradient and the Hessian of the objective function, Hennig [15] deduces closed-form equations for the mean of the (also Gaussian) posterior distribution. Sohl-Dickstein, Poole, and Ganguli [28] develop a quasi-Newton method for sums of functions where they keep a quadratic approximation for each of these functions, subsampling them at each iteration and updating the respective quadratic approximation with new gradient information. They use a shared, low-dimensional subspace to avoid the cost of storing a Hessian matrix for each function. Byrd et al. [7] propose a stochastic quasi-Newton method comprising a two-loop scheme where the outer loop updates the curvature pairs in an L-BFGS scheme (with large sample sizes), and the inner loop performs preconditioned stochastic gradient descent (with lower sample sizes). Thus, the method proposed by Byrd et al. [7] does not incorporate the information from the noisy gradients of the inner loops in the Hessian approximation. Moritz, Nishihara, and Jordan [20] improve the stochastic quasi-Newton method proposed by Byrd et al. [7] by combining it with a control variate method to estimate the gradient at each iteration, achieving linear convergence.

Gower, Goldfarb, and Richtárik [13] introduce a block BFGS method to update the Hessian approximation each iteration using a sample independent of the sample used to compute the gradient. The proposed block update fits the Hessian inverse to a sketched Hessian computed from a sub-sample of the data, but requires a Cholesky decomposition each iteration. Wang et al. [30] propose a damped quasi-Newton method for stochastic optimization that uses all noisy gradient evaluations to approximate the Hessian inverse. They provide a convergence proof for nonconvex problems, provided the Hessian eigenvalues have uniform lower bound. Li [18] preconditions the stochastic gradient estimates using a matrix computed to fit the secant equation while avoiding amplifying gradient estimate noise. Bollapragada et al. [2] use standard L-BFGS equations in stochastic optimization by increasing sample sizes to keep the relative error for the gradients uniformly bounded. Bollapragada et al. [2] recommend a stochastic line-search procedure to be used with their method. Goldfarb, Ren, and Bahamou [12] develop a quasi-Newton method for deep neural networks that approximates Hessian inverses as block matrices such that each block is decomposed as the kronecker products of smaller matrices. Lately, Shi et al. [27] build on top of their analysis in [31] to build a quasi-Newton method for stochastic optimization with a noise-resilient Wolfe–Armijo line search. Their method, however, requires the gradient estimate to be bounded above in norm by a known constant, which they use in their line search. We refer the interested reader to [4,19] for extensive reviews of quasi-Newton methods for stochastic optimization.

Most approaches to tackle the condition number issue in stochastic optimization use equations from deterministic quasi-Newton methods and somehow circumvent the randomness. Here, we propose a Bayesian framework to build a posterior distribution for the Hessian matrix from noisy gradient observations. In Bayesian inference, a posterior distribution is built from prior and likelihood distributions. We model the prior distribution for the Hessian matrix with a probability density function (pdf) that decreases exponentially as the Hessian moves farther away (in the Frobenius norm sense) from the initial Hessian estimate. This approach is motivated by classic deterministic quasi-Newton methods, where a Frobenius norm regularizer centred at the current Hessian is used to circumvent ill-posedness [10]. Moreover, we impose constraints on the Hessian matrix eigenvalues extremes, enforcing positive-definiteness and gradient smoothness. We model the likelihood using the secant equation, i.e. how likely it is to obtain the observed gradient differences for a given Hessian. The Hessian that maximizes the posterior distribution, i.e. the Maximum a Posteriori (MAP), is defined as the Hessian approximation.

The choice of prior and likelihood distributions guarantees that the negative log posterior is strongly-convex. Therefore, we opt to use the Newton-CG method to find the MAP. Although it is possible to compute the second-order derivatives of the negative log posterior of the Hessian directly, it could be expensive and memory intensive. Keep in mind that, for a d-dimensional problem, the Hessian matrix is a $d\times d$ matrix, thus, the second-order derivatives of the log posterior must be gathered in a fourth-order tensor with $d^4$ components. Instead of assembling this tensor, we compute its action over any symmetric matrix of adequate dimensions and then use the conjugate gradient method to find the Newton direction. Moreover, to enforce the stability of Newton-CG, we use the central-path interior point method to impose the prior conditions on the eigenvalues of the Hessian approximation, i.e. we start with a large penalization parameter and decrease it in steps, keeping the problem well-conditioned. If properly tuned, our approach not only guarantees that Newton-CG converges to the MAP as it also results in a slow decay for the Hessian approximation's smallest eigenvalue, thus improving the stability for the SGD preconditioned with our Hessian approximation.

Other works address the problem of approximating a Hessian from noisy gradient observations. Hennig [15] also uses Bayesian inference; however, to build a posterior distribution of the objective function as a Gaussian process. The main difference between the proposed approach and that of Hennig [15] is in the way he models the prior and likelihood distributions. Moreover, Hennig [15] assumes independence between the components of the gradient noise in his analysis, and also that they have the same variance. In the same spirit of the present work, Li [18] also solves an optimization sub-problem each iteration to find a Hessian approximation from noisy gradients. His approach differs from ours in the way he uses the available information to estimate the preconditioning matrix; here, we use a Bayesian approach to build a posterior distribution of the Hessian matrix while [18] obtains a point estimate by minimizing a different objective function.

To illustrate the effectiveness of our approach, we solve two stochastic optimization problems. The first is minimizing a quadratic function, and the second is training an ${\ell }_2$ regularized logistic regression model. In all cases, we start with SGD preconditioned with a diagonal matrix and update the preconditioning matrix in intervals using the proposed Bayesian approach. We study cases with up to 300 design variables and condition numbers of orders up to ${10}^9$. We present results combining the proposed Bayesian Hessian approach with different variations of SGD, namely, vanilla SGD, SGD-MICE [8], SVRG [16], and SARAH [22].

1.1. Problem definition

The stochastic optimization problem is stated as follows. Let $\mathit{\xi }$ be the design variable in dimension $d_{\mathit{\xi }}$ and $\mathit{\theta }$ be a vector-valued random variable in dimension $d_{\mathit{\theta }}$ equipped with a probability density function $\varrho :{\mathbb{R}}^{d_{\mathit{\theta }}}\to {\mathbb{R}}^{+}$. Here, $\mathbb{E}[\cdot |\mathit{\xi }]$ and $\mathbb{V}[\cdot |\mathit{\xi }]$ are the expectation and variance operators conditioned on $\mathit{\xi }$, respectively. Aiming at minimizing the expectations conditioned on $\mathit{\xi }$, we state the problem as (1) ${\mathit{\xi }}^{\ast }=\underset{\mathit{\xi }\in \mathrm{\Xi }}{\arg min}\mathbb{E}\left[f\right(\mathit{\xi },\mathit{\theta }\left)|\mathit{\xi }\right],$(1) where $\mathrm{\Xi }\subset {\mathbb{R}}^{d_{\mathit{\xi }}}$ is a feasible set and $f:\mathrm{\Xi }\times {\mathbb{R}}^{d_{\mathit{\theta }}}\to \mathbb{R}$. Through what follows, let $F\left(\mathit{\xi }\right):=\mathbb{E}\left[f\right(\mathit{\xi },\mathit{\theta }\left)|\mathit{\xi }\right]$ where $F:\mathrm{\Xi }\to \mathbb{R}$.

In minimizing (1(1) ${\mathit{\xi }}^{\ast }=\underset{\mathit{\xi }\in \mathrm{\Xi }}{\arg min}\mathbb{E}\left[f\right(\mathit{\xi },\mathit{\theta }\left)|\mathit{\xi }\right],$(1) ) with respect to the design variable $\mathit{\xi }\in \mathrm{\Xi }$, vanilla SGD is constructed with the update rule (2) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\mathit{\upsilon }}_k,$(2) where ${\eta }_k$ is the step size at iteration k, and ${\mathit{\upsilon }}_k$ is an unbiased estimator of the gradient of F at ${\mathit{\xi }}_k$. For instance, an unbiased estimator ${\mathit{\upsilon }}_k$ of the gradient of F at ${\mathit{\xi }}_k$ may be constructed using a Monte Carlo estimator, (3) ${\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_k\right)=\mathbb{E}\left[f\right({\mathit{\xi }}_k,\mathit{\theta }\left)|{\mathit{\xi }}_k\right]\approx {\mathit{\upsilon }}_k:={\displaystyle \frac{1}{b}}\sum _{i=1}^b{\mathrm{\nabla }}_{\mathit{\xi }}f({\mathit{\xi }}_k,{\mathit{\theta }}_i),$(3) with b independent and identically distributed (iid) random variables ${\mathit{\theta }}_i\sim \varrho$. The estimator (3(3) ${\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_k\right)=\mathbb{E}\left[f\right({\mathit{\xi }}_k,\mathit{\theta }\left)|{\mathit{\xi }}_k\right]\approx {\mathit{\upsilon }}_k:={\displaystyle \frac{1}{b}}\sum _{i=1}^b{\mathrm{\nabla }}_{\mathit{\xi }}f({\mathit{\xi }}_k,{\mathit{\theta }}_i),$(3) ) is a random variable and its use in optimization algorithms gives rise to the so-called Stochastic Optimizers.

First-order methods are the most popular choice in stochastic optimization due to their robustness to noise. However, these methods might suffer from slow convergence when the condition number of the problem is large. Here, we are interested in finding a matrix ${\hat{\mathit{B}}}_k$ that approximates ${\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)$ so that its inverse can be used as a preconditioner for the stochastic gradient, (4) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k,$(4) thus diminishing the effect of the condition number on the optimization convergence. In the deterministic setting, quasi-Newton methods circumvent the sensibility to the condition number by using the gradient observations to construct the Hessian matrix approximation at the current iterate. However, in the stochastic optimization setting, we first have to address the following questions:

Q1:	How does one obtain a Hessian approximation from noisy gradient observations?
Q2:	How do we avoid amplifying the noise in gradient estimates with our preconditioner matrix ${\hat{\mathit{B}}}_k^{-1}$?
Q3:	To mitigate the effect of the condition number, is it enough to have an accurate Hessian approximation, however, acting on noisy gradient estimates?

To address Q1, consider the curvature pair ${\mathit{y}}_k:={\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_k\right)-{\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_{k-1}\right)$ and ${\mathit{s}}_k:={\mathit{\xi }}_k-{\mathit{\xi }}_{k-1}$. Thus, we first aim to find the Hessian approximation that minimizes the residue of the secant equation, that is, (5) $\underset{\mathit{B}\in {\mathbb{R}}^{d_{\mathit{\xi }}}\times {\mathbb{R}}^{d_{\mathit{\xi }}}}{min}\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\mathit{y}}_{\ell }{\Vert }^2.$(5) However, we do not have access to ${\mathit{y}}_k:=\mathbb{E}\left[{\mathrm{\nabla }}_{\mathit{\xi }}f\right({\mathit{\xi }}_k,\mathit{\theta }\left)|{\mathit{\xi }}_k\right]-\mathbb{E}\left[{\mathrm{\nabla }}_{\mathit{\xi }}f\right({\mathit{\xi }}_{k-1},\mathit{\theta }\left)|{\mathit{\xi }}_{k-1}\right]$, instead, we may sample finitely many times ${\hat{\mathit{y}}}_k\left(\mathit{\theta }\right):={\mathrm{\nabla }}_{\mathit{\xi }}f({\mathit{\xi }}_k,\mathit{\theta })-{\mathrm{\nabla }}_{\mathit{\xi }}f({\mathit{\xi }}_{k-1},\mathit{\theta })$. In Section 2, we discuss how to use Bayesian inference to find a matrix ${\hat{\mathit{B}}}_k$ given the noisy curvature pairs $\hat{\mathit{y}}$ and $\mathit{s}$.

One of the main problems faced by preconditioning in stochastic optimization is the amplification of the noise in gradient estimates. This amplification is proportional to the smallest eigenvalue of matrix ${\hat{\mathit{B}}}_k$, (6) $\mathbb{E}\left[{\Vert {\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k-\mathbb{E}\left[{\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k|{\hat{\mathit{B}}}_k^{-1}\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right]\le {\Vert {\hat{\mathit{B}}}_k^{-1}\Vert }^2\mathbb{E}\left[{\Vert {\mathit{\upsilon }}_k-\mathbb{E}\left[{\mathit{\upsilon }}_k\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right],$(6) by the matrix-norm compatibility. Thus, to tackle Q2, we impose a lower bound constraint on the eigenvalues of ${\hat{\mathit{B}}}_k$, i.e. ${\hat{\mathit{B}}}_k⪰\tilde{\mu }\mathit{I}$.

To address Q3 in a numerical fashion, we investigate and compare stochastic gradients preconditioned with our Hessian inverse approximation with and without control in the relative statistical error.

Regarding notation, in this work, we use bold-face lower-case letters as vectors and bold-face upper-case letters as matrices. Moreover, we refer to first-order derivatives of functions with respect to matrices as gradients, even though they are also matrices.

2. Approximating the Hessian and its inverse using Bayesian inference

Important to what follows is the Monte Carlo estimator of ${\mathit{y}}_k$ (7) ${\overline{\mathit{y}}}_k:=\frac{1}{b_k}\sum _{i=1}^{b_k}{\hat{\mathit{y}}}_k\left({\mathit{\theta }}_i\right),$(7) with $b_k$ iid random variables ${\mathit{\theta }}_i\sim \varrho$.

2.1. Bayesian formulation of the Hessian inference problem

For the pair of sets $Y:=\{{\overline{\mathit{y}}}_{\ell }{\}}_{\ell =1}^k$ and $S:=\{{\mathit{s}}_{\ell }{\}}_{\ell =1}^k$, we employ the Bayes' formula to build a posterior distribution of the Hessian $\mathit{B}$ given Y and S, (8) $\pi (\mathit{B}|S,Y)\propto p(S,Y|\mathit{B})\pi \left(\mathit{B}\right),$(8) where the likelihood distribution $p(S,Y|\mathit{B})$ and the prior distribution $\pi \left(\mathit{B}\right)$ will be discussed in Sections 2.1.1 and 2.1.2, respectively. We opt to use as a Hessian approximation the matrix that minimizes the negative log posterior subject to the eigenvalues constraints (9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) where ${\lambda }_{min}$ and ${\lambda }_{max}$ are operators that return the minimum and maximum eigenvalues, respectively; L is the smoothness constant of F; and $\tilde{\mu }>0$ is an arbitrary parameter. If F is μ-convex, $\tilde{\mu }$ can be chosen as μ, however, a larger $\tilde{\mu }$ might be desirable to avoid amplifying the noise in the search direction. Since computing ${\mathit{B}}_k^{\ast }$ is not possible in general, we use an approximation ${\hat{\mathit{B}}}_k$, obtained numerically, as a Hessian approximation.

2.1.1. Likelihood of observing curvature pairs given a Hessian

We use the secant equations to derive the likelihood of observing the curvature pairs $(S,Y)$ given a Hessian $\mathit{B}$. Let the difference in gradients be as in (7(7) ${\overline{\mathit{y}}}_k:=\frac{1}{b_k}\sum _{i=1}^{b_k}{\hat{\mathit{y}}}_k\left({\mathit{\theta }}_i\right),$(7) ). From the observations ${\hat{\mathit{y}}}_k$, we build a sample covariance matrix ${\mathbf{\Sigma }}_k$. Then, the covariance matrix of the estimator ${\overline{\mathit{y}}}_k$ is given by (10) ${\overline{\mathbf{\Sigma }}}_k:=\mathrm{Cov}\left[{\overline{\mathit{y}}}_k\right]=\frac{{\mathbf{\Sigma }}_k}{b_k}.$(10) We are interested in the inverse of ${\overline{\mathbf{\Sigma }}}_k$. However, to avoid singularity issues, we use a precision matrix approximation ${\mathit{P}}_{\ell }$ instead of ${\overline{\mathbf{\Sigma }}}_{\ell }^{-1}$ for each ${\overline{\mathit{y}}}_{\ell }\in Y$. We define ${\mathit{P}}_{\ell }$ as (11) ${\mathit{P}}_{\ell }:={({\overline{\mathbf{\Sigma }}}_{\ell }+{\sigma }_P(max\left({\overline{\mathbf{\Sigma }}}_{\ell }\right)\left){\mathit{I}}_{d_{\xi }}\right)}^{-1},$(11) where ${\sigma }_P$ is a regularization parameter and ${\mathit{I}}_{d_{\xi }}$ is the identity matrix of size $d_{\xi }$. Then, we define the negative log likelihood of curvature pair ${\mathit{s}}_k$ and ${\overline{\mathit{y}}}_k$ given a matrix $\mathit{B}$ as (12) $-\log \left(p\right({\mathit{s}}_k,{\overline{\mathit{y}}}_k|\mathit{B}\left)\right)=\frac{1}{2\nu }\Vert \mathit{B}{\mathit{s}}_k-{\overline{\mathit{y}}}_k{\Vert }_{{\mathit{P}}_k}^2+C_k,$(12) thus, for S and Y generated by a first-order stochastic gradient descent method, (13) $-\log \left(p\right(S,Y|\mathit{B}\left)\right)=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+C.$(13) In the above equations, $C=\sum _{i=\ell }^k{C}_{\ell }$ are constants and ν is used to normalize the negative log likelihood, (14) $\nu =\sum _{\ell =1}^k\Vert {\mathit{P}}_{\ell }({\hat{\mathit{B}}}_{k-1}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\Vert \Vert {\mathit{s}}_{\ell }\Vert ,$(14) being ${\hat{\mathit{B}}}_{k-1}$ a Hessian approximation at iteration k−1.

2.1.2. Prior distribution of the Hessian matrix

The prior distribution must encode the knowledge available about the Hessian before the observations of S and Y. In spirit of quasi-Newton methods, we want to discourage a Hessian that differs much from the current estimate. For this reason, we adopt a negative log prior of $\mathit{B}$ as the squared weighted Frobenius norm of the difference between $\mathit{B}$ and ${\hat{\mathit{B}}}_{k-1}$,

for a nonsingular symmetric matrix of weights $\mathit{W}$. Another information available about the matrix $\mathit{B}$ is that its eigenvalues must lay between $\tilde{\mu }$ and L, as described in the constraints of (9(9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) ). To enforce these eigenvalue constraints on $\mathit{B}$, we define the prior distribution from its negative logarithm as (16) $-\log \left(\pi \right(\mathit{B}\left)\right)=\{\begin{array}{ll}\frac{1}{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2+C& \mathrm{if}\hat{\mu }{\mathit{I}}_{d_{\xi }}⪯\mathit{B}⪯\hat{L}{\mathit{I}}_{d_{\xi }}\\ +\mathrm{\infty }& \mathrm{otherwise},\end{array}$(16) where again C is a normalizing constant. Here, we opt to impose the eigenvalues constraints on the negative log prior as logarithmic barriers of form (17) $-\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))-\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B})),$(17) where (18) $\hat{\mu }:=\frac{\tilde{\mu }}{\alpha },\hat{L}:=\mathrm{L\alpha },$(18)

and $\alpha \ge 1$ is a relaxation parameter. Hence, we propose a family of priors with logarithmic barriers with hyperparameter β as (19) $\begin{array}{rl}-\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C,\end{array}$(19) being supported on the set of symmetric positive-definite matrices with eigenvalues between $\hat{\mu }$ and $\hat{L}$.

The parameter $\rho >0$ weights the Frobenius regularization term and $\text{β}>0$ weights the logarithmic barrier terms. Note that $\log \left(\pi \right(\mathit{B};\text{β}\left)\right)\to \log \left(\pi \right(\mathit{B}\left)\right)$ as $\text{β}\downarrow 0$ with $\rho =1$.

Some options for $\mathit{W}$ are ${\mathit{I}}_{d_{\xi }}$ and . According to [10], the Frobenius norm as in (15) with weights measures the relative error of $\mathit{B}$ with respect to ${\hat{\mathit{B}}}_{k-1}$. In practice, estimating can be expensive. In these cases, we recommend using $\mathit{W}={\mathit{I}}_{d_{\xi }}$.

2.1.3. Posterior distribution of the Hessian matrix

Substituting (13(13) $-\log \left(p\right(S,Y|\mathit{B}\left)\right)=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+C.$(13) ) and (19(19) $\begin{array}{rl}-\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C,\end{array}$(19) ) in (8(8) $\pi (\mathit{B}|S,Y)\propto p(S,Y|\mathit{B})\pi \left(\mathit{B}\right),$(8) ), we obtain the posterior distribution of the Hessian matrix $\mathit{B}$ given curvature pairs S and Y. The negative log posterior with parameter β is (20) $\begin{array}{rl}{\mathcal{L}}_k(\mathit{B};\text{β})& :=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C.\end{array}$(20) The Equation (20(20) $\begin{array}{rl}{\mathcal{L}}_k(\mathit{B};\text{β})& :=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C.\end{array}$(20) ) can be interpreted as the sum of the secant equation residues for each curvature pair (in the metrics induced by ${\mathit{P}}_{\ell }$), subject to the logarithmic barrier terms that impose the eigenvalues constraints, and a Frobenius norm regularization term.

Since we are interested in minimizing the negative log posterior in (20(20) $\begin{array}{rl}{\mathcal{L}}_k(\mathit{B};\text{β})& :=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C.\end{array}$(20) ), it is our interest to know if it is strongly convex. The negative log-likelihood is convex but not strictly-convex until enough linearly independent curvature pairs $({\mathit{s}}_{\ell },{\mathit{y}}_{\ell })$ are observed. Regarding the logarithmic barrier terms of the prior, the negative logarithm of the determinant defined on symmetric positive-definite matrices is a convex function [5, Section 3.1.5]; however, not a strictly-convex one. The Frobenius norm regularization term with $\mathit{W}={\mathit{I}}_{d_{\xi }}$ is ρ-convex, therefore the negative log posterior is strongly convex with parameter ρ. Here, the regularization term works as a Tikhonov regularization term ($\Vert \mathit{A}{\Vert }_F\equiv \Vert \mathrm{vec}\left(\mathit{A}\right){\Vert }_2$), guaranteeing uniqueness of solution. Properly tuning β and ρ guarantees that the problem of minimizing (20(20) $\begin{array}{rl}{\mathcal{L}}_k(\mathit{B};\text{β})& :=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C.\end{array}$(20) ) is well-conditioned. A discussion on how to choose β is presented in Section 2.2.

2.2. Finding the posterior maximizer

Finding the ${\mathit{B}}_k^{\ast }$ that solves the problem in (9(9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) ) is not a trivial task. Since the negative log posterior is ρ-convex, local search methods like the Newton method can be used to find its minimizer. However, the Newton method requires assembling the fourth-order tensor of second-order derivatives of the negative log posterior, which can be expensive in terms of processing and memory allocation. Instead, we opt to derive the action of this fourth-order tensor of second-order derivatives over any given matrix and use the conjugate gradient method to find the Newton direction, a procedure known as the Newton-CG method. Since we have a constraint in (9(9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) ) that $\mathit{B}$ should be symmetric, the first-order optimality condition is (21) $\text{sym}\left({\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k\right(\mathit{B};\text{β}\left)\right)=0,$(21) where the $\text{sym}(\cdot )$ operator returns the symmetric part of its argument, e.g. $\text{sym}\left(\mathit{A}\right)=\frac{1}{2}(A+A^T)$. Thus, to use the Newton-CG method to solve this problem, every iteration we need to find the direction $\mathit{D}$ that satisfies (22) $\text{sym}\left({\delta }_{\mathit{D}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k\right(\mathit{B};\text{β}\left)\right)=-\text{sym}\left({\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k\right(\mathit{B};\text{β}\left)\right),$(22) where ${\delta }_{\mathit{D}}$ is the directional derivative with respect to $\mathit{B}$ in the direction of $\mathit{D}$, (23) ${\delta }_{\mathit{D}}g\left(\mathit{B}\right)=\underset{\vartheta \downarrow 0}{lim}\frac{g(\mathit{B}+\vartheta \mathit{D})-g\left(\mathit{B}\right)}{\vartheta }.$(23) To avoid ill-conditioning, we need a large enough regularization parameter ρ, but not too large since it introduces bias in the posterior maximizer. Also, the bias given by the log-barriers grows with β, the logarithmic barrier parameter; but a small β negatively impacts the smoothness of the posterior distribution. We opt to keep ρ fixed and use an interior-point central-path method [5, Section 11.2.2] to decrease the parameter β; the solution ${\hat{\mathit{B}}}_k\left(\text{β}\right)$ of the approximated problem (24) ${\hat{\mathit{B}}}_k\left(\text{β}\right)=\arg \underset{\mathit{B}}{min}\left({\mathcal{L}}_k\right(\mathit{B};\text{β}\left)\right)$(24) converges to the solution ${\mathit{B}}_k^{\ast }$ of (9(9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) ) as $\text{β}\downarrow 0$. Thus, we start with a large β and a large tolerance tol and perform optimization until the norm of the gradient of the log posterior is below tol. Then, we decrease both β and tol by a factor γ and repeat the process. This procedure is repeated until the desired tolerance is achieved with the desired β. To minimize the negative log posterior for a given β, we use the Newton-CG method with backtracking-Armijo line search; we use the conjugate gradient method to find the Newton direction and then the backtracking-Armijo to find a step size. Therefore, to use the Newton-CG, we need the gradient of ${\mathcal{L}}_k$ with respect to $\mathit{B}$ and the directional derivative of the gradient of ${\mathcal{L}}_k$ with respect to $\mathit{B}$ in the direction of any arbitrary symmetric matrix $\mathit{V}\in {\mathbb{R}}^{d_{\mathit{\xi }}}\times {\mathbb{R}}^{d_{\mathit{\xi }}}$. The pseudo-code for the minimization of the negative log posterior is presented in Algorithm 1. We use the central-path method to decrease β and tol every step outer loop, indexed as i in Algorithm 1. In lines 6 to 10, we have the Newton-CG method, where, each iteration, the conjugate gradient method is used to find the approximate Newton direction $\mathit{D}$ and a backtracking-Armijo line-search is used to find the step size.

Next, we will present how to compute ${\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k({\hat{\mathit{B}}}_k;\text{β})$ and ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k({\hat{\mathit{B}}}_k;\text{β})$ for any symmetric $\mathit{V}\in {\mathbb{R}}^{d_{\mathit{\xi }}}\times {\mathbb{R}}^{d_{\mathit{\xi }}}$. Note that, since $\mathit{B}$ is a matrix, ${\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k$ and ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k$ are also symmetric matrices of the same dimensions of $\mathit{B}$.

2.3. Gradient of the negative log posterior

To compute ${\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k$, we need to compute the gradients of the log prior and of the log likelihood, (25) ${\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=-{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right).$(25) The gradient of the negative log likelihood can be calculated from (13(13) $-\log \left(p\right(S,Y|\mathit{B}\left)\right)=\frac{1}{2\nu }\sum _{\ell =1}^k\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2+C.$(13) ) as (26) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)& =\frac{1}{2\nu }\sum _{\ell =1}^k{\mathrm{\nabla }}_{\mathit{B}}\Vert \mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell }{\Vert }_{{\mathit{P}}_{\ell }}^2\end{array}$(26) (27) $\begin{array}{rl}& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }(\mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }.\end{array}$(27) The gradient of the negative log prior can be derived from (19(19) $\begin{array}{rl}-\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\frac{\rho }{2}\Vert \mathit{B}-{\hat{\mathit{B}}}_{k-1}{\Vert }_{F,\mathit{W}}^2-\text{β}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))+C,\end{array}$(19) ) as (28) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\frac{\rho }{2}{\mathrm{\nabla }}_{\mathit{B}}\Vert {\mathit{W}}^{\frac{1}{2}}(\mathit{B}-{\hat{\mathit{B}}}_{k-1}){\mathit{W}}^{\frac{1}{2}}{\Vert }_F^2-\text{β}{\mathrm{\nabla }}_{\mathit{B}}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}{\mathrm{\nabla }}_{\mathit{B}}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B})).\end{array}$(28) The gradient of the Frobenius norm regularizer is (29) $\frac{1}{2}{\mathrm{\nabla }}_{\mathit{B}}\Vert {\mathit{W}}^{\frac{1}{2}}(\mathit{B}-{\hat{\mathit{B}}}_{k-1}){\mathit{W}}^{\frac{1}{2}}{\Vert }_F^2=\mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W},$(29) whereas the gradients of the logarithmic barrier terms are derived using the Jacobi formula, (30) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))& =-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\end{array}$(30) (31) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))& =(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(31) Substituting (29(29) $\frac{1}{2}{\mathrm{\nabla }}_{\mathit{B}}\Vert {\mathit{W}}^{\frac{1}{2}}(\mathit{B}-{\hat{\mathit{B}}}_{k-1}){\mathit{W}}^{\frac{1}{2}}{\Vert }_F^2=\mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W},$(29) ), (30(30) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))& =-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\end{array}$(30) ), and (31(31) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}))& =(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(31) ) in (28(28) $\begin{array}{rl}-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\frac{\rho }{2}{\mathrm{\nabla }}_{\mathit{B}}\Vert {\mathit{W}}^{\frac{1}{2}}(\mathit{B}-{\hat{\mathit{B}}}_{k-1}){\mathit{W}}^{\frac{1}{2}}{\Vert }_F^2-\text{β}{\mathrm{\nabla }}_{\mathit{B}}\log (det(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}))\\ & -\text{β}{\mathrm{\nabla }}_{\mathit{B}}\log (det(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B})).\end{array}$(28) ), (32) $-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)=\rho \mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}-\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.$(32) Finally, substituting the gradient of the negative log likelihood (27(27) $\begin{array}{rl}& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }(\mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }.\end{array}$(27) ) and the gradient of the negative log prior (32(32) $-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)=\rho \mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}-\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.$(32) ) in (25(25) ${\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=-{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right).$(25) ), we obtain the gradient of the negative log posterior (33) $\begin{array}{rl}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }(\mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }+\rho \mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}\\ & -\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(33)

2.4. Directional derivative of the gradient of the negative log posterior

As discussed in Section 2.2, to find the posterior distribution maximizer using Newton-CG, we need to compute the directional derivative of the gradient of the negative log posterior in a given direction. Thus, for any symmetric matrix $\mathit{V}\in {\mathbb{R}}^{d_{\mathit{\xi }}}\times {\mathbb{R}}^{d_{\mathit{\xi }}}$, we derive ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k$ as (34) ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right).$(34) Following from (27(27) $\begin{array}{rl}& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }(\mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }.\end{array}$(27) ), (35) $\begin{array}{rl}& -{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)\\ & =-\underset{\vartheta \downarrow 0}{lim}\frac{{\mathrm{\nabla }}_{\mathit{B}}\log p(S,Y|\mathit{B}+\vartheta \mathit{V})-{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)}{\vartheta }\end{array}$(35) (36) $\begin{array}{rl}& =\frac{1}{\nu }\underset{\vartheta \downarrow 0}{lim}\frac{\sum _{\ell =1}^k{\mathit{P}}_{\ell }\left(\right(\mathit{B}+\vartheta \mathit{V}){\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }-\sum _{\ell =1}^k{\mathit{P}}_{\ell }(\mathit{B}{\mathit{s}}_{\ell }-{\overline{\mathit{y}}}_{\ell })\otimes {\mathit{s}}_{\ell }}{\vartheta }\end{array}$(36) (37) $\begin{array}{rl}& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }\mathit{V}{\mathit{s}}_{\ell }\otimes {\mathit{s}}_{\ell }.\end{array}$(37) The directional derivative of the gradient of the negative log prior is derived from (32(32) $-{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)=\rho \mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}-\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.$(32) ), (38) $\begin{array}{rl}-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\rho {\delta }_{\mathit{V}}\left(\mathit{W}\right(\mathit{B}-{\hat{\mathit{B}}}_{k-1}\left)\mathit{W}\right)\\ & -\text{β}{\delta }_{\mathit{V}}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}{\delta }_{\mathit{V}}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(38) The Frobenius regularizer term from (38(38) $\begin{array}{rl}-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\rho {\delta }_{\mathit{V}}\left(\mathit{W}\right(\mathit{B}-{\hat{\mathit{B}}}_{k-1}\left)\mathit{W}\right)\\ & -\text{β}{\delta }_{\mathit{V}}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}{\delta }_{\mathit{V}}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(38) ) is derived as (39) $\begin{array}{rl}{\delta }_{\mathit{V}}\left(\mathit{W}\right(\mathit{B}-{\hat{\mathit{B}}}_{k-1}\left)\mathit{W}\right)& =\underset{\vartheta \downarrow 0}{lim}\frac{\mathit{W}(\mathit{B}+\vartheta \mathit{V}-{\hat{\mathit{B}}}_{k-1})\mathit{W}-\mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}}{\vartheta }\\ & =\mathit{WVW}.\end{array}$(39)

Following the same rationale for the log-barrier terms, (40) $\begin{array}{rl}& {\delta }_{\mathit{V}}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\mathit{B}+\vartheta \mathit{V})-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}\left)\right({\mathit{I}}_{d_{\xi }}+\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}){)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}+\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}{)}^{-1}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}-\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}+\mathcal{O}\left({\vartheta }^2\right)\left)\right(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1},\end{array}$(40) and, (41) $\begin{array}{rl}{\delta }_{\mathit{V}}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}& =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\hat{L}{\mathit{I}}_d-(\mathit{B}+\vartheta \mathit{V})){)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\hat{L}{\mathit{I}}_d-\mathit{B}\left)\right({\mathit{I}}_{d_{\xi }}-\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}){)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}-\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}{)}^{-1}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}+\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}+\mathcal{O}\left({\vartheta }^2\right)\left)\right(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}.\end{array}$(41) Substituting (39(39) $\begin{array}{rl}{\delta }_{\mathit{V}}\left(\mathit{W}\right(\mathit{B}-{\hat{\mathit{B}}}_{k-1}\left)\mathit{W}\right)& =\underset{\vartheta \downarrow 0}{lim}\frac{\mathit{W}(\mathit{B}+\vartheta \mathit{V}-{\hat{\mathit{B}}}_{k-1})\mathit{W}-\mathit{W}(\mathit{B}-{\hat{\mathit{B}}}_{k-1})\mathit{W}}{\vartheta }\\ & =\mathit{WVW}.\end{array}$(39) ), (40(40) $\begin{array}{rl}& {\delta }_{\mathit{V}}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\mathit{B}+\vartheta \mathit{V})-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}\left)\right({\mathit{I}}_{d_{\xi }}+\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}){)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}+\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}{)}^{-1}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}-\vartheta (\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}+\mathcal{O}\left({\vartheta }^2\right)\left)\right(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}}{\vartheta }\\ & =-(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1},\end{array}$(40) ), and (41(41) $\begin{array}{rl}{\delta }_{\mathit{V}}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}& =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\hat{L}{\mathit{I}}_d-(\mathit{B}+\vartheta \mathit{V})){)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{\left(\right(\hat{L}{\mathit{I}}_d-\mathit{B}\left)\right({\mathit{I}}_{d_{\xi }}-\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}){)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}-\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}{)}^{-1}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =\underset{\vartheta \downarrow 0}{lim}\frac{({\mathit{I}}_{d_{\xi }}+\vartheta (\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}+\mathcal{O}\left({\vartheta }^2\right)\left)\right(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}-(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}}{\vartheta }\\ & =(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_d-\mathit{B}{)}^{-1}.\end{array}$(41) ) back into (38(38) $\begin{array}{rl}-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)& =\rho {\delta }_{\mathit{V}}\left(\mathit{W}\right(\mathit{B}-{\hat{\mathit{B}}}_{k-1}\left)\mathit{W}\right)\\ & -\text{β}{\delta }_{\mathit{V}}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}{\delta }_{\mathit{V}}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(38) ), we get the directional derivative of the gradient of the negative log prior, (42) $\begin{array}{rl}& -{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)=\rho \mathit{WVW}+\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\\ & +\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(42) Finally, substituting (37(37) $\begin{array}{rl}& =\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }\mathit{V}{\mathit{s}}_{\ell }\otimes {\mathit{s}}_{\ell }.\end{array}$(37) ) and (42(42) $\begin{array}{rl}& -{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right)=\rho \mathit{WVW}+\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\\ & +\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(42) ) into (34(34) ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(p\right(S,Y|\mathit{B}\left)\right)-{\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}\log \left(\pi \right(\mathit{B};\text{β}\left)\right).$(34) ), we obtain the directional derivative of the gradient of the negative log posterior, (43) $\begin{array}{rl}& {\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }\mathit{V}{\mathit{s}}_{\ell }\otimes {\mathit{s}}_{\ell }+\rho \mathit{WVW}\\ & +\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(43)

2.5. Hessian approximation inverse

To use the Newton's method to solve (9(9) ${\mathit{B}}_k^{\ast }:=\{\begin{array}{rl}\arg \underset{\mathit{B}}{min}& \left({\mathcal{L}}_k\right(\mathit{B}):=-\log (\pi (\mathit{B}|S,Y)\left)\right)\\ \mathrm{subjectto}& {\lambda }_{min}\left(\mathit{B}\right)\ge \tilde{\mu }\\ & {\lambda }_{max}\left(\mathit{B}\right)\le L\\ & \mathit{B}={\mathit{B}}^{\mathrm{\top }},\end{array}$(9) ), one needs to either solve the linear system ${\mathit{B}}_k{\mathit{d}}_k={\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_k\right)$, or find the inverse of the Hessian approximation ${\mathit{H}}_k={\mathit{B}}_k^{-1}$, thus ${\mathit{d}}_k={\mathit{H}}_k{\mathrm{\nabla }}_{\mathit{\xi }}F\left({\mathit{\xi }}_k\right)$.

To avoid inverting ${\hat{\mathit{B}}}_k$, we propose to use the Newton–Raphson method to find an approximation ${\hat{\mathit{H}}}_k$ such that $\Vert {\hat{\mathit{B}}}_k{\hat{\mathit{H}}}_k-{\mathit{I}}_{d_{\xi }}{\Vert }_F\le \mathrm{to}{l}_{\mathit{H}}$; we can start the Newton–Raphson from the previous Hessian inverse approximation, thus saving a significant computational time. The Newton–Raphson method for this problem is given by (44) ${\hat{\mathit{H}}}_k^{i+1}={\hat{\mathit{H}}}_k^i-\left(R^{\prime }\right({\hat{\mathit{H}}}_k^i\left){)}^{-1}R\right({\hat{\mathit{H}}}_k^i),$(44) where $R\left(\mathit{H}\right)={\mathit{B}}_k\mathit{H}-{\mathit{I}}_{d_{\xi }}$ and $R^{\prime }\left(\mathit{H}\right)={\mathit{B}}_k$, thus, (45) ${\hat{\mathit{H}}}_k^{i+1}=(2{\mathit{I}}_{d_{\xi }}-{\hat{\mathit{H}}}_k^i{\mathit{B}}_k){\hat{\mathit{H}}}_k^i.$(45) The recursive procedure is repeated until $\Vert R\left({\hat{\mathit{H}}}_k^i\right){\Vert }_F<\mathrm{to}{l}_{\mathit{H}}$. Then, we set ${\hat{\mathit{H}}}_k:={\hat{\mathit{H}}}_k^i$ as an approximation of ${\mathit{H}}_k$ in optimization. The convergence of the Newton–Raphson method depends on the initial guess of the inverse, ${\hat{\mathit{H}}}_k^0$; the method converges if the spectral norm of $R\left({\hat{\mathit{H}}}_k^0\right)$ is less than one [23]. We start from the inverse estimate in the previous iteration, ${\hat{\mathit{H}}}_k^0={\hat{\mathit{H}}}_{k-1}$. If the method does not converge, that is, if the residue increases, we reset ${\hat{\mathit{H}}}_k^i=2/(\hat{L}+\hat{\mu }){\mathit{I}}_{d_{\xi }}$. This starting point is similar to the one suggested by [23] for Hermitian positive definite ${\hat{\mathit{B}}}_k$, ${\hat{\mathit{H}}}_k^i=2/({\lambda }_{\max }+{\lambda }_{\min }){\mathit{I}}_{d_{\xi }}$, where ${\lambda }_{\max }$ and ${\lambda }_{\min }$ are, respectively, the largest and the smallest eigenvalues of ${\hat{\mathit{B}}}_k$.

3. Preconditioned stochastic gradient descent methods

Instead of the common update of SGD as presented in (2(2) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\mathit{\upsilon }}_k,$(2) ), we investigate the case where the stochastic gradient estimate is preconditioned with our Hessian inverse approximation, (46) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\hat{\mathit{H}}}_k{\mathit{\upsilon }}_k.$(46) If the function to be minimized is strongly convex, SGD with fixed step size converges linearly to a neighbourhood of the optimum [14].

Here, we expect the preconditioned SGD to converge to a neighbourhood of the optimum faster than vanilla SGD. One possible side effect of the preconditioning of the stochastic gradient estimate is that the noise in the gradient might be amplified [18]. Therefore, it is fundamental to control the error in the gradient estimates using variance control techniques, e.g. MICE [8], SVRG [16], SARAH [22]. Note that the variance control methods mentioned before are specializations of (2(2) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\mathit{\upsilon }}_k,$(2) ), only differing in the way that they compute ${\mathit{\upsilon }}_k$. Moreover, acknowledging the possibly prohibitive cost of updating the Hessian estimate every iteration, we update the Hessian after a given number of iterations or gradient evaluations.

Some variance control methods are based on computing differences between gradients at known points and are thus well-suited to be combined with our approach. The MICE algorithm [8] builds a hierarchy of iterates, then progressively computes gradient differences between subsequent iterates in this hierarchy to control the gradient estimator error. Since both the differences between the iterates and between the gradients are available, curvature pairs $\mathit{s}$ and $\overline{\mathit{y}}$ can be obtained without any extra cost in terms of memory or gradient evaluations. If we refer to the algorithm in (2(2) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\mathit{\upsilon }}_k,$(2) ) using MICE to obtain ${\mathit{\upsilon }}_k$ as SGD-MICE, preconditioning the gradient estimates by our inverse Hessian approximation as in (46(46) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\hat{\mathit{H}}}_k{\mathit{\upsilon }}_k.$(46) ) furnishes the SGD-MICE-Bay method. For the finite sum case, other control variates methods are available, e.g. SVRG [16], SARAH [22], SAGA [9], SAG [25]. However, SVRG and SARAH have an advantage over SAGA and SAG here in that curvature pairs can be obtained directly from the algorithms without extra memory or processing costs. We refer to SVRG and SARAH preconditioned by our inverse Hessian approximation as SVRG-Bay and SARAH-Bay, respectively.

3.1. Convergence of preconditioned stochastic gradient descent

Assume that the relative statistical error of the gradient estimator is uniformly bounded above by a constant ϵ and that an arbitrary lower-bound $\tilde{\mu }$ is imposed on the eigenvalues of the Hessian approximation, i.e. $\tilde{\mu }=\mathrm{ϵ\mu }$ with $ϵ\ge 1$. On the one hand, a large value of ε decreases the noise amplification, as previously discussed, and at the same time improves the stability of the preconditioned stochastic gradient descent. On the other hand, if μ is indeed the smallest eigenvalue of the true Hessian matrix of F, imposing the constraint that the Hessian approximation must have a smallest eigenvalue $\tilde{\mu }\gg \mu$ causes the Hessian approximation to differ from the true Hessian. In our analysis, we seek to consider the effect of the choice of parameters ε, ϵ, and the step size η.

Assumption 3.1

Regularity and convexity of the objective function

The objective function F has Lipschitz continuous second-order derivatives and is μ-convex. Namely, for all $\mathit{x},\mathit{y}\in \mathrm{\Xi }$, (47) $\begin{array}{rl}\Vert {\mathrm{\nabla }}^{2}F\left(\mathit{x}\right)-{\mathrm{\nabla }}^{2}F\left(\mathit{y}\right)\Vert & \le M\Vert \mathit{x}-\mathit{y}\Vert \end{array}$(47) (48) $\begin{array}{rl}F\left(\mathit{y}\right)& \ge F\left(\mathit{x}\right)+\mathrm{\nabla F}\left(\mathit{x}\right)\cdot (\mathit{y}-\mathit{x})+\frac{\mu }{2}{\Vert \mathit{y}-\mathit{x}\Vert }^2.\end{array}$(48)

Assumption 3.2

Unbiasedness and bounded statistical error of gradient estimator

The gradient estimator ${\mathit{\upsilon }}_k$ is an unbiased estimator of $\mathrm{\nabla F}\left({\mathit{\xi }}_k\right)$ with an upper bound ϵ on its relative statistical error, (49) $\begin{array}{rl}\mathbb{E}\left[{\mathit{\upsilon }}_k|{\mathit{\xi }}_k\right]& =\mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\end{array}$(49) (50) $\begin{array}{rl}\mathbb{E}\left[{\Vert {\mathit{\upsilon }}_k-\mathbb{E}\left[{\mathit{\upsilon }}_k|{\mathit{\xi }}_k\right]\Vert }^2|{\mathit{\xi }}_k\right]& \le {ϵ}^2{\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2.\end{array}$(50)

Assumption 3.3

Bounded error of Hessian approximation

The approximation $\widehat{{\mathit{B}}_k}$ of ${\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)$ has a lower-bound $\tilde{\mu }:=\mathrm{ϵ\mu }$ on its eigenvalues with $ϵ\ge 1$ and the product of the true Hessian by the approximation inverse $\widehat{{\mathit{H}}_k}$ results in the identity matrix added by an error that is bounded above in norm by a linear function of ε, (51) $\begin{array}{rl}\widehat{{\mathit{B}}_k}& ⪰\tilde{\mu }{\mathit{I}}_{d_{\mathit{\xi }}}:=\mathrm{ϵ\mu },ϵ\ge 1,\end{array}$(51) (52) $\begin{array}{rl}{\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}& =\mathit{I}+{\mathcal{E}}_k,\end{array}$(52) (53) $\begin{array}{rl}\Vert {\mathcal{E}}_k\Vert & \le C_{ϵ}(ϵ-1),\end{array}$(53) where $C_{ϵ}>0$ is a constant.

Theorem 3.1

Convergence of preconditioned stochastic gradient descent

If Assumptions 3.1, 3.2, and 3.3 are satisfied and (54) $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert <\frac{2{ϵ}^2{\mu }^2(1-ϵ-C_{ϵ}(1+ϵ)(ϵ-1)-\frac{r}{\eta })}{\mathrm{M\eta }(1+{ϵ}^2)}$(54) holds for all k>0, then the gradients of the iterates generated by the preconditioned stochastic gradient descent as described in (46(46) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\hat{\mathit{H}}}_k{\mathit{\upsilon }}_k.$(46) ) converge linearly wirh rate 1−r, (55) $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert \le (1-r{)}^{k+1}\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_0\right)\Vert .$(55)

Proof.

From Assumption 3.1, we obtain an upper bound for the norm of the gradient at ${\mathit{\xi }}_{k+1}$, (56) $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert \le \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)+{\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)({\mathit{\xi }}_{k+1}-{\mathit{\xi }}_k)\Vert +\frac{M}{2}{\Vert {\mathit{\xi }}_{k+1}-{\mathit{\xi }}_k\Vert }^2.$(56) Using the preconditioned SGD update in (46(46) ${\mathit{\xi }}_{k+1}={\mathit{\xi }}_k-{\eta }_k{\hat{\mathit{H}}}_k{\mathit{\upsilon }}_k.$(46) ), one can further develop the equation above as (57) $\begin{array}{rl}\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert & \le \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)-\eta {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\left(\widehat{{\mathit{H}}_k}{\mathit{\upsilon }}_k\right)\Vert +\frac{M{\eta }^2}{2}{\Vert \widehat{{\mathit{H}}_k}{\mathit{\upsilon }}_k\Vert }^2\\ & \le \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)-\eta {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}\mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\Vert \eta {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}({\mathit{\upsilon }}_k-\mathrm{\nabla F}({\mathit{\xi }}_k\left)\right)\Vert +\frac{M{\eta }^2}{2}{\Vert \widehat{{\mathit{H}}_k}{\mathit{\upsilon }}_k\Vert }^2\end{array}$(57) (58) $\begin{array}{rl}& \le \Vert {\mathit{I}}_{d_{\mathit{\xi }}}-\eta {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}\Vert \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\eta \Vert {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}\Vert \Vert {\mathit{\upsilon }}_k-\mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert +\frac{M{\eta }^2}{2}{\Vert \widehat{{\mathit{H}}_k}\Vert }^2{\Vert {\mathit{\upsilon }}_k\Vert }^2.\end{array}$(58) Taking the expectation conditioned on ${\mathit{\xi }}_k$ and using Assumption 3.2, (59) $\begin{array}{rl}& \mathbb{E}\left[\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert |{\mathit{\xi }}_k\right]\le \Vert {\mathit{I}}_{d_{\mathit{\xi }}}-\eta {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}\Vert \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\mathrm{\eta ϵ}\Vert {\mathrm{\nabla }}^{2}F\left({\mathit{\xi }}_k\right)\widehat{{\mathit{H}}_k}\Vert \Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert +\frac{M{\eta }^2}{2}(1+{ϵ}^2){\Vert \widehat{{\mathit{H}}_k}\Vert }^2{\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2.\end{array}$(59) Using Assumption 3.3, (60) $\begin{array}{rl}\mathbb{E}\left[\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert |{\mathit{\xi }}_k\right]& \le (\Vert {\mathit{I}}_{d_{\mathit{\xi }}}-\eta ({\mathit{I}}_{d_{\mathit{\xi }}}+{\mathcal{E}}_k)\Vert +\mathrm{\eta ϵ}\Vert {\mathit{I}}_{d_{\mathit{\xi }}}+{\mathcal{E}}_k\Vert )\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{\tilde{\mu }}^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2\\ & \le ((1-\eta )+\eta \Vert {\mathcal{E}}_k\Vert +\mathrm{\eta ϵ}(1+\Vert {\mathcal{E}}_k\Vert ))\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{\tilde{\mu }}^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2\\ & \le (1-\eta (1-ϵ)+\eta (1+ϵ)C_{ϵ}(ϵ-1))\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{ϵ}^2{\mu }^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2.\end{array}$(60) Letting the right hand side of the above inequality be strictly less than $(1-r)\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert$ results in (54(54) $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert <\frac{2{ϵ}^2{\mu }^2(1-ϵ-C_{ϵ}(1+ϵ)(ϵ-1)-\frac{r}{\eta })}{\mathrm{M\eta }(1+{ϵ}^2)}$(54) ). Taking expectation and unrolling the recursion complete the proof.

Theorem 3.1 shows that increasing the parameter ε that controls the smallest eigenvalue of the Hessian approximation improves the stability of preconditioned stochastic gradient descent. Moreover, controlling ε and the step size η to satisfy (54(54) $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert <\frac{2{ϵ}^2{\mu }^2(1-ϵ-C_{ϵ}(1+ϵ)(ϵ-1)-\frac{r}{\eta })}{\mathrm{M\eta }(1+{ϵ}^2)}$(54) ) for a given r guarantees global convergence. By letting $ϵ=0$, thus resulting in deterministic gradient observations; $ϵ=1$, meaning that $\tilde{\mu }=\mu$; and unitary step size $\eta =1$, one recovers the quadratic convergence of the Newton method with the condition that $\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert <2{\mu }^2/M$ (see (60(60) $\begin{array}{rl}\mathbb{E}\left[\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_{k+1}\right)\Vert |{\mathit{\xi }}_k\right]& \le (\Vert {\mathit{I}}_{d_{\mathit{\xi }}}-\eta ({\mathit{I}}_{d_{\mathit{\xi }}}+{\mathcal{E}}_k)\Vert +\mathrm{\eta ϵ}\Vert {\mathit{I}}_{d_{\mathit{\xi }}}+{\mathcal{E}}_k\Vert )\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{\tilde{\mu }}^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2\\ & \le ((1-\eta )+\eta \Vert {\mathcal{E}}_k\Vert +\mathrm{\eta ϵ}(1+\Vert {\mathcal{E}}_k\Vert ))\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{\tilde{\mu }}^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2\\ & \le (1-\eta (1-ϵ)+\eta (1+ϵ)C_{ϵ}(ϵ-1))\Vert \mathrm{\nabla F}\left({\mathit{\xi }}_k\right)\Vert \\ & +\frac{M{\eta }^2}{2{ϵ}^2{\mu }^2}(1+{ϵ}^2){\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2.\end{array}$(60) )).

Since our analysis depends on quantities that are hard to characterize in practice, we do not use it to control ε and η. Instead, we start optimization with $ϵ=(L+\mu )/\left(2\mu \right)$ and control the parameter ρ of the Frobenius regularization of the prior; the regularization term penalizes large changes in the Hessian approximation. Thus, at every Hessian approximation update, ε progressively decays. Figure 1 illustrates this behaviour for the Example in Section 4.1.

Figure 1. Value of $\tilde{\mu }$ versus iteration for SGD-MICE preconditioned with our Hessian inverse approximation applied to the minimization of the quadratic problem in example in Section 4.1.

3.2. Time complexity

Given that the interval between Hessian updates is set to be $\mathcal{O}\left(d_{\mathit{\xi }}\right)$, then the extra cost of finding the Hessian approximation using our Bayesian approach is $\mathcal{O}\left(d_{\mathit{\xi }}^2\right)$. The number of Newton iterations and conjugate gradient iterations do not depend on the problem's dimensionality but on the tolerance set for each algorithm. The cost of each conjugate gradient iteration is the cost of evaluating ${\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B},\text{β})$ in (43(43) $\begin{array}{rl}& {\delta }_{\mathit{V}}{\mathrm{\nabla }}_{\mathit{B}}{\mathcal{L}}_k(\mathit{B};\text{β})=\frac{1}{\nu }\sum _{\ell =1}^k{\mathit{P}}_{\ell }\mathit{V}{\mathit{s}}_{\ell }\otimes {\mathit{s}}_{\ell }+\rho \mathit{WVW}\\ & +\text{β}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}\mathit{V}(\mathit{B}-\hat{\mu }{\mathit{I}}_{d_{\xi }}{)}^{-1}+\text{β}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}\mathit{V}(\hat{L}{\mathit{I}}_{d_{\xi }}-\mathit{B}{)}^{-1}.\end{array}$(43) ), which is $\mathcal{O}\left(d_{\mathit{\xi }}^3\right)$ due to matrix multiplications. If, instead of using the conjugate gradient method to compute the Newton direction, one builds the fourth-order tensor of second-order derivatives of the log posterior with respect to the Hessian, the time complexity increases to $\mathcal{O}\left(d_{\mathit{\xi }}^4\right)$.

In each iteration, we have to compute $P_k$, the precision matrix for each curvature pair. If a full matrix is computed, we have a complexity of $\mathcal{O}\left(d_{\mathit{\xi }}^2\right)$ each iteration. If, however, we opt to ignore the off-diagonal terms, this complexity is improved to $\mathcal{O}\left(d_{\mathit{\xi }}\right)$. We can further decrease the memory overhead of our method by storing only a scalar, the trace of the covariance matrix, for each curvature pair. The reasoning behind this approach is that the trace of the covariance matrix works as a weight for the curvature pairs; in the curvature pairs where the gradient difference is better estimated, the respective curvature pair has a larger weight when computing the likelihood. When our Bayesian Hessian approximation is used with the MICE gradient estimator, the trace of the covariance matrix is available without any extra cost. In numerical tests, we observed that keeping the trace of the covariance matrix resulted in the best option regarding the balance between accuracy and computational time.

We do not discuss here the cost to compute ${\mathit{\upsilon }}_k$ satisfying condition (50(50) $\begin{array}{rl}\mathbb{E}\left[{\Vert {\mathit{\upsilon }}_k-\mathbb{E}\left[{\mathit{\upsilon }}_k|{\mathit{\xi }}_k\right]\Vert }^2|{\mathit{\xi }}_k\right]& \le {ϵ}^2{\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^2.\end{array}$(50) ), as it depends on the method used to estimate it. For example, if ${\mathit{\upsilon }}_k$ is a Monte Carlo estimator with adaptive batch sizes, the cost of evaluating ${\mathit{\upsilon }}_k$ is $\mathcal{O}\left({ϵ}^{-2}{\Vert \mathrm{\nabla }F\left({\mathit{\xi }}_k\right)\Vert }^{-2}\right)$. If the conditions of Theorem 1 are satisfied, the gradient sampling cost is $\mathcal{O}\left({ϵ}^{-2}\right(1-r{)}^{-2k})$.

4. Numerical examples

In this section, we solve some numerical examples with and without our Bayesian approach to validate its performance. In all cases, we present the convergence of the optimality gap, $F\left({\mathit{\xi }}_k\right)-F\left({\mathit{\xi }}^{\ast }\right)$, and in some cases, the convergence of the central-path Newton-CG used to find the Hessian approximations. The parameters of the Bayesian approach are kept fixed for all the numerical examples, namely, the logarithmic barrier parameter, $\text{β}={10}^{-2}$; the Frobenius norm regularizer parameter, $\rho ={10}^{-2}$; the tolerance on the norm of the gradient of the negative log posterior, $\mathrm{tol}={10}^{-6}$; the constraint on the lower bound of the eigenvalues of $\mathit{B}$, $\tilde{\mu }$, is set as the strong-convexity parameter μ; the relaxation parameter, $\alpha =1.05$; the number of central-path steps is six; and the factor of decrease of both $\mathrm{to}{l}_i$ and ${\text{β}}_i$ for each ith central-path step is $\gamma =2$. Moreover, to keep our Bayesian approach competitive in terms of computational time and memory allocation, we do not compute the full matrix $P_{\ell }$ for each curvature pair, using the inverse of the trace of the covariance matrix instead. Note that the trace of the covariance matrix of the difference between gradients is available from MICE since it is used to control the relative error of the gradient [8].

The results presented in this section are obtained using a python implementation of the proposed method whose source code can be found at GitHub.¹ Moreover, the source code of the MICE gradient estimator, also used in this section, can be found at Bitbucket.² A repository with all the code necessary to reproduce the results presented in this section is available at GitHub.³

4.1. Quadratic function

The first numerical test is a quadratic function with noise added to its gradient. We want to minimize ${\mathbb{E}}_{\mathit{\theta }}\left[f\right(\mathit{\xi },\mathit{\theta }\left)|\mathit{\xi }\right]$ with (61) $f(\mathit{\xi },\mathit{\theta })=\frac{1}{2}\mathit{\xi }\cdot \mathit{A}\mathit{\xi }-\mathit{b}\cdot \mathit{\xi }+\mathit{\xi }\cdot \sqrt{\mathbf{\Sigma }}\mathit{\theta },$(61) where $\mathit{A}$ is a square symmetric positive-definite matrix with $1⪯\mathit{A}⪯\kappa$, $\mathbf{\Sigma }$ is a symmetric and positive definite covariance matrix, $\mathit{\theta }$ is a vector of independent standard normal distributed random variables, and the vector $\mathit{b}$ is defined as $\mathit{b}={\mathbf{1}}_{d_{\mathit{\xi }}}$, where ${\mathbf{1}}_{d_{\mathit{\xi }}}$ is a vector of ones of dimension $d_{\mathit{\xi }}$. The eigenvalues of $\mathit{A}$ are sampled uniformly between 1 and κ, with necessarily one eigenvalue being 1 and another being κ. The gradient of f with respect to $\mathit{\xi }$ is (62) $\mathrm{\nabla f}(\mathit{\xi },\mathit{\theta })=\mathit{A}\mathit{\xi }-\mathit{b}+\sqrt{\mathbf{\Sigma }}\mathit{\theta },$(62) thus, the covariance matrix of $\mathrm{\nabla f}$ is $\mathbf{\Sigma }$.

We solve this problem for two different values of the condition number κ, namely ${10}^3$ and ${10}^6$, and for a fixed number of dimensions $d_{\mathit{\xi }}=10$. We compare the performances of different methods in the solution of this problem: SGD, SGD-MICE, SGD-Bay with decreasing step size, SGD-Bay with fixed step size, and SGD-MICE-Bay. In all cases, we run the optimization until ${10}^7$ gradients are evaluated. For SGD-Bay, since the number of gradient evaluations per iteration is fixed, we know the number of iterations a priori and thus choose the interval between Hessian updates to have 15 equally distributed updates. For SGD-MICE-Bay, we opt to have an update every $d_{\mathit{\xi }}$ iterations. The step sizes and batch sizes used for each method are presented in Table 1, where, for the cases when MICE is used, we show the tolerance ϵ on the relative statistical error of the gradient.

Table 1. Step sizes and batch sizes used for Example 4.1.

Method	Step size	Batch size
SGD	$\frac{1}{L\sqrt{k}}$	1
SGD-MICE	$\frac{2}{(L+\mu )(1+{ϵ}^2)}$	$ϵ=1$
SGD-Bay (decr.)	$\frac{1}{\sqrt{k}}$	10
SGD-Bay (fixed)	1	1000
SGD-MICE-Bay	$1/(1+{ϵ}^2)$	$ϵ=0.5$

Notes: Here, L is the Lipschitz constant of $\mathrm{\nabla f}$, μ is the strong-convexity constant, k is the iteration, and ϵ is the tolerance on the relative statistical error. 

In Figure 2, we present the optimality gap versus iteration for both $\kappa ={10}^3$ (left) and $\kappa ={10}^6$ (right) for the previously mentioned methods. In the cases where our Bayesian approach is used, the Hessian updates are illustrated as coloured squares with black edges. For the case where $\kappa =1000$, SGD-MICE-Bay achieved, in 39 iterations, a lower optimality gap than SGD in ${10}^7$ iterations. Moreover, SGD presents oscillations around the optimum, making it difficult to have a reliable optimum estimate. In the case of SGD-Bay, with and without the decreasing step size, we observe oscillations of larger amplitude due to the preconditioning with our Hessian approach. This behaviour is due to the lack of variance control of SGD; not even the decreasing step size is enough to avoid the oscillations. In the case of an even larger condition number, $\kappa ={10}^6$, SGD-MICE-Bay far outperforms the other methods while requiring 178 iterations. Both SGD and SGD-MICE, being first-order methods, cannot improve the optimality gap significantly once close to the optimum. As for SGD-Bay, the same oscillating behaviour is observed for both the fixed and decreasing step cases. From these results, it is clear that, in large condition numbers cases, using second-order information is necessary for an efficient convergence. Moreover, we observe that controlling the variance of the gradient estimator is of central importance to maintaining convergence in the stochastic setting, avoiding the amplification of noise due to the preconditioning.

Figure 2. Convergence of the optimality gap versus iteration with $\kappa =1000$ (left) and $\kappa ={10}^6$ (right) for Example 4.1. The Hessian updates are marked with coloured squares with black edges. In both cases, SGD-MICE-Bay reduced the optimality gap more than the other methods and in a much smaller number of iterations. Moreover, SGD-MICE-Bay did not oscillate around the optimum like both SGD-Bay cases.

One crucial aspect of our method is the cost to update the Hessian approximation, i.e. the cost of the Newton-CG method in this setting. For SGD-MICE-Bay in the $\kappa =1000$ case, we present in Figure 3 the Newton-CG convergence of the gradient norm $\Vert {\mathrm{\nabla }}_{\mathit{B}}\mathcal{L}\Vert$ per Newton iteration. The black numbers on top of each update denote the number of conjugate gradient steps needed to find the Newton direction. Each colour represents a different central-path step, i.e. at each colour change, both β and tol are decreased by a factor $\gamma =2$, with each tol denoted by a dash-dotted line. For this problem, the runtime of each of the Hessian updates is at fractions of seconds. Given that $d_{\mathit{\xi }}=10$, the Hessian matrix has 100 components. Still, in the worst case, 51 conjugate gradient iterations were needed to find the Newton direction, showing that the problem of finding $\hat{\mathit{B}}$ is well-conditioned. Moreover, the total number of Newton iterations in the three cases did not exceed 25.

Figure 3. For each of the three Hessian updates of SGD-MICE-Bay with $\kappa =1000$ in Example 4.1, we present the convergence of the gradient norm for the Newton-CG method. The number of curvature pairs available and the time taken are presented for each Hessian update. Each colour represents a different step of the central-path method with a different logarithmic barrier parameter β and a different residue tolerance. The tolerance for each central-path step is presented as a dash-dotted line. The number above each Newton iteration represents the number of conjugate gradient iterations needed to find the Newton direction.

Controlling the eigenvalues maxima and minima is essential to keep the stability and efficiency of preconditioned SGD. To motivate why our approach is better suited to stochastic optimization than using the BFGS formulae to approximate the Hessian, we, for every Hessian update of SGD-MICE-Bay, also compute a Hessian approximation using BFGS. In Figure 4, we present, for both $\kappa =1000$ (left) and $\kappa ={10}^6$ (right), the extreme eigenvalues for our Bayesian approach, for BFGS, and the true extremes. The same curvature pairs are used for both our Bayesian approach and BFGS. For both condition number cases, the BFGS method can keep the minimum eigenvalue above the actual value; however, the largest eigenvalue far exceeds the true largest eigenvalue. In contrast, our Bayesian approach keeps the eigenvalues of the Hessian approximation between the extreme values due to the logarithmic barrier constraints.

Figure 4. Maximum and minimum eigenvalues for both our Bayesian approach and BFGS for Example 4.1 with $\kappa =1000$ (left) and $\kappa ={10}^6$ (right). We use the same curvature pairs for both cases, obtained from SGD-MICE-Bay. The true extremes of the eigenvalues are presented as black lines. Dashed lines represent the largest eigenvalues, and dash-dotted lines represent the smallest eigenvalues. It is clear that our approach gets closer to the actual extreme eigenvalues than BFGS, which is farther from the smallest eigenvalue than our method and exceeds the value of the largest eigenvalue.

As motivated in (6(6) $\mathbb{E}\left[{\Vert {\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k-\mathbb{E}\left[{\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k|{\hat{\mathit{B}}}_k^{-1}\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right]\le {\Vert {\hat{\mathit{B}}}_k^{-1}\Vert }^2\mathbb{E}\left[{\Vert {\mathit{\upsilon }}_k-\mathbb{E}\left[{\mathit{\upsilon }}_k\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right],$(6) ) and Theorem 3.1, increasing $\hat{\mu }$ decreases the noise amplification of preconditioning the gradients by our Bayesian Hessian approximation. Figure 5 presents a comparison between SGD-Bay with a fixed step size for different values of $\hat{\mu }$. One thing to notice is the reduction of the amplitude of the oscillations as $\hat{\mu }$ increases. Moreover, we observe a clear linear convergence for the more conservative cases of $\hat{\mu }=1000$ and $\hat{\mu }=\mathrm{10,000}$. Comparing these results with the ones presented in Figure 2 for the same condition number $\kappa ={10}^6$, it is clear that increasing $\hat{\mu }$ improves the performance of SGD-Bay even without controlling the gradient statistical error or decreasing the step size. However, in the general convex case, SGD-Bay with a fixed step size will converge to a neighbourhood of the optimum and oscillate around it with a fixed amplitude; increasing $\hat{\mu }$ or decreasing the step size may reduce the amplitude of the oscillations but not circumvent them.

Figure 5. Optimality gap versus iteration for SGD-Bay with fixed step size on Example 4.1 for different values of $\hat{\mu }$. The coloured squares with black edges represent the Hessian updates of each method. The different values of $\hat{\mu }$ result in different amplitudes for the oscillations around the optimum of SGD-Bay; larger values of $\hat{\mu }$ result in smaller amplitudes of oscillations (cf. (6(6) $\mathbb{E}\left[{\Vert {\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k-\mathbb{E}\left[{\hat{\mathit{B}}}_k^{-1}{\mathit{\upsilon }}_k|{\hat{\mathit{B}}}_k^{-1}\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right]\le {\Vert {\hat{\mathit{B}}}_k^{-1}\Vert }^2\mathbb{E}\left[{\Vert {\mathit{\upsilon }}_k-\mathbb{E}\left[{\mathit{\upsilon }}_k\right]\Vert }^2|{\hat{\mathit{B}}}_k^{-1}\right],$(6) )).

4.2. Logistic regression

Training logistic regression models for classification of finite data is an empirical risk minimization problem, thus a problem of minimizing a sum of functions. When training a model on a large dataset, it is common to use SGD to reduce the iteration cost by subsampling the data points used in each iteration. Since inexact information about the gradient is used in each iteration, deterministic quasi-Newton methods usually fail to recover the Hessian and its inverse. Here, we use our Bayesian framework to recover the Hessian from noisy gradient observations and precondition the gradient estimates with the Hessian inverse. We compare plain SGD, SGD-MICE, SGD-MICE-Bay, SVRG, SVRG-Bay, SARAH, and SARAH-Bay on the task of training a binary classification logistic regression model with ${\ell }_2$ regularization on five different datasets, mushrooms, ijcnn1, w8a [17], cod-rna [29], and HIGGS [1].⁴

The ${\ell }_2$-regularized log loss function we seek to minimize is (63) $F\left(\mathit{\xi }\right)=\frac{1}{N}\sum _{i=1}^{N}f(\mathit{\xi },{\mathit{\theta }}_i=({\mathit{x}}_i,y_i\left)\right)=\frac{1}{N}\sum _{i=1}^N\log (1+\exp (-y_i\mathit{\xi }\cdot {\mathit{x}}_i\left)\right)+\frac{\lambda }{2}\Vert \mathit{\xi }{\Vert }^2,$(63) where each data point $({\mathit{x}}_i,y_i)$ is such that ${\mathit{x}}_i\in {\mathbb{R}}^{d_{\mathit{\xi }}}$ and $y_i\in \{-1,1\}$. For all datasets, we use a regularization parameter $\lambda ={10}^{-5}$.

For SGD, only a singleton sample is used, i.e. the minibatch size is set to 1. Also, as in the example in Section 4.1, a decreasing schedule is used for the step size. For SGD-MICE, the sample size is adaptively controlled to keep the relative error of the gradient estimates below $ϵ=0.5$ with a step size of $\eta =2/\left(\right(L+\mu \left)\right(1+{ϵ}^2\left)\right)$ [8]. For the HIGGS dataset, the choice of $ϵ=0.5$ was too conservative, thus we used $ϵ=0.8$. As for SGD-MICE-Bay, we used the same values of ϵ, however, since we precondition the gradient estimates with the inverse of the Hessian, we use a step size of $\eta =(1+{ϵ}^2{)}^{-1}$. Therefore, with the initial Hessian guess with diagonal entries $B_{\mathrm{ii}}=\frac{L+\mu }{2}$, SGD–MICE-Bay recovers exactly SGD-MICE with the optimal step size. Similarly, for SVRG and SARAH, we use a fixed step size of ${\eta }_k=0.1/L$, whereas for SVRG-Bay and SARAH-Bay we use ${\eta }_k=0.1$. For SVRG, SVRG-Bay, SARAH, and SARAH-Bay, we start with a full batch, iterate using a mini-batch size of five, and restart every two epochs.

In Table 2, we present the sizes, number of features, and condition numbers of the resulting optimization problems for the five datasets. The used datasets differ significantly in their characteristics.

Table 2. For each of the datasets used in Example 4.2, we present its data size, number of features, and condition number of the optimization problem.

dataset	size	num. of features	κ
mushrooms	8124	112	$5.25\times {10}^5$
ijcnn1	$\mathrm{49,990}$	22	$3.76\times {10}^4$
w8a	$\mathrm{49,749}$	$\mathbf{300}$	$2.91\times {10}^5$
cod-rna	$\mathrm{59,535}$	8	$\mathbf{6.48}\mathbf{\times }{\mathbf{10}}^{\mathbf{9}}$
HIGGS	$\mathbf{11,000,000}$	28	$2.50\times {10}^4$

In Figure 6, we present the convergence for the optimality gap versus iterations (left) and versus runtime in seconds (right) for the mushrooms dataset. For SGD, SGD-MICE, SVRG, and SARAH, we present their convergences as solid lines and their respective preconditioned cases as dashed lines. For all the methods, preconditioning the gradient estimates using our Bayesian approach improved the final optimality gap in at least one order of magnitude. In terms of convergence of the optimality gap versus runtime, we observe that computing the Hessian approximation adds an overhead to the optimization process; however, the improvement in convergence compensates for this overhead. In this example, SGD-MICE-Bay achieved an optimality gap ten times smaller than that found by SGD-MICE and more than two orders of magnitude smaller than the one found by SGD. Preconditioning the gradient estimates of SVRG and SARAH with our Bayesian Hessian inverse approximation reduced the final optimality gap found by more than one order of magnitude.

Figure 6. Convergence of the optimality gap versus iteration (left) and runtime (right) for Example 4.2, mushrooms dataset. The Hessian updates are marked with coloured squares with black edges.

In Figure 7, the convergence of the Newton-CG method is presented for each of the Hessian updates of SGD-MICE-Bay for the mushrooms dataset. Note that a relatively small number of Newton steps are needed for each central-path step; 8 in the worst case. Moreover, the number of conjugate gradient steps does not exceed 18, which is a remarkably small number of iterations considering that we have a problem with $112\times 112$ variables. This result indicates that the problem is well-conditioned.

Figure 7. Convergence of the Newton-CG method for each of the Hessian updates of SGD-MICE-Bay in Example 4.2, mushrooms dataset. The number of curvature pairs available and the time taken are presented for each Hessian update. Each colour represents a different step of the central-path method with a different logarithmic barrier parameter β and a different residue tolerance. The tolerance for each central-path step is presented as a dash-dotted line. The number above each Newton iteration represents the number of conjugate gradient iterations needed to find the Newton direction.

Figure 8 presents the convergence for the ijcnn1 dataset. Here, as in the case of the mushrooms dataset, our Bayesian approach improved the convergence in all cases. However, the most remarkable improvement is the one of SVRG-Bay, which achieved an optimality gap of more than two orders of magnitude smaller than the vanilla SVRG with just one Hessian update. Both SVRG-Bay and SGD-MICE-Bay achieved very similar values of optimality gap, with SVRG-Bay requiring less runtime than SGD-MICE-Bay. In Figure 9, we present the convergence of Newton-CG for the one Hessian update of SGD-MICE-Bay for the ijcnn1 dataset. Note that the time spent computing the Hessian approximation is 0.33 seconds; for comparison, computing the true Hessian took 0.54 seconds in an average of 100 evaluations.

Figure 8. Convergence of the optimality gap versus iteration (left) and runtime (right) for Example 4.2, ijcnn1 dataset. The Hessian updates are marked with coloured squares with black edges.

Figure 9. Convergence of the Newton-CG method for each of the Hessian updates of SGD-MICE-Bay in Example 4.2, ijcnn1 dataset. The number of curvature pairs available and the time taken are presented for each Hessian update. Each colour represents a different step of the central-path method with a different logarithmic barrier parameter β and a different residue tolerance. The tolerance for each central-path step is presented as a dash-dotted line. The number above each Newton iteration represents the number of conjugate gradient iterations needed to find the Newton direction.

The w8a dataset is the case with the largest number of features, 300, meaning that the Hessian to be approximated by our Bayesian approach has $300\times 300$ components. The convergence for the w8a dataset is presented in Figure 10. Coupling SVRG, SARAH, and SGD-MICE with our Bayesian approach improved the convergence in all cases, but it was more significant for SVRG and SARAH. However, in terms of runtime, SGD-MICE was able to perform better than its competitors, and SGD-MICE-Bay had a small edge over SGD-MICE. Figure 11 presents the convergence of Newton-CG for SGD-MICE-Bay on the w8a dataset. As in the case of the mushrooms dataset, a small number of Newton steps (a maximum of 24) and conjugate gradient steps (a maximum of 34) were needed despite the large dimensionality of the problem at $\mathrm{90,000}$ variables. Even with the large Hessian, the time spent to find an approximation (7.32 s in the first update and 15.9 s in the second) is competitive with the time of computing the true Hessian, of 8.52 s in an average of 100 evaluations. Remember that our approach has stability advantages over using the true Hessian inverse as a preconditioner due to the control of the noise amplification given by the logarithmic barrier on the smallest eigenvalue.

Figure 10. Convergence of the optimality gap versus iteration (left) and runtime (right) for Example 4.2, w8a dataset. The Hessian updates are marked with coloured squares with black edges.

Figure 11. Convergence of the Newton-CG method for each of the Hessian updates of SGD-MICE-Bay in Example 4.2, w8a dataset. The number of curvature pairs available and the time taken are presented for each Hessian update. Each colour represents a different step of the central-path method with a different logarithmic barrier parameter β and a different residue tolerance. The tolerance for each central-path step is presented as a dash-dotted line. The number above each Newton iteration represents the number of conjugate gradient iterations needed to find the Newton direction.

Training the logistic regression model to the cod-rna dataset is the problem with the largest condition number of those studied here, at $6.48\times {10}^9$. We expect our Bayesian approach to be able to assist the optimization methods in converging to the true solution despite the large condition number of the problem. The convergence of the optimality gap for the cod-rna dataset is presented in Figure 12. Note that SVRG, SARAH, and SGD-MICE get stuck after some iterations, improving very little for most of the optimization process. Their preconditioned counterparts, however, perform much better, being able to decrease the optimality gap continuously. When looking at the convergence of the optimality gap versus runtime for SGD-MICE-Bay, we observe a sharp decrease after the first Hessian update that follows until its stop. Figure 13 presents the convergence of Newton-CG for the four Hessian updates of SGD-MICE-Bay. In the second Hessian update, the third central-path step, the Newton-CG method got temporarily stuck for 12 iterations, leading us to believe that the parameters β and/or ρ were not properly set for this example. Still, the Newton-CG method could resume the Hessian approximation in a reasonable time, at 1.89 s, showing the robustness of our approach.

Figure 12. Convergence of the optimality gap versus iteration (left) and runtime (right) for Example 4.2, cod-rna dataset. The Hessian updates are marked with coloured squares with black edges.

Figure 13. Convergence of the Newton-CG method for each of the Hessian updates of SGD-MICE-Bay in Example 4.2, cod-rna dataset. The number of curvature pairs available and the time taken are presented for each Hessian update. Each colour represents a different step of the central-path method with a different logarithmic barrier parameter β and a different residue tolerance. The tolerance for each central-path step is presented as a dash-dotted line. The number above each Newton iteration represents the number of conjugate gradient iterations needed to find the Newton direction.

The HIGGS dataset contains $11\times {10}^6$ data points, being the largest studied here. Thus, since our Hessian approximation is independent of the data size, we expect to have a time advantage compared to computing the true Hessian. Figure 14 presents the convergence of the optimality gap for the HIGGS dataset. Here we observe that, in contrast to the previous results, our Bayesian approach worsened the performance of SARAH and SVRG. For SGD-MICE, however, our method improved the final optimality gap by more than one order of magnitude. Also, concerning runtime, SGD-MICE-Bay performed better than the other methods, despite doing 9 Hessian updates. The longest Hessian update was the last one, lasting 2.36 s. In comparison, the average time taken to compute a Hessian for this example is 130.58 s. One possible interpretation for the bad performance of SVRG-Bay and SARAH-Bay is that, contrary to SGD-MICE, SVRG and SARAH do not control the relative statistical error in the gradient estimates.

Figure 14. Convergence of the optimality gap versus iteration (left) and runtime (right) for Example 4.2, HIGGS dataset. The Hessian updates are marked with coloured squares with black edges.

5. Conclusion

We presented a Bayesian approach to approximate Hessians of functions given noisy observations of their gradients. This problem is of great importance in stochastic optimization, where pre-conditioning gradient estimates with the inverse of the Hessian matrix can improve the convergence of stochastic gradient descent methods. The proposed Bayesian approach takes into consideration not only the secant equations as in deterministic quasi-Newton methods but also the noise in the gradients. To mitigate the known effect of amplification of the gradient noise, we control the smallest eigenvalue of the Hessian approximation. To maximize the log-posterior of the Hessian matrix, we use a Newton-CG method with a central-path approach to impose the eigenvalue constraints. The numerical results presented show that our approach is not only interesting from the theoretical perspective but also results in practical advantage. In both a stochastic quadratic equation and in the training of a logistic regression model with ${\ell }_2$ regularization, using our Bayesian approximation of the Hessian matrix resulted in better convergence rates without a significant increase in runtime. For future research, we suggest using low-rank or diagonal approximations of the Hessians.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partially supported by the KAUST Office of Sponsored Research (OSR) under Award number OSR-2019-CRG8-4033 in the KAUST Competitive Research Grants Program Round 8, the Alexander von Humboldt Foundation, and the Flexible Interdisciplinary Research Collaboration Fund at the University of Nottingham Project ID 7466664.

Notes on contributors

André Gustavo Carlon

Dr André Gustavo Carlon completed his Bachelor's degree in Civil Engineering at the Universidade Federal de Santa Catarina (UFSC) in Brazil in 2013. He then pursued his Master's and PhD degrees at the same institution, defending his Master's thesis on structural optimization in 2015 and his Ph.D. dissertation on uncertainty quantification in 2019. Dr. Carlon has since been a postdoctoral fellow in the stochastic numerics group at King Abdullah University of Science and Technology (KAUST) in Saudi Arabia. His research focuses on stochastic optimization, structural optimization, reliability engineering, and optimal experimental design.

Luis Espath

Dr Luis Espath obtained his Engineering Diploma from the Pontifical Catholic University of Rio Grande do Sul (PUCRS) in Brazil in 2007. He then completed both his Master's and Doctoral degrees at the Federal University of Rio Grande do Sul (UFRGS) in Brazil, concluding his studies in 2013. Dr Espath served as a postdoctoral fellow at King Abdullah University of Science and Technology (KAUST) until 2019. He then joined RWTH Aachen University in Germany as a Research Scientist, where he defended his Habilitation thesis in Mathematics in the field of Theoretical Mechanics in 2022. His research interests include theoretical and computational mechanics, uncertainty quantification, stochastic optimization, and machine learning.

Raúl Tempone

Dr Raúl Tempone received his PhD from the Royal Institute of Technology in Stockholm, Sweden, in 2002. He subsequently held postdoctoral positions at the Institute for Computational Engineering and Sciences (ICES) at the University of Texas at Austin and served as an Assistant Professor at Florida State University in Tallahassee. Since 2009, Dr. Tempone has been a Full Professor at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. In 2018, he was awarded the Alexander von Humboldt Professorship at RWTH Aachen University in Germany. Dr. Tempone's research focuses on numerical analysis, specifically the development and analysis of efficient and robust numerical methods for problems involving stochastic models and differential equations in various fields, including computational mechanics, quantitative finance, biological and chemical modeling, and wireless communications.

Notes

1 https://github.com/agcarlon/bayhess

2 https://bitbucket.org/agcarlon/mice

3 https://github.com/agcarlon/bayhess_numerics

4 Obtained from LibSVM datasets at https://www.csie.ntu.edu.tw/cjlin/libsvmtools/datasets/binary.html, except for the HIGGS dataset that was obtained from https://archive.ics.uci.edu [11].