Stochastic L-S-BFGS Optimization Methods

• Stochastic L-S-BFGS is a limited-memory quasi-Newton method that uses noisy curvature and variance-reduced gradients to tackle stochastic, finite-sum, and composite optimization problems.
• It incorporates strategies like overlap sampling, Hessian-vector products, and damping techniques to ensure stability, positive definiteness, and controlled spectral properties.
• Empirical results demonstrate enhanced convergence in ill-conditioned and large-scale scenarios, balancing computational efficiency with robust curvature approximation.

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 (Moritz et al., 2015, Zhao et al., 2017). 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 (Mokhtari et al., 2014, Carlon et al., 10 Jul 2025).

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 (Mokhtari et al., 2014, Zhao et al., 2017, Song et al., 2024). 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 (Şimşekli et al., 2018, Moritz et al., 2015, Song et al., 2024).

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 (Li et al., 2018, Tong et al., 2020). 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

$v_t$0

with

$v_t$1

or, in loopless formulations,

$v_t$2

(Zhao et al., 2017, Song et al., 2024). In composite settings, the Euclidean proximal step is replaced by a metric proximal step involving $v_t$3 or $v_t$4: $v_t$5 (Song et al., 2024).

The literature uses closely related naming conventions. Some papers speak explicitly of “stochastic L-BFGS” (Moritz et al., 2015, Zhao et al., 2017), some of stochastic damped L-BFGS (Li et al., 2018, Lotfi et al., 2020), some of stochastic quasi-Newton or SQN-VR (Tong et al., 2020), and some of full-matrix stochastic BFGS methods that are presented as conceptual precursors to limited-memory versions (Mokhtari et al., 2014). 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 $v_t$6. Directly forming $v_t$7 from two unrelated stochastic gradients is widely treated as unstable. Early regularized stochastic BFGS work therefore enforced same-sample finite-difference curvature: $v_t$8 with the explicit requirement that both gradients use the same mini-batch $v_t$9, because the apparently natural alternative $H_t$0 does not guarantee the positive curvature condition needed for the update (Mokhtari et al., 2014). This same concern reappears in multi-batch L-BFGS, where the stable construction uses the overlap of consecutive batches,

$H_t$1

rather than inconsistent full-batch differences (Berahas et al., 2016).

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: $H_t$2 (Moritz et al., 2015). The improved convergence-rate framework uses the same idea with

$H_t$3

(Zhao et al., 2017). In single-loop proximal quasi-Newton methods, the correction pair is updated only every $H_t$4 iterations via

$H_t$5

(Song et al., 2024). This substitution is motivated by the claim that naive stochastic-gradient differencing yields biased or noisy curvature, whereas Hessian-vector products are more stable (Moritz et al., 2015, Song et al., 2024).

A third class uses smoothed second-order surrogates for gradient differences. In the Hessian/Fisher framework, the stochastic secant vector is replaced by

$H_t$6

where $H_t$7 is a mini-batch Hessian approximation or Fisher/generalized Gauss-Newton matrix (Liu et al., 2018). The paper distinguishes LBFGS-H, using $H_t$8, from LBFGS-F, using Fisher information; the latter is highlighted as especially convenient for least-squares and cross-entropy losses (Liu et al., 2018).

A more recent reinterpretation abandons exact secant fitting altogether. In the Bayesian formulation, each pair $H_t$9 is treated as noisy evidence about the inverse Hessian, and the update is a MAP estimator solving

$$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$0

so the method explicitly balances prior proximity to $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$1 against consistency with the noisy curvature observation (Carlon et al., 10 Jul 2025). 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 $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$2 curvature pairs and applies the standard two-loop recursion to compute $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$3 in $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$4 time and memory, rather than storing a $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$5 matrix (Li et al., 2018, Tong et al., 2020). In stochastic damped variants, the stored pair is $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$6 or $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$7 rather than $$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$8, reflecting the fact that the quasi-Newton information is modified before entering the memory (Lotfi et al., 2020, Li et al., 2018).

The limited-memory mechanism has been generalized in several directions. Stochastic Block BFGS replaces scalar secant conditions with a sketched block relation

$$x^\star := \arg\min_x \; \mathbb{E}_\theta[f(x,\theta)] = \arg\min_x F(x),$$9

