---
title: Stochastic L-S-BFGS Optimization Methods
url: https://www.emergentmind.com/topics/stochastic-l-bfgs-l-s-bfgs
type: topic
---

# Stochastic L-S-BFGS Optimization Methods

Stochastic L-BFGS, often written L-S-BFGS, designates a family of limited-memory quasi-Newton methods for stochastic, finite-sum, and mini-batch optimization in which the search direction is preconditioned by an inverse-Hessian approximation learned from noisy curvature information rather than exact full-batch secant pairs. In representative formulations, the iterate update is written as \(x_{t+1}=x_t-\eta H_t v_t\), where \(v_t\) is a stochastic or variance-reduced gradient estimator and \(H_t\) is built from a short memory of curvature information [1508.02087][1704.00116]. Across the literature, the defining technical issue is that stochastic gradients make the classical secant equation brittle; as a result, modern stochastic L-BFGS methods differ chiefly in how they estimate curvature, preserve positive definiteness, control the spectrum of the metric, and reduce gradient noise [1401.7625][2507.07729].

## 1. Problem formulations and defining structure

The core problem classes are the stochastic program
\[
x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),
\]
the finite-sum empirical-risk problem
\[
\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),
\]
and the composite convex model
\[
\min_{x\in\mathbb{R}^d} F(x):= f(x)+h(x), \qquad f(x)=\frac{1}{n}\sum_{i=1}^n f_i(x),
\]
where \(h\) is convex and possibly nonsmooth [1401.7625][1704.00116][2409.16971]. These formulations cover regularized empirical risk minimization, logistic regression, ridge regression, support-vector machines, matrix completion, MAP estimation, and regularized logistic regression with \(\ell_1\) or elastic-net penalties [1806.02617][1508.02087][2409.16971].

Classical deterministic L-BFGS maintains a limited-memory inverse-Hessian approximation from curvature pairs
\[
s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),
\]
and computes \(H_k g_k\) via the two-loop recursion rather than forming a dense matrix explicitly [1805.02338][2011.00667]. Stochastic L-BFGS inherits the same low-memory architecture, but replaces exact gradients by stochastic or variance-reduced estimators. In SVRG-type variants, a typical update is
\[
x_{s,t+1} = x_{s,t} - \eta H_r v_{s,t},
\]
with
\[
v_{s,t} = \nabla f_{B_{s,t}}(x_{s,t}) - \nabla f_{B_{s,t}}(x^s) + g_s,
\]
or, in loopless formulations,
\[
v_k = \nabla f_{i_k}(x_k)-\nabla f_{i_k}(w_k)+\nabla f(w_k)
\]
[1704.00116][2409.16971]. In composite settings, the Euclidean proximal step is replaced by a metric proximal step involving \(B_k\) or \(H_k=B_k^{-1}\):
\[
x_{k+1} = \operatorname{prox}^{B_k}_{\eta_k,h}\!\left(x_k-\eta_k H_k v_k\right)
\]
[2409.16971].

The literature uses closely related naming conventions. Some papers speak explicitly of “stochastic L-BFGS” [1508.02087][1704.00116], some of stochastic damped L-BFGS [1805.02338][2012.05783], some of stochastic quasi-Newton or SQN-VR [2011.00667], and some of full-matrix stochastic BFGS methods that are presented as conceptual precursors to limited-memory versions [1401.7625]. The common invariant is limited-memory quasi-Newton preconditioning under stochastic curvature information.

## 2. Stochastic curvature acquisition

The central methodological divergence from deterministic L-BFGS lies in the construction of \(y_k\). Directly forming \(y_k\) from two unrelated stochastic gradients is widely treated as unstable. Early regularized stochastic BFGS work therefore enforced same-sample finite-difference curvature:
\[
r_t = g(x_{t+1},\Theta_t)-g(x_t,\Theta_t), \qquad \tilde r_t = r_t-\delta s_t,
\]
with the explicit requirement that both gradients use the same mini-batch \(\Theta_t\), because the apparently natural alternative \(g(x_{t+1},\Theta_{t+1})-g(x_t,\Theta_t)\) does not guarantee the positive curvature condition needed for the update [1401.7625]. This same concern reappears in multi-batch L-BFGS, where the stable construction uses the overlap of consecutive batches,
\[
O_k = S_k \cap S_{k+1}, \qquad y_{k+1}=g_{O_k}(w_{k+1})-g_{O_k}(w_k),
\]
rather than inconsistent full-batch differences [1605.06049].

