---
title: Stochastic Proximal Point Methods
url: https://www.emergentmind.com/topics/stochastic-proximal-point-spp
type: topic
---

# Stochastic Proximal Point Methods

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 [1810.05633], [2405.14255], [2605.20805].

## 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 [1810.05633], [2301.03125].

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 [2308.09310].

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 [2405.14255].

The saddle-point specialization replaces minimization of a scalar loss by a proximal step on a convex-concave component. For
\[
f(x,y)=\frac1n\sum_{i=1}^n f_i(x,y),
\]
with \(f_i\) \(\mu\)-strongly convex in \(x\), \(\mu\)-strongly concave in \(y\), and \(L\)-smooth, the operator
\[
g_i(x,y)=\begin{bmatrix}\partial_x f_i(x,y)\\ -\partial_y f_i(x,y)\end{bmatrix}
\]
is used, and the Point-SAGA/SPP update is
\[
(x^{k+1},y^{k+1})=\prox_{f_j}^{\gamma}(p_j^k,q_j^k),
\]
with memory variables \(\{x_i^k,y_i^k\}\) and a sampled index \(j\) [1909.06946].

In nonlinear metric settings, Hilbert-space proximal mappings are replaced by metric-space resolvents. In a separable Hadamard space \((X,d)\), for a normal convex integrand \(f:E\times X\to(-\infty,+\infty]\),
\[
J_{\lambda,e}(x)=\arg\min_{y\in X}\Bigl\{f(e,y)+\tfrac1{2\lambda}d^2(x,y)\Bigr\},
\]
and the stochastic proximal-point iteration is
\[
x_{n+1}=J_{\lambda_n,\xi_{n+1}}(x_n).
\]
This is the natural generalization of stochastic proximality to geodesic metric spaces of nonpositive curvature [2605.20805].

## 2. Main analytical regimes and convergence guarantees

For constrained stochastic convex optimization with simple sampled constraint sets \(X_S\), a two-step SPP update first computes a proximal point
\[
y^k=z_{\mu_k}(x^k;S_k)=\arg\min_{z\in\mathbb R^n}\Bigl\{f(z;S_k)+\frac1{2\mu_k}\|z-x^k\|^2\Bigr\},
\]
and then projects onto \(X_{S_k}\). Under convex Lipschitz continuity and linear regularity of the constraints, averaged iterates satisfy an expected value-function gap of order \(\mathcal O(1/k^{1/2})\); under smooth strongly convex objectives, expected quadratic distance to the optimal solution can be of order \(\mathcal O(1/k)\); and a restarting variant overcomes step-size restrictions [1706.06297].

A related convex theory based on weak linear regularity analyzes
\[
x^{k+1}:=\arg\min_{z\in\mathbb R^n}F(z;\xi_k)+\frac1{2\mu_k}\|z-x^k\|^2.
\]
Under weak linear regularity, SPP satisfies a recurrence
\[
E[\dist_{X^*}^2(x^{k+1})]\le (1-\mu_k\sigma_{F,\mu_0})E[\dist_{X^*}^2(x^k)]+(S_F^*+2\beta)\mu_k^2,
\]
which yields \(O(1/k)\) distance convergence for \(\mu_k=O(1/k)\) and linear convergence under the interpolation assumption \(0\in\partial F(x^*;\xi)\) for all \(x^*\in X^*\) and \(\xi\) [1901.08663].

In convex composite risk minimization, the minibatch variant M-SPP studies
\[
F_t(w_t)\le \min_{w\in\mathcal W}\Bigl\{R_{S_t}(w)+(\gamma_t/2)\|w-w_{t-1}\|^2\Bigr\}+\varepsilon_t,
\]
with weighted output
\[
\bar w_T=\frac{2}{T(T+1)}\sum_{t=1}^T t\,w_t.
\]
Under \(L\)-smooth losses, a convex Lipschitz regularizer, and quadratic growth with constant \(\lambda\), Theorem 1(a) gives
\[
E[R(\bar w_T)-R^*]\le \frac{4\rho[R(w_0)-R^*]}{T^2}+\frac{2^9L}{\lambda\rho nT}R^*.
\]
The decomposition into a \(1/T^2\) bias term and a \(1/(nT)\) variance term is central to the later stability-based literature [2301.03125].