where $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$0 is a sketch matrix with $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$1 (Gower et al., 2016). The limited-memory implementation stores block triples

$$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$2

and computes the inverse-Hessian action with a block two-loop recursion (Gower et al., 2016). 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 $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$3, with

$$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$4

and computes $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$5 by a manifold version of the two-loop recursion (Roychowdhury, 2017). The method retains constant step sizes and uses variance reduction rather than Wolfe line search (Roychowdhury, 2017).

The Bayesian L-S-BFGS formulation introduces a different memory object: it stores the most recent curvature triples

$$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$6

where $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$7 is a curvature precision parameter. The paper states that for a problem of dimension $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$8, the iteration cost of full S-BFGS is $$\min_{x\in \mathbb{R}^d} f(x)=\frac1n\sum_{i=1}^n f_i(x),$$9 and the cost of L-S-BFGS is $$\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),$$0; it also remarks that the literal nested-loop implementation is $$\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),$$1, which is effectively linear in $$\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),$$2 for fixed small memory $$\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),$$3 (Carlon et al., 10 Jul 2025).

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 (Berahas et al., 2019). Although that work is not primarily stochastic, it explicitly notes applicability whenever curvature pairs are generated from sampled or noisy gradients (Berahas et al., 2019). 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 (Mokhtari et al., 2014, Moritz et al., 2015, Song et al., 2024).

Paper	Setting	Main guarantee
(Mokhtari et al., 2014)	Convex stochastic objective with diminishing stepsizes	Almost sure convergence and $$\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),$$4 expected objective gap
(Moritz et al., 2015)	Strongly convex smooth finite sums	Linear convergence in expectation with constant stepsize
(Zhao et al., 2017)	Strongly convex smooth finite sums	Improved linear convergence, almost sure linear convergence, and total complexity $$\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),$$5
(Gower et al., 2016)	Block BFGS with SVRG	Global linear convergence under uniform subsampled-Hessian bounds
(Song et al., 2024)	Composite convex proximal quasi-Newton	Global linear convergence
(Lotfi et al., 2020)	Nonconvex variance-reduced damped L-BFGS	Almost sure convergence to a stationary point and a sublinear complexity bound
(Şimşekli et al., 2018)	Asynchronous nonconvex L-BFGS MCMC	Ergodic rate $$\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),$$6 and linear speedup under conditions
(Tong et al., 2020)	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

$$\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),$$7

and bounded metric eigenvalues

$$\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),$$8

it proves

$$\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),$$9

for an explicit $h$0 determined by $h$1 (Moritz et al., 2015). The later coordinate-transformation analysis sharpens both variance and spectral bounds, obtaining the rate

$h$2

together with an almost sure $h$3-linear convergence statement and the complexity bound

$h$4

(Zhao et al., 2017). That paper explicitly states that this bound is still worse than the classical SVRG complexity $h$5, identifying the extra $h$6 factor as an open issue (Zhao et al., 2017).

For purely stochastic convex objectives, regularized stochastic BFGS establishes almost sure convergence under diminishing stepsizes satisfying

$h$7

and proves

$h$8

provided $h$9 (Mokhtari et al., 2014). For composite nonsmooth convex optimization, the single-loop proximal quasi-Newton method proves global linear convergence using a Lyapunov function that couples $\ell_1$0, the function gap, and the loopless variance-reduction state (Song et al., 2024).

Nonconvex theory is weaker but now substantial. VARCHEN proves

$\ell_1$1

and, under an extra bounded-second-moment assumption on stochastic gradients,

$\ell_1$2

(Lotfi et al., 2020). The asynchronous SG-MCMC-based as-L-BFGS derives

$\ell_1$3

which yields the stated ergodic rate $\ell_1$4 for $\ell_1$5 (Şimşekli et al., 2018). In shared-memory asynchronous SQN-VR, the convergence factor is degraded by the delay parameter $\ell_1$6, but remains linear if the step size is chosen small enough (Tong et al., 2020).

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 $\ell_1$7 by $\ell_1$8 enforces the lower eigenvalue bound

$\ell_1$9

which implies

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$0

and, with the added identity bias,

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$1