A second large class replaces stochastic gradient differencing by stochastic Hessian-vector products. In the linearly convergent finite-sum algorithm, curvature pairs are built from averaged iterates and subsampled Hessians:
\[
s_r = u_r-u_{r-1}, \qquad y_r = \nabla^2 f_{\mathcal T_r}(u_r)\,s_r
\]
[1508.02087]. The improved convergence-rate framework uses the same idea with
\[
s_r = \bar x_r-\bar x_{r-1}, \qquad y_r=\nabla^2 f_{T_r}(\bar x_r)\,s_r
\]
[1704.00116]. In single-loop proximal quasi-Newton methods, the correction pair is updated only every \(r\) iterations via
\[
\bar{x}^t:=\frac{1}{r}\sum_{j=k-r+1}^{k}x_j, \qquad s^t=\bar{x}^t-\bar{x}^{t-1}, \qquad y^t=\nabla^2 f_{\mathcal{S}^t}(\bar{x}^t)s^t
\]
[2409.16971]. This substitution is motivated by the claim that naive stochastic-gradient differencing yields biased or noisy curvature, whereas Hessian-vector products are more stable [1508.02087][2409.16971].

A third class uses smoothed second-order surrogates for gradient differences. In the Hessian/Fisher framework, the stochastic secant vector is replaced by
\[
y_k = B_k^{S_k}(w_{k+1}-w_k),
\]
where \(B_k^{S_k}\) is a mini-batch Hessian approximation or Fisher/generalized Gauss-Newton matrix [1807.05328]. The paper distinguishes LBFGS-H, using \(\nabla^2 F^{S_k}(w_k)\), from LBFGS-F, using Fisher information; the latter is highlighted as especially convenient for least-squares and cross-entropy losses [1807.05328].

A more recent reinterpretation abandons exact secant fitting altogether. In the Bayesian formulation, each pair \((s_k,y_k)\) is treated as noisy evidence about the inverse Hessian, and the update is a MAP estimator solving
\[
H_{k+1} = \arg\min_{H=H^\top}\left(\frac{1}{\rho}\|s_k-Hy_k\|_{W_l}^2+\frac12\|H-H_k\|_{F,W_{pr}}^2\right),
\]
so the method explicitly balances prior proximity to \(H_k\) against consistency with the noisy curvature observation [2507.07729]. This marks a conceptual shift from deterministic secant satisfaction to noise-aware curvature assimilation.

## 3. Limited-memory recursion and generalized update families

In its canonical form, stochastic L-BFGS keeps only the most recent \(M\) curvature pairs and applies the standard two-loop recursion to compute \(H_k v\) in \(O(Md)\) time and memory, rather than storing a \(d\times d\) matrix [1805.02338][2011.00667]. In stochastic damped variants, the stored pair is \((s_k,\hat y_k)\) or \((s_k,\bar y_k)\) rather than \((s_k,y_k)\), reflecting the fact that the quasi-Newton information is modified before entering the memory [2012.05783][1805.02338].

The limited-memory mechanism has been generalized in several directions. Stochastic Block BFGS replaces scalar secant conditions with a sketched block relation
\[
H_t \nabla^2 f(x_t) D_t = D_t,
\]
where \(D_t\in \mathbb{R}^{d\times q}\) is a sketch matrix with \(q\ll d\) [1603.09649]. The limited-memory implementation stores block triples
\[
(D_i,Y_i,\Delta_i), \qquad i=t+1-M,\dots,t,
\]
and computes the inverse-Hessian action with a block two-loop recursion [1603.09649]. In this sense, block/sketched stochastic quasi-Newton methods generalize stochastic L-BFGS by matching curvature on a subspace rather than a single direction.