For weakly convex, nonsmooth, nonconvex objectives \(F(x)=E_{z\sim P}[f(x,z)]\), SPP is analyzed through proximal subproblems
\[
\arg\min_y\Bigl\{F(y)+\tfrac{1}{2\lambda}\|y-x\|^2\Bigr\},
\]
the Moreau envelope
\[
M_\lambda(x)=\min_y\Bigl\{F(y)+\tfrac1{2\lambda}\|y-x\|^2\Bigr\},
\]
and the stationarity measure
\[
\|\nabla M_\lambda(x)\|=\tfrac1\lambda\|x-\prox_{\lambda F}(x)\|.
\]
The proximally guided stochastic subgradient method of Davis and Grimmer is an inexact proximal-point iteration whose total oracle complexity is \(O(\varepsilon^{-2})\) for obtaining
\[
E\|x_R-\prox_{\gamma F}(x_R)\|^2\le \varepsilon,
\]
matching the \(O(\varepsilon^{-2})\) rate of smooth nonconvex SGD under the same stationarity measure [1707.03505].

In Hadamard spaces, weak rather than strong convergence is the generic conclusion. Under the Robbins–Monro condition
\[
\sum_{n=0}^\infty \lambda_n=+\infty,\qquad \sum_{n=0}^\infty \lambda_n^2<+\infty,
\]
a generalized-Lipschitz growth condition, and \(\arg\min F\neq\emptyset\), the stochastic proximal-point iterates are almost surely bounded, every weak cluster point lies in \(\arg\min F\), \(x_n\) converges weakly almost surely to an \(\arg\min F\)-valued random variable, and \(F(x_n)\to\min F\) almost surely [2605.20805]. Under strong monotonicity and an additional second-moment bound on Yosida approximates, a later metric-space result establishes
\[
E[d^2(x_n,x^*)]\le u/(n+2)
\]
for \(\lambda_n=1/(\bar\mu(n+2))\), together with explicit almost-sure tail bounds [2510.10697].

## 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),
\]
with \(e^k\) chosen so that \(v^k=\nabla f_{i_k}(x^k)-e^k\) is an unbiased estimator of \(\nabla F(x^k)\). Proximal SVRG uses
\[
e^k=\nabla f_{i_k}(\tilde x^s)-\nabla F(\tilde x^s),
\]
while proximal SAGA stores historical points \(y_i^k\) and takes
\[
e^k=\nabla f_{i_k}(y_{i_k}^k)-\tfrac1n\sum_{i=1}^n\nabla f_i(y_i^k).
\]
Under a generic variance-control recursion, the smooth convex case admits an \(O(1/k)\) rate for averaged iterates, and under the Polyak–Łojasiewicz condition the method admits a global linear rate [2308.09310].

In strongly convex settings without smoothness, an even broader SPP-LC template uses a correction vector \(h_k\) and a control state \(\phi_k\). With the Lyapunov function
\[
\Psi_k=\|x_k-x_\star\|^2+\alpha\,\sigma_k^2,
\]
Theorem 4.1 yields
\[
\Psi_k\le \theta^k\Psi_0+\frac{\zeta}{1-\theta}.
\]
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 \((1+\gamma^2\delta^2)/(1+\gamma\mu)^2\), and Point-SAGA attains exact linear convergence under the corresponding parameter specialization [2405.15941].