(Mokhtari et al., 2014). The purpose is not merely formal definiteness; it is to prevent explosive stochastic behavior when curvature estimates become unreliable (Mokhtari et al., 2014).

Damping is the dominant safeguard in nonconvex limited-memory variants. Stochastic damped L-BFGS uses a damped curvature vector $$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$2 to enforce

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$3

thereby keeping every update well-defined and positive definite (Lotfi et al., 2020). The PyTorch implementation of SdLBFGS modifies the original method in two practical ways: it sets $$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$4 at each step, and it normalizes the computed search direction by its $$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$5 norm when no line search is used (Li et al., 2018). 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 (Li et al., 2018).

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

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$6

or

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$7

otherwise the update is skipped (Şimşekli et al., 2018, Berahas et al., 2016). The Hessian/Fisher framework uses the same idea in nonconvex settings,

$$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$8

to keep the approximate inverse Hessians bounded above and away from zero (Liu et al., 2018).

Conditioning control extends beyond positive definiteness. VARCHEN argues that SdLBFGS may preserve definiteness but still allow $$s_k = x_k - x_{k-1}, \qquad y_k = \nabla F(x_k)-\nabla F(x_{k-1}),$$9 to become ill-conditioned. Its remedy is to estimate lower and upper eigenvalue bounds for $H_k g_k$0, clip the initial scaling

$H_k g_k$1

and delete stored curvature pairs when the estimated spectrum leaves prescribed limits (Lotfi et al., 2020). 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 $H_k g_k$2, uses a likelihood scale $H_k g_k$3, and accepts pairs only under the stronger curvature threshold

$H_k g_k$4

(Carlon et al., 10 Jul 2025). 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 $H_k g_k$5 is induced by the intersection of worker sets across consecutive iterations (Berahas et al., 2016). 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 (Berahas et al., 2016).

AsySQN pushes the architecture further in shared memory. Each thread reads a stale version of $H_k g_k$6, computes a variance-reduced stochastic gradient,

$H_k g_k$7

applies the two-loop recursion, and writes back the update without global synchronization (Tong et al., 2020). 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 (Tong et al., 2020). The analysis models delay by $H_k g_k$8, $H_k g_k$9, and maximum delay $v_t$00, and shows that asynchrony worsens the rate constant but does not destroy linear convergence under bounded delay and suitable stepsizes (Tong et al., 2020).

A distinct asynchronous line embeds L-BFGS inside SG-MCMC. The as-L-BFGS method uses workers with local L-BFGS memories, stale variables $v_t$01, and updates

$v_t$02

combined with Gaussian noise and momentum (Şimşekli et al., 2018). The paper emphasizes that local memory is necessary because asynchronous workers cannot reliably share a single synchronized L-BFGS memory (Şimşekli et al., 2018). 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 (Şimşekli et al., 2018).

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 $v_t$03, and under the diagonal-output-Hessian condition the communication cost per distributed round is

$v_t$04

(Liu et al., 2018). 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 (Gower et al., 2016).

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 (Mokhtari et al., 2014). In the support-vector-machine experiment with squared hinge loss, it reaches the same objective value after about $v_t$05 samples where SGD needs about $v_t$06 for $v_t$07, and for $v_t$08 SGD becomes very slow while RES remains practical (Mokhtari et al., 2014).

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 (Moritz et al., 2015). 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 (Zhao et al., 2017). 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 (Gower et al., 2016).

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 $v_t$09 test accuracy on CIFAR10 versus about $v_t$10 for SGD and Adagrad, and achieves roughly $v_t$11 on MNIST, while PyTorch LBFGS performs poorly and the original SdLBFGS0 becomes unstable and crashes after a few epochs (Li et al., 2018). 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 (Lotfi et al., 2020).

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 $v_t$12 sooner than both rSVRG and VR-PCA on the last three datasets (Roychowdhury, 2017). 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,

$v_t$13

and is evaluated on six Atari 2600 games with batch sizes $v_t$14 and memory sizes $v_t$15 (Rafati et al., 2018).

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 (Mokhtari et al., 2014, Zhao et al., 2017, Lotfi et al., 2020). 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 (Moritz et al., 2015, Liu et al., 2018, Carlon et al., 10 Jul 2025).