On Riemannian manifolds, the correction pairs must be transported between tangent spaces. The stochastic Riemannian L-BFGS method forms transported pairs \((z_r,y_r)\), with
\[
y_r = \tilde{\nabla}f(u_r)-\Gamma_\gamma \tilde{\nabla}f(u_{r-1}),
\]
and computes \(H_r\nu_i^{t+1}\) by a manifold version of the two-loop recursion [1704.01700]. The method retains constant step sizes and uses variance reduction rather than Wolfe line search [1704.01700].

The Bayesian L-S-BFGS formulation introduces a different memory object: it stores the most recent curvature triples
\[
\{(s_\ell,y_\ell,p_\ell)\}_{\ell=k-r+1}^k,
\]
where \(p_\ell\) is a curvature precision parameter. The paper states that for a problem of dimension \(d\), the iteration cost of full S-BFGS is \(O(d^2)\) and the cost of L-S-BFGS is \(O(d)\); it also remarks that the literal nested-loop implementation is \(O(dr^2)\), which is effectively linear in \(d\) for fixed small memory \(r\) [2507.07729].

Memory compression itself has become an object of study. Displacement aggregation shows that if one stored step lies in the span of later steps, an old pair can be removed and its information redistributed into the remaining gradient displacements without changing the BFGS matrix exactly [1903.03471]. Although that work is not primarily stochastic, it explicitly notes applicability whenever curvature pairs are generated from sampled or noisy gradients [1903.03471]. A plausible implication is that future stochastic L-BFGS schemes need not restrict limited memory to “keep the newest pairs” if dependence structure among steps can be exploited.

## 4. Convergence theory and complexity

The convergence theory of stochastic L-BFGS is heterogeneous because it tracks the problem class: strongly convex finite sums, stochastic strongly convex objectives, composite convex models, nonconvex smooth objectives, and asynchronous stale-gradient regimes all admit different guarantees. The most common assumptions are uniform lower and upper bounds on subsampled Hessians, bounded second moments for stochastic gradients, and uniform spectral bounds on the quasi-Newton metric [1401.7625][1508.02087][2409.16971].

| Paper | Setting | Main guarantee |
|---|---|---|
| [1401.7625] | Convex stochastic objective with diminishing stepsizes | Almost sure convergence and \(O(1/t)\) expected objective gap |
| [1508.02087] | Strongly convex smooth finite sums | Linear convergence in expectation with constant stepsize |
| [1704.00116] | Strongly convex smooth finite sums | Improved linear convergence, almost sure linear convergence, and total complexity \(O((n+\kappa\kappa_H)d\log(1/\epsilon))\) |
| [1603.09649] | Block BFGS with SVRG | Global linear convergence under uniform subsampled-Hessian bounds |
| [2409.16971] | Composite convex proximal quasi-Newton | Global linear convergence |
| [2012.05783] | Nonconvex variance-reduced damped L-BFGS | Almost sure convergence to a stationary point and a sublinear complexity bound |
| [1806.02617] | Asynchronous nonconvex L-BFGS MCMC | Ergodic rate \(O(1/\sqrt N)\) and linear speedup under conditions |
| [2011.00667] | Asynchronous parallel SQN-VR | Global linear convergence under strong convexity |

The 2015 linearly convergent algorithm is the canonical strong-convexity result. Under uniform Hessian bounds
\[
\lambda I \preceq \nabla^2 f_{\mathcal T}(w) \preceq \Lambda I
\]
and bounded metric eigenvalues
\[
\gamma I \preceq H_r \preceq \Gamma I,
\]
it proves
\[
\mathbb E[f(w_k)-f(w_*)] \le \alpha^k \mathbb E[f(w_0)-f(w_*)]
\]
for an explicit \(\alpha<1\) determined by \(\eta,m,\gamma,\Gamma,\lambda,\Lambda\) [1508.02087]. The later coordinate-transformation analysis sharpens both variance and spectral bounds, obtaining the rate
\[
\mathbb{E}[f(x^s)-f(x^*)] \le \rho^s\,(f(x^0)-f(x^*))
\]
together with an almost sure \(R\)-linear convergence statement and the complexity bound
\[
O\!\left((n+\kappa\kappa_H)d\log(1/\epsilon)\right)
\]
[1704.00116]. That paper explicitly states that this bound is still worse than the classical SVRG complexity \(O((n+\kappa)d\log(1/\epsilon))\), identifying the extra \(\kappa_H\) factor as an open issue [1704.00116].