For finite-sum strongly convex-concave saddle-point problems, the Point-SAGA/SPP method combines a component-wise proximal step with memory. With
\[
\gamma=\frac{\sqrt{(n-1)^2\mu^2+4L^2n}-(n-1)\mu}{2L^2n},
\]
and Lyapunov function
\[
T^k=c\frac1n\sum_i\|g_i(x_i^k,y_i^k)-g_i(x^*,y^*)\|^2+\Bigl\|\begin{bmatrix}x^k-x^*\\ y^k-y^*\end{bmatrix}\Bigr\|^2,
\]
the method satisfies
\[
E[T^{k+1}]\le \alpha T^k,\qquad \alpha=\frac1{1+\mu\gamma}<1.
\]
Its iteration complexity is
\[
k=O\bigl((n+\kappa\sqrt n)\log(1/\varepsilon)\bigr),
\]
while vanilla SAGA/SVRG need \(O((n+\kappa^2)\log(1/\varepsilon))\) and Catalyst-accelerated SAGA/SVRG require \(O((n+\kappa\sqrt n)\log\kappa\log(1/\varepsilon))\) [1909.06946].

The Bregman generalization replaces Euclidean quadratic regularization by
\[
D_h(x,y)=h(x)-h(y)-\langle \nabla h(y),x-y\rangle,
\]
and updates
\[
x_{k+1}=\arg\min_{x\in H}\Bigl\{f_{i_k}(x)+\frac1{\eta_k}D_h(x,x_k)\Bigr\}.
\]
Variance-reduced Bregman variants include BSAPA and BLSVRP, with sublinear \(O(1/K)\) convergence in the convex case and linear \(O(q^k)\) convergence in the relatively strongly convex case. These results are also presented as recovering variance-reduced Bregman SGD in a unified way [2510.16655].

## 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
\[
\hat z_k=\arg\min_{x\in\mathbb R^d}\Bigl\{\phi(x)+\frac1{2\lambda}\|x-\bar z_{k-1}\|^2\Bigr\},
\]
then calls a subproblem solver PSS \(n\) times, and finally applies a probability booster PB. Proposition 3.4 shows that PSS yields
\[
E[\phi^\lambda(y_{I+1})-\phi^\lambda(\hat x)]\le \alpha^I(\sigma D+(LD^2)/2)+\sigma^2/I,
\]
and Theorem 6.3 gives overall stochastic-gradient complexity
\[
O\!\Bigl(\max\{\kappa\log(\kappa/\varepsilon),(\kappa\sigma^2)/(\mu\varepsilon)\}\cdot \log(1/p)\cdot \log(1/\varepsilon)\Bigr)
\]
for a high-probability guarantee under bounded variance alone [2402.08992].

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 \(\mathcal V(\xi)=0\) by generalized Jacobian steps
\[
(\mathcal W(\xi^j)+\eta_j I)d^j=-\mathcal V(\xi^j),
\]
together with an Armijo line search on a strongly convex merit function. In the weakly convex case it achieves
\[
\frac1{m(S+1)}\sum_{s=0}^S\sum_{k=0}^{m-1}\alpha_k\,E\|F_{\rm nat}(x^{s,k})\|^2=O\!\Bigl(\frac1{m(S+1)}\Bigr),
\]
and in the strongly convex case it has linear convergence in expectation [2204.00406].

A related constrained minimax variant studies the augmented Lagrangian
\[
L(x,y,\lambda)=\varphi(x)+g(x)+f(x,y)-h(y)-\psi(y)+\langle \lambda,Ax+By+c\rangle
\]
and performs stochastic implicit proximal-point updates for \(x\), \(y\), and \(\lambda\). 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 \(q\)-linear convergence of the primal iterates and global \(r\)-linear convergence of the multipliers in expectation [2605.23488].

Inexactness is also analyzed abstractly. Under \(\phi\)-smoothness,
\[
\|\nabla f_\xi(x)-\nabla f_\xi(y)\|\le \phi(\|x-y\|,\|\nabla f_\xi(y)\|)\|x-y\|,
\]
exact SPPM satisfies
\[
\mathbb E[f(\hat x_k)]-f(x_\star)\le \frac{\tfrac2\gamma+\phi(\|x_0-x_\star\|,\tfrac1\gamma\|x_0-x_\star\|)}{2k}\|x_0-x_\star\|^2,
\]
and if \(f\) is \(\mu\)-strongly convex then
\[
\mathbb E\|x_k-x_\star\|^2\le \Bigl(1-\frac{\mu}{\tfrac2\gamma+\phi(\|x_0-x_\star\|,\tfrac1\gamma\|x_0-x_\star\|)}\Bigr)^k\|x_0-x_\star\|^2.
\]
For SPPM-inexact, the same iteration complexity is retained up to a constant-factor loss [2502.03401].

