Stochastic Proximal Point Methods

• Stochastic Proximal Point (SPP) is a family of methods that solve implicit proximal subproblems using sample-wise regularization instead of explicit gradient steps.
• The approach enhances numerical stability and robustness, making it resilient to step-size variations and compatible with variance reduction, minibatching, and momentum.
• SPP extends to convex composites, saddle-point problems, and non-Euclidean settings, yielding scalable performance and improved convergence rates in both convex and weakly convex scenarios.

Stochastic Proximal Point (SPP) denotes a family of implicit stochastic optimization methods in which each iteration solves a sample-wise proximal or resolvent subproblem instead of taking an explicit stochastic gradient step. In Euclidean form, a basic update is

$$x_{k+1}=\arg\min_x\Bigl\{f(x;\zeta_k)+\frac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$

while in monotone-operator form it is

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$

The framework appears in stochastic convex optimization, convex composite risk minimization, finite-sum saddle-point problems, monotone inclusions, weakly convex nonconvex optimization, and Hadamard-space optimization. Across these settings, it is studied for numerical stability, robustness against imperfect tuning, and compatibility with variance reduction, minibatching, momentum, and inexact inner solves (Asi et al., 2018, Sadiev et al., 2024, Pischke, 20 May 2026).

1. Canonical formulations and operator viewpoints

The canonical stochastic proximal-point update for convex minimization samples a fresh realization and solves a regularized one-sample problem. In the finite-sample or incremental convex setting, the exact update is

$$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$

and in composite risk minimization one writes

$$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$

These forms make the method implicit: the new iterate is defined through a strongly regularized subproblem rather than through an explicit gradient step (Asi et al., 2018, Yuan et al., 2023).

In convex composite finite-sum minimization, the proximity operator is written as

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$

and variance-reduced SPP schemes use proximal updates of the form

$x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$

where $e^k$ is a correction term (Traoré et al., 2023).

For monotone inclusions, the basic object is the resolvent. If $B:X\rightrightarrows X$ is maximally monotone, then

$J_{\gamma B}:=(I+\gamma B)^{-1},$

and a stochastic proximal-point iteration applies one resolvent call of a randomly sampled operator: $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$ This formulation covers minimization, saddle-point, and equilibria problems, and it remains meaningful even when the sampled operators are set-valued (Sadiev et al., 2024).

The saddle-point specialization replaces minimization of a scalar loss by a proximal step on a convex-concave component. For

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$0

with $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$1 $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$2-strongly convex in $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$3, $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$4-strongly concave in $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$5, and $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$6-smooth, the operator

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$7

is used, and the Point-SAGA/SPP update is

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$8

with memory variables $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$9 and a sampled index $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$0 (Luo et al., 2019).

In nonlinear metric settings, Hilbert-space proximal mappings are replaced by metric-space resolvents. In a separable Hadamard space $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$1, for a normal convex integrand $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$2,

$$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$3

and the stochastic proximal-point iteration is

$$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$4

This is the natural generalization of stochastic proximality to geodesic metric spaces of nonpositive curvature (Pischke, 20 May 2026).

2. Main analytical regimes and convergence guarantees

For constrained stochastic convex optimization with simple sampled constraint sets $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$5, a two-step SPP update first computes a proximal point

$$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$6

and then projects onto $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$7. Under convex Lipschitz continuity and linear regularity of the constraints, averaged iterates satisfy an expected value-function gap of order $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$8; under smooth strongly convex objectives, expected quadratic distance to the optimal solution can be of order $$x_{k+1}=\underset{x\in\mathbb R^n}{\arg\min}\;\Bigl\{\,f_{i_k}(x)+\tfrac{1}{2\lambda_k}\|x-x_k\|^2\Bigr\},$$9; and a restarting variant overcomes step-size restrictions (Patrascu et al., 2017).

A related convex theory based on weak linear regularity analyzes

$$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$0

Under weak linear regularity, SPP satisfies a recurrence

$$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$1

which yields $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$2 distance convergence for $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$3 and linear convergence under the interpolation assumption $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$4 for all $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$5 and $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$6 (Patrascu, 2019).

In convex composite risk minimization, the minibatch variant M-SPP studies

$$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$7

with weighted output

$$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$8

Under $$w_t=\arg\min_{w\in\mathcal W}\Bigl\{\ell(w;z_t)+r(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}.$$9-smooth losses, a convex Lipschitz regularizer, and quadratic growth with constant $\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$0, Theorem 1(a) gives

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$1

The decomposition into a $\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$2 bias term and a $\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$3 variance term is central to the later stability-based literature (Yuan et al., 2023).

For weakly convex, nonsmooth, nonconvex objectives $\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$4, SPP is analyzed through proximal subproblems

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$5

the Moreau envelope

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$6

and the stationarity measure

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$7

The proximally guided stochastic subgradient method of Davis and Grimmer is an inexact proximal-point iteration whose total oracle complexity is $\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$8 for obtaining

$\prox_{\alpha\,f}(z)=\arg\min_{x\in H}\Bigl\{\,f(x)+\tfrac1{2\alpha}\|x-z\|^2\Bigr\},$9

matching the $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$0 rate of smooth nonconvex SGD under the same stationarity measure (Davis et al., 2017).

In Hadamard spaces, weak rather than strong convergence is the generic conclusion. Under the Robbins–Monro condition

$x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$1

a generalized-Lipschitz growth condition, and $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$2, the stochastic proximal-point iterates are almost surely bounded, every weak cluster point lies in $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$3, $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$4 converges weakly almost surely to an $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$5-valued random variable, and $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$6 almost surely (Pischke, 20 May 2026). Under strong monotonicity and an additional second-moment bound on Yosida approximates, a later metric-space result establishes

$x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$7

for $x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$8, together with explicit almost-sure tail bounds (Pischke, 12 Oct 2025).

3. Variance reduction, sampling, and acceleration

Variance reduction enters SPP through correction vectors that preserve the implicit character of the update while controlling stochastic error. A unified proximal variance-reduction framework considers

$x^{k+1}=\prox_{\alpha\,f_{i_k}}(x^k+\alpha e^k),$9

with $e^k$0 chosen so that $e^k$1 is an unbiased estimator of $e^k$2. Proximal SVRG uses

$e^k$3

while proximal SAGA stores historical points $e^k$4 and takes

$e^k$5

Under a generic variance-control recursion, the smooth convex case admits an $e^k$6 rate for averaged iterates, and under the Polyak–Łojasiewicz condition the method admits a global linear rate (Traoré et al., 2023).

In strongly convex settings without smoothness, an even broader SPP-LC template uses a correction vector $e^k$7 and a control state $e^k$8. With the Lyapunov function

$e^k$9

Theorem 4.1 yields

$B:X\rightrightarrows X$0

This theorem recovers plain SPP, arbitrary-sampling SPP, “Star” SPP, SPPM-GC, loopless SVRP, and Point-SAGA as special cases. In particular, “Star” SPP has exact linear convergence, SPPM-GC has rate $B:X\rightrightarrows X$1, and Point-SAGA attains exact linear convergence under the corresponding parameter specialization (Richtárik et al., 2024).

For finite-sum strongly convex-concave saddle-point problems, the Point-SAGA/SPP method combines a component-wise proximal step with memory. With

$B:X\rightrightarrows X$2

and Lyapunov function

$B:X\rightrightarrows X$3

the method satisfies

$B:X\rightrightarrows X$4

Its iteration complexity is

$B:X\rightrightarrows X$5

while vanilla SAGA/SVRG need $B:X\rightrightarrows X$6 and Catalyst-accelerated SAGA/SVRG require $B:X\rightrightarrows X$7 (Luo et al., 2019).

The Bregman generalization replaces Euclidean quadratic regularization by

$B:X\rightrightarrows X$8

and updates

$B:X\rightrightarrows X$9

Variance-reduced Bregman variants include BSAPA and BLSVRP, with sublinear $J_{\gamma B}:=(I+\gamma B)^{-1},$0 convergence in the convex case and linear $J_{\gamma B}:=(I+\gamma B)^{-1},$1 convergence in the relatively strongly convex case. These results are also presented as recovering variance-reduced Bregman SGD in a unified way (Traoré et al., 18 Oct 2025).

4. Implicit realizations, minibatching, inexactness, and momentum

A recurring practical issue is that the proximal subproblem may not admit a closed form. One line of work analyzes this directly. In convex composite stochastic optimization, the SPP outer loop of SPPM defines the exact proximal subproblem

$J_{\gamma B}:=(I+\gamma B)^{-1},$2

then calls a subproblem solver PSS $J_{\gamma B}:=(I+\gamma B)^{-1},$3 times, and finally applies a probability booster PB. Proposition 3.4 shows that PSS yields

$J_{\gamma B}:=(I+\gamma B)^{-1},$4

and Theorem 6.3 gives overall stochastic-gradient complexity

$J_{\gamma B}:=(I+\gamma B)^{-1},$5

for a high-probability guarantee under bounded variance alone (Liang, 2024).

In weakly convex composite optimization, semismooth Newton methods are used to implement the implicit SPP step. The semismooth Newton stochastic proximal-point algorithm with variance reduction solves a nonsmooth system $J_{\gamma B}:=(I+\gamma B)^{-1},$6 by generalized Jacobian steps

$J_{\gamma B}:=(I+\gamma B)^{-1},$7

together with an Armijo line search on a strongly convex merit function. In the weakly convex case it achieves

$J_{\gamma B}:=(I+\gamma B)^{-1},$8

and in the strongly convex case it has linear convergence in expectation (Milzarek et al., 2022).

A related constrained minimax variant studies the augmented Lagrangian

$J_{\gamma B}:=(I+\gamma B)^{-1},$9

and performs stochastic implicit proximal-point updates for $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$0, $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$1, and $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$2. The resulting SNmMSPP method combines SVRG-type variance reduction with semismooth Newton and Armijo line search, and under strong convexity-concavity, semismoothness, and full-row-rank assumptions it has global $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$3-linear convergence of the primal iterates and global $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$4-linear convergence of the multipliers in expectation (Zhu et al., 22 May 2026).

Inexactness is also analyzed abstractly. Under $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$5-smoothness,

$x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$6

exact SPPM satisfies

$x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$7

and if $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$8 is $x^{k+1}=J_{\gamma_k A_{\xi^k}}(x^k).$9-strongly convex then

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$00

For SPPM-inexact, the same iteration complexity is retained up to a constant-factor loss (Tovmasyan et al., 5 Feb 2025).

Momentum can also be inserted into the implicit step. In SPPAM,

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$01

and Theorem 1 yields a two-step recursion whose spectral radius determines contraction. SPPA is unconditionally stable in the sense that arbitrarily large $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$02 only improves the one-step factor $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$03, while SPPAM admits a two-dimensional stable region in the $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$04-plane and can have a strictly smaller one-step factor than SPPA under the stated acceleration condition (Kim et al., 2021).

5. Applications and empirical behavior

Policy evaluation is one of the best-developed SPP applications in saddle-point form. After reducing the empirical mean-squared projected Bellman error to

$x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$05

each proximal-point step of Point-SAGA can be implemented in $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$06 via a rank-two Woodbury update. Experiments on Mountain Car features with $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$07, $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$08 or $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$09 samples, and regularizers $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$10 evaluate the primal optimality gap $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$11 versus epochs and wall-clock time, and report that SPP consistently converges in fewer epochs and less time than SVRG, SAGA, and SVRG+Catalyst, especially as the condition number grows (Luo et al., 2019).

In statistical learning, M-SPP is instantiated for Lasso regression and logistic regression. Numerical evidence on simulated Lasso data confirms that larger $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$12 leads to faster early convergence through the $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$13 bias term, that convergence slows as the noise level $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$14 grows, and that two-phase M-SPP significantly outperforms vanilla M-SPP when $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$15. On gisette and covtype, M-SPP and M-SPP-TP converge faster and more stably than minibatch-SGD with the same total passes over data, especially with large minibatches (Yuan et al., 2023).

The stability advantage of implicitness is also reported for generalized linear models. Kim et al. test SPPA, SPPAM, SGD, and SGDM on linear and Poisson regression, sweeping $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$16 over several orders of magnitude with $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$17. They report that SGD and SGDM converge only in narrow bands of $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$18, that SPPA converges for a very wide range of $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$19, and that SPPAM converges $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$20 faster while retaining the wide stability region of SPPA (Kim et al., 2021).

The semismooth-Newton literature emphasizes sparse and nonsmooth models. SNSPP is evaluated on $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$21-regularized sparse logistic regression on MNIST, Gisette, Sido0, Covtype, Higgs, and Madelon, and on sparse Student–t regression on synthetic data and on Sido0 features. The reported metrics include objective gap, natural residual, test loss, runtime per epoch, and total gradient evaluations; the experiments state that SNSPP tolerates much larger step-sizes without divergence and is competitive with or faster than Prox-SVRG, SAGA, and AdaGrad (Milzarek et al., 2022).

Constrained stochastic minimax experiments provide an additional large-scale application. SNmMSPP is tested on adversarial network flow and constrained linear regression; the reported findings are that it outperforms deterministic multiplier-gradient-descent and several heuristics on random graphs with Gaussian cost-noise, that it is robust to step-size choice, and that small inner-iteration counts together with high-accuracy Newton solves give the best CPU-versus-accuracy balance (Zhu et al., 22 May 2026).

6. Relation to neighboring methods, strengths, limitations, and extensions

The most common comparison is with SGD and its variance-reduced descendants. SGD uses the explicit update $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$22, whereas SPPM uses the implicit update $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$23. Multiple papers describe SPP methods as far less sensitive to step-size selection than SGD, numerically stable, and robust against imperfect tuning; in several settings they tolerate much larger constant step-sizes than gradient-based counterparts (Yuan et al., 2023, Traoré et al., 2023, Tovmasyan et al., 5 Feb 2025).

Relative to deterministic proximal point methods, SPP replaces the full objective or mean operator by a sampled component. In monotone-inclusion language, deterministic PPM applies $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$24 to the expectation operator $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$25, whereas stochastic versions apply one resolvent of a random $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$26 per iteration. This makes the method highly scalable, and variance reduction restores exact linear convergence under strong monotonicity and expected similarity (Sadiev et al., 2024).

Several contributions formalize a central strength of the method: variance reduction can remove the need for vanishing step-sizes. Vanilla SPPA or BSPPA typically requires diminishing or vanishing $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$27 to converge to the exact minimizer, whereas variance-reduced SPP schemes allow constant step-sizes and achieve $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$28 or linear convergence, depending on the structural assumptions (Traoré et al., 18 Oct 2025, Traoré et al., 2023).

The principal limitation repeatedly stated in the literature is access to the proximal operator or resolvent of each sampled component. In the saddle-point paper, this is explicit: SPP requires access to the proximal operator of each $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$29, and while many applications such as AUC maximization and policy evaluation admit closed-form or efficient rank-one solves, other settings may not. A closely related limitation is that practical performance may depend on the efficiency of the inner solver, especially when the implicit step is realized by semismooth Newton or by an inner stochastic method (Luo et al., 2019, Milzarek et al., 2022).

The assumption sets vary substantially across subfields. Smooth finite-sum convex analysis invokes $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$30-smoothness, strong convexity, quadratic growth, or PL; monotone-inclusion results invoke strong monotonicity and expected similarity; weakly convex nonconvex results rely on prox-regularization and Moreau-envelope stationarity; Hadamard-space analyses use geodesic convexity, generalized-Lipschitz growth, and weak convergence notions. This suggests that “Stochastic Proximal Point” is best understood as a method family unified by implicit sample-wise regularization rather than by a single theorem or a single rate (Davis et al., 2017, Richtárik et al., 2024, Pischke, 20 May 2026).

Extensions stated in the literature include non-uniform sampling, arbitrary sampling, non-smooth components with sublinear $x^{k+1}=(I+\gamma_k A_{\xi^k})^{-1}(x^k).$31 rates, heterogeneous strong convexity/concavity scales via variable rescaling, Bregman geometry, sampling without replacement, momentum, and potentially non-convex/non-concave local saddle-point finding. In that sense, SPP has developed from a stochastic analogue of the classical proximal point method into a broad design principle for implicit stochastic optimization (Luo et al., 2019, Traoré et al., 18 Oct 2025, Yuan et al., 2023).