For purely stochastic convex objectives, regularized stochastic BFGS establishes almost sure convergence under diminishing stepsizes satisfying
\[
\sum_{t=0}^\infty \epsilon_t=\infty, \qquad \sum_{t=0}^\infty \epsilon_t^2<\infty,
\]
and proves
\[
\mathbb{E}[F(x_t)]-F(x^\star)\le \frac{C_0}{T_0+t}
\]
provided \(2\epsilon_0T_0\Gamma>1\) [1401.7625]. For composite nonsmooth convex optimization, the single-loop proximal quasi-Newton method proves global linear convergence using a Lyapunov function that couples \(\|x_k-x^\ast\|_{B_k}^2\), the function gap, and the loopless variance-reduction state [2409.16971].

Nonconvex theory is weaker but now substantial. VARCHEN proves
\[
\liminf \|\nabla f(x_k)\| = 0 \quad \text{with probability } 1,
\]
and, under an extra bounded-second-moment assumption on stochastic gradients,
\[
\lim \|\nabla f(x_k)\| = 0 \quad \text{with probability } 1
\]
[2012.05783]. The asynchronous SG-MCMC-based as-L-BFGS derives
\[
\bigl|\mathbb E\hat U_N-U^\star\bigr| = O\!\left(\frac{1}{Nh}+\max(l_{\max}h,h)+\frac1\beta\right),
\]
which yields the stated ergodic rate \(O(1/\sqrt N)\) for \(h\sim 1/\sqrt N\) [1806.02617]. In shared-memory asynchronous SQN-VR, the convergence factor is degraded by the delay parameter \(\tau\), but remains linear if the step size is chosen small enough [2011.00667].

## 5. Stabilization, damping, and conditioning control

Stochastic L-BFGS is often described as if positive definiteness alone were the decisive property. The literature is more stringent. Regularized stochastic BFGS makes the point explicitly: replacing \(r_t\) by \(\tilde r_t=r_t-\delta s_t\) enforces the lower eigenvalue bound
\[
H_{t+1}\succeq \delta I,
\]
which implies
\[
H_t^{-1}\preceq \frac{1}{\delta}I
\]
and, with the added identity bias,
\[
\Gamma I \preceq H_t^{-1}+\Gamma I \preceq \left(\Gamma+\frac{1}{\delta}\right)I
\]
[1401.7625]. The purpose is not merely formal definiteness; it is to prevent explosive stochastic behavior when curvature estimates become unreliable [1401.7625].

Damping is the dominant safeguard in nonconvex limited-memory variants. Stochastic damped L-BFGS uses a damped curvature vector \(\hat y_k\) to enforce
\[
s_k^\top \hat y_k \ge \eta \, s_k^\top B_{k+1}^0 s_k > 0,
\]
thereby keeping every update well-defined and positive definite [2012.05783]. The PyTorch implementation of SdLBFGS modifies the original method in two practical ways: it sets \(H_{k,0}=I\) at each step, and it normalizes the computed search direction by its \(\ell_2\) norm when no line search is used [1805.02338]. The paper attributes improved stability to these choices and reports that the original SdLBFGS0 can crash or produce NaNs whereas the modified implementation is more robust [1805.02338].

Several papers impose a cautious skip rule. In asynchronous nonconvex L-BFGS and in multi-batch L-BFGS, the curvature pair is accepted only if
\[
y_i^\top s_i \ge \epsilon \|s_i\|^2
\]
or
\[
y_k^\top s_k \ge \epsilon \|s_k\|^2,
\]
otherwise the update is skipped [1806.02617][1605.06049]. The Hessian/Fisher framework uses the same idea in nonconvex settings,
\[
y_k^Ts_k < \epsilon \|s_k\|^2 \;\Rightarrow\; H_{k+1}=H_k,
\]
to keep the approximate inverse Hessians bounded above and away from zero [1807.05328].