Momentum can also be inserted into the implicit step. In SPPAM,
\[
\begin{cases}
x^+_{t+1}=x_t-\gamma\nabla f(x^+_{t+1})+\beta(x_t-x_{t-1}),\\
x_{t+1}=x^+_{t+1}-\gamma\varepsilon_{t+1},
\end{cases}
\]
and Theorem 1 yields a two-step recursion whose spectral radius determines contraction. SPPA is unconditionally stable in the sense that arbitrarily large \(\gamma\) only improves the one-step factor \(1/(1+2\gamma\mu)\), while SPPAM admits a two-dimensional stable region in the \((\gamma,\beta)\)-plane and can have a strictly smaller one-step factor than SPPA under the stated acceleration condition [2111.06171].

## 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
\[
f_i(x,y)=\rho/2\|x\|^2-y^\top A_i x-\tfrac12 y^\top(C_i+\lambda I)y+y^\top b_i,
\]
each proximal-point step of Point-SAGA can be implemented in \(O(d)\) via a rank-two Woodbury update. Experiments on Mountain Car features with \(d=400\), \(n=5k\) or \(20k\) samples, and regularizers \(\rho=\lambda\in 10^{-5\ldots -7}\) evaluate the primal optimality gap \(\ell(x^k)-\ell(x^*)\) 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 [1909.06946].

In statistical learning, M-SPP is instantiated for Lasso regression and logistic regression. Numerical evidence on simulated Lasso data confirms that larger \(T\) leads to faster early convergence through the \(1/T^2\) bias term, that convergence slows as the noise level \(\sigma\) grows, and that two-phase M-SPP significantly outperforms vanilla M-SPP when \(T\ll n\). 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 [2301.03125].

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 \(\gamma\) over several orders of magnitude with \(\beta=0.9\). They report that SGD and SGDM converge only in narrow bands of \(\gamma\), that SPPA converges for a very wide range of \(\gamma\), and that SPPAM converges \(3\text{–}10\times\) faster while retaining the wide stability region of SPPA [2111.06171].

The semismooth-Newton literature emphasizes sparse and nonsmooth models. SNSPP is evaluated on \(\ell_1\)-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 [2204.00406].

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 [2605.23488].

## 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}=x_k-\gamma\nabla f_{\xi_k}(x_k)\), whereas SPPM uses the implicit update \(x_{k+1}=x_k-\gamma\nabla f_{\xi_k}(x_{k+1})\). 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 [2301.03125], [2308.09310], [2502.03401].

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 \(J_{\gamma A}\) to the expectation operator \(A\), whereas stochastic versions apply one resolvent of a random \(A_\xi\) per iteration. This makes the method highly scalable, and variance reduction restores exact linear convergence under strong monotonicity and expected similarity [2405.14255].

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 \(\eta_k\) to converge to the exact minimizer, whereas variance-reduced SPP schemes allow constant step-sizes and achieve \(O(1/k)\) or linear convergence, depending on the structural assumptions [2510.16655], [2308.09310].

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 \(f_i\), 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 [1909.06946], [2204.00406].

The assumption sets vary substantially across subfields. Smooth finite-sum convex analysis invokes \(L\)-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 [1707.03505], [2405.15941], [2605.20805].

Extensions stated in the literature include non-uniform sampling, arbitrary sampling, non-smooth components with sublinear \(O(1/K)\) 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 [1909.06946], [2510.16655], [2301.03125].

Source: https://www.emergentmind.com/topics/stochastic-proximal-point-spp