Conditioning control extends beyond positive definiteness. VARCHEN argues that SdLBFGS may preserve definiteness but still allow \(H_k\) to become ill-conditioned. Its remedy is to estimate lower and upper eigenvalue bounds for \(H_k\), clip the initial scaling
\[
H_{k+1}^0 = \left( \max(\underline\gamma_{k+1}, \min(\gamma_{k+1}, \overline\gamma_{k+1})) \right) I,
\]
and delete stored curvature pairs when the estimated spectrum leaves prescribed limits [2012.05783]. This directly contradicts the common simplification that damping alone solves the stability problem.

The Bayesian formulation introduces a different stabilization vocabulary. Instead of forcing the secant equation exactly, it regularizes every denominator by \(\rho/p_k\), uses a likelihood scale \(\rho\), and accepts pairs only under the stronger curvature threshold
\[
y_k^\top s_k \ge m\|s_k\|^2
\]
[2507.07729]. This suggests a general principle visible across the field: reliable stochastic quasi-Newton behavior typically requires not only limited memory, but explicit mechanisms for rejecting, damping, or statistically discounting noisy curvature.

## 6. Parallel, distributed, and asynchronous realizations

Stochastic L-BFGS is closely tied to systems questions because mini-batch curvature extraction is often used precisely to make quasi-Newton methods viable at scale. Multi-batch L-BFGS was designed for large batches and distributed computation. A master broadcasts the current iterate, workers compute local gradients on data blocks, and the master aggregates returned gradients; if some workers are slow or fail, their contributions are omitted, and the overlap set \(O_k\) is induced by the intersection of worker sets across consecutive iterations [1605.06049]. The paper reports MPI-based strong-scaling and weak-scaling behavior and argues that the method introduces minimal extra communication beyond large-batch gradient evaluation [1605.06049].

AsySQN pushes the architecture further in shared memory. Each thread reads a stale version of \(x\), computes a variance-reduced stochastic gradient,
\[
v_t = \nabla f_{i_t}(x_{D(t)}) - \nabla f_{i_t}(w_k) + \nabla f(w_k),
\]
applies the two-loop recursion, and writes back the update without global synchronization [2011.00667]. The paper’s novelty claim is that it parallelizes the whole stochastic L-BFGS procedure asynchronously, not merely gradient computation or the two-loop recursion [2011.00667]. The analysis models delay by \(D(t)\), \(D'(t)\), and maximum delay \(\tau\), and shows that asynchrony worsens the rate constant but does not destroy linear convergence under bounded delay and suitable stepsizes [2011.00667].

A distinct asynchronous line embeds L-BFGS inside SG-MCMC. The as-L-BFGS method uses workers with local L-BFGS memories, stale variables \((\theta_{n-l_n},u_{n-l_n})\), and updates
\[
\Delta \theta_{n+1}=H_{n+1}(\theta_{n-l_n})\,u_{n-l_n},
\]
combined with Gaussian noise and momentum [1806.02617]. The paper emphasizes that local memory is necessary because asynchronous workers cannot reliably share a single synchronized L-BFGS memory [1806.02617]. It also states that increasing the number of workers from 1 to 10 yields near-linear speedup in wall-clock time in the real distributed speedup experiment, while too many workers can make the method unstable unless the step size is reduced [1806.02617].

Block/sketched variants also have systems implications. The vector-free recursion in LBFGS-H/F rewrites the two-loop recursion in terms of dot products and a small matrix \(M\), and under the diagonal-output-Hessian condition the communication cost per distributed round is
\[
\mathcal{O}\big(d\log(\tau)+m^2\big)
\]
[1807.05328]. Stochastic Block BFGS similarly targets large-scale logistic regression and highlights sketch selection mechanisms such as previous search directions and self-conditioning as practically effective [1603.09649].

## 7. Applications and empirical behavior

Empirical work on stochastic L-BFGS is unusually broad. In convex stochastic optimization, regularized stochastic BFGS consistently beats SGD on synthetic quadratics, with the advantage growing as the condition number grows; for ill-conditioned problems, it can reduce the number of processed samples by about an order of magnitude [1401.7625]. In the support-vector-machine experiment with squared hinge loss, it reaches the same objective value after about \(315\) samples where SGD needs about \(1.74\times 10^3\) for \(n=4\), and for \(n=40\) SGD becomes very slow while RES remains practical [1401.7625].

Large-scale linear models are a recurring benchmark. The 2015 linearly convergent algorithm is evaluated on ridge regression on Millionsong, SVM on RCV1, and matrix completion on Netflix, and is reported to work well over a wide range of step sizes, often differing by several orders of magnitude [1508.02087]. The improved convergence-rate paper evaluates logistic regression on `rcv1` and ridge regression on `E2006-tf`, concluding that geometric sampling or averaging performs as well as the last-iterate rule and better than uniform selection, that subsampled gradient stabilization improves early-stage convergence, and that all stochastic L-BFGS variants outperform SVRG on those datasets [1704.00116]. Stochastic Block BFGS is tested on large-scale logistic regression problems from LIBSVM, including `gisette-scale`, `covtype`, `HIGGS`, `SUSY`, `epsilon`, `rcv1`, and `url-combined`, with the `prev` sketch described as the most robust overall [1603.09649].

Deep-learning-oriented implementations reach more mixed conclusions, but they remain consistent about stabilization. The PyTorch SdLBFGS implementation reports that modified SdLBFGS is stable on CIFAR10 and MNIST, achieves about \(66\%\) test accuracy on CIFAR10 versus about \(62\%\) for SGD and Adagrad, and achieves roughly \(98\%\) on MNIST, while PyTorch LBFGS performs poorly and the original SdLBFGS0 becomes unstable and crashes after a few epochs [1805.02338]. VARCHEN, which adds eigenvalue control to SdLBFGS-VR, is reported to be more robust than SdLBFGS-VR and SVRG on modified DavidNet on CIFAR-10, while performing comparably on MNIST logistic regression and a nonconvex SVM on RCV1 [2012.05783].

The application range extends beyond Euclidean supervised learning. Stochastic Riemannian L-BFGS with variance reduction is evaluated on Karcher mean computation for SPD matrices and leading-eigenvalue computation on the sphere, where it is reported to converge faster than Riemannian SVRG and, in the eigenvalue experiment, to reach errors around \(10^{-15}\) sooner than both rSVRG and VR-PCA on the last three datasets [1704.01700]. In deep reinforcement learning, a multi-batch line-search L-BFGS method replaces SGD/RMSProp in DQN-style training, uses overlap batches for the curvature pair,
\[
y_k = \nabla (w_{k+1})^{(O_k)} - \nabla (w_k)^{(O_k)},
\]
and is evaluated on six Atari 2600 games with batch sizes \(b\in\{512,1024,2048,4096,8192\}\) and memory sizes \(m\in\{20,40,80\}\) [1811.02693].

These experiments do not support a single universal verdict. They do, however, support a consistent pattern: stochastic L-BFGS is typically strongest when ill-conditioning is significant, when variance reduction or overlap-based curvature estimation is available, and when the quasi-Newton metric is aggressively stabilized against stochastic curvature noise [1401.7625][1704.00116][2012.05783]. A common misconception is that stochastic L-BFGS is merely “SGD with a two-loop recursion.” The empirical record indicates otherwise: performance depends materially on how curvature pairs are constructed, when updates are skipped, how the initial metric is scaled, whether eigenvalues are controlled, and whether gradient noise is reduced by SVRG-type mechanisms, overlap sampling, or Bayesian weighting [1508.02087][1807.05328][2507.07729].

Source: https://www.emergentmind.com/topics/stochastic-l-bfgs-l-s-bfgs