A NOVEL CONVERGENCE ANALYSIS FOR THE STOCHASTIC PROXIMAL POINT ALGORITHM
Year: 2021

Abstract:
In this paper, we study the stochastic proximal point algorithm (SPPA) for general empirical risk minimization (ERM) problems as well as deep learning problems. We present an efficient implementation of SPPA with minor modification for different problem definitions and we observe that efficiently implemented SPPA has faster and more stable convergence than the celebrated stochastic gradient descent (SGD) algorithm, and its many variations, for both convex and non-convex problems. Due to the fact that per-iteration update of SPPA is defined abstractly and has long been considered expensive, its convergence proof has not been wellstudied until recently. In this paper, we close the theoretical gap by providing its convergence for convex problems. Our proof technique is different from some of the recent attempts. As a result, we present a surprising result that SPPA for convex problems may converge arbitrarily fast, depending on how the step sizes are chosen. As a second contribution, we also show that for some of the canonical ERM problems and deep learning problems, each iteration of SPPA can be efficiently calculated either in closed form or closed to closed form via bisection—the resulting complexity is exactly the same as that of SGD. Real data experiments showcase its effectiveness in terms of convergence compared to SGD and its variants.

heading: 1 INTRODUCTION
text: It has been widely accepted that when training large-scale machine learning models, the training algorithm should act in a sample-by-sample manner in order to reduce computational and memory overhead—the size of the data set may be too large that calculating the full gradient information is too costly. Moreover, most machine learning problems does not have to be solved with very high accuracy, since the ultimate goal is not to fit the training data but rather to generalize the algorithm well such that the performance is decent on unseen data.
Most existing stochastic algorithms are based upon the stochastic gradient descent (SGD) framework (Bottou et al., 2018). SGD is extremely easy to implement and provides asymptotic convergence, although the convergence rate is generally slow (and subject to careful choice of step sizes). Various approaches have been proposed to accelerate the plain vanilla SGD. Reducing the variance of the stochastic gradient and introducing adaptive learning schemes are two main lines of research. SVRG (Johnson & Zhang, 2013) and SAGA (Defazio et al., 2014) (and their follow-up works such as (Defazio, 2016)) focus on reducing the variance of the stochastic gradient descent, at the cost of increased time or memory complexities (to the order of the entire data set). On the other hand, AdaGrad (Duchi et al., 2011) and Adam (Diederik P. Kingma, 2014) introduce adaptive learning schemes and effectively keep the algorithm fully stochastic and light-weight. Besides these practical improvements, theoretical progress has been made by on quantifying the best possible convergence rate using first-order information (Lei et al., 2017; Allen-Zhu, 2017; 2018a;b).

heading: 1.1 STOCHASTIC PROXIMAL POINT ALGORITHM (SPPA)
text: In this paper, we explore a different type of stochastic algorithm called the stochastic proximal point algorithm (SPPA), also known as stochastic proximal iterations (Ryu & Boyd, 2014) or incremental proximal point method (Bertsekas, 2011a;b). Consider the following optimization problem with the
objective function in the form of a finite sum of component functions
minimize w∈R3
1
= =∑ 8=1 58 (w) = (w). (1)
SPPA takes the following simple form:
Algorithm 1 Stochastic proximal point algorithm (SPPA) 1: repeat 2: randomly draw 8C uniformly from {1, . . . , =} 3: wC+1 ← argminw _C 58C (w) + (1/2)‖w − wC ‖2 = Prox_C 58C (wC ) 4: until convergence
Line 3 of Algorithm 1 calculates the proximal operator of the function _C 58 evaluated at wC , denoted as Prox_C 58 (wC ). This is the stochastic version of the proximal point algorithm, which dates back to Rockafellar (1976).
SPPA has an abstract per-iteration update rule and it acquires more information from the problem than solely the first order derivatives, which makes it not as universally applicable as SGD. Yet, thanks to the more information inquired, it is possible to obtain faster and more robust convergence guarantees. While some ‘accelerated’ versions of SGD demand additional time or space overhead to go over the entire data set, SPPA does not have any overhead and it is suitable to be used in a completely online setting, it performs well even the data samples are not revisited again.
To the best of our knowledge, people started studying the convergence behavior of SPPA only recently (Bertsekas, 2011a; Ryu & Boyd, 2014; Bianchi, 2016; Pătraşcu, 2020; Toulis et al., 2021). Somewhat surprisingly, the convergence analysis of SPPA draws little resemblance to its deterministic counterpart, the proximal point algorithm. This is unlike the case for SGD, of which the convergence analysis follows almost line-by-line to that of the subgradient method. Moreover, existing analysis of SPPA shows no improvement in terms of convergence rate, which seems counter-intuitive due to the nature of the updates. Most authors also accept the premise that the proximal operator is sometimes difficult to evaluate, and thus proposed variations to the plain vanilla version to handle more complicated problem structures (Wang & Bertsekas, 2013; Duchi & Ruan, 2018; Asi & Duchi, 2019; Davis & Drusvyatskiy, 2019).

heading: 1.2 CONTRIBUTIONS
text: The main contribution of this paper is to provide a completely novel convergence analysis of SPPA for general convex problems. This contribution, together with the efficient implementation strategies discussed in the appendix, results in great practical results in almost all of the classical empirical risk minimization (ERM) problems in large-scale machine learning. While there exists some convergence results of SPPA for convex problems, the novel analysis provided in this work requires minimal assumptions (nothing but the convexity of the loss functions). The new convergence analysis also shows great resemblance to its deterministic counterpart, which is not the case from other works.
In the appendix we will discuss how to efficiently compute the abstract per-iteration update rule for SPPA with complexity that is comparable to SGD-type methods. An interesting observation is that in a lot of ERM examples the resulting update has a similar form like SGD, but with a smartly chosen step size derived from the proximal update. For a class of smooth risks such as the logistic loss, we appeal to reformulate the optimization problem such that we get close to closed form solution via bisection. We also briefly discuss how to modify the algorithm to have the ability to handle nonsmooth regularization terms such as the !1 norm, using the stochastic Douglas-Rachford splitting approach.
Finally, we apply SPPA with efficient implementations to a large variety of classification and regression problems. Numerical experiments are conducted not only for linear classification and regression (in the appendix) models, but also nonconvex deep neural network problems. Although the convergence analysis provided in this paper does not cover nonconvex cases, empirical results suggest that it is still worth treating SPPA as an effective alternative for deep learning.

heading: 2 CONVERGENCE ANALYSIS
text: In this section we provide convergence analysis of SPPA for general convex loss functions equation 1 to a global minimum in expectation. In recent years there have been some work tackling the same problem, e.g., Bertsekas (2011a); Pătraşcu (2020); Toulis et al. (2021). In this paper, however, we provide new convergence analysis that is much easier to understand while requiring nearly no assumptions other than convexity. For this reason we believe the theoretical contribution is significant enough as it broadens the applicability of SPPA.
There is a well-known resemblance between proximal methods and gradient descent, assuming the loss function is differentiable: while a full gradient descent step takes the form wC+1 = wC − _C∇ (wC ), the definition of a full proximal step guarantees that _C∇ (wC+1) = wC+1 − wC , meaning that wC+1 = wC − _C∇ (wC+1). Therefore, one might expect that a well-established convergence analysis of SGD for nonconvex problems, for example (Bottou et al., 2018, §4.3), can be seamlessly applied to SPPA. However, when applied in a stochastic fashion, the situation is a little more complicated.
Consider E[∇ 58 (wC ) | wC ], where the expectation is taken over the sampling procedure conditioned on wC , for SGD we typically require ∇ 58 (wC ) to be an unbiased estimator of ∇ (wC ), which is easy to satisfy if 8C is uniformly sampled from {1, . . . , =} given wC . For SPPA then one needs to consider E[∇ 58 (wC+1) | wC ], again over the sampling procedure conditioned on wC . This is in fact difficult to quantify because the update wC+1 depends on the sample that is drawn from the data set. The equation E[∇ 58 (wC+1) | wC ] = ∇ (wC+1) does not make sense because wC+1 on the right-hand-side is still random conditioned on wC . It is for this reason that existing analysis of SPPA differs drastically from its deterministic counterpart PPA (Bertsekas, 2011a; Ryu & Boyd, 2014; Bianchi, 2016).
What we can show, however, is that the distribution of 8C is still uniform conditioned on wC+1 instead of wC , as formalized in the following proposition. Lemma 1. At every iteration of SPPA (Algorithm 1), we have the conditional probability
?(8C |wC+1) = 1/=.
The proof of Lemma 1 is relegated to the supplementary material. What Lemma 1 suggests is that given the current iterate wC+1, without knowing its predecessors wC and beyond, every component function 58 is equally likely to be picked as the proximal update that leads to wC+1. In other words, conditioned on the current iterate without knowing the past, we do not gain additional information about which component function 58 is more likely to have been selected. As it turns out, Lemma 1 significantly simplifies the following convergence analysis while showing resemblance to that of deterministic PPA. What is more, the resulting analysis requires no assumptions other than convexity: typical assumptions such as Lipschitz smoothness or bounded variance are not required.
Remark. Lemma 1 holds for C = 0 if the initialization w0 is random, and the distribution from which it is drawn has the sample space the same as the domain of the objective function equation 1. In practice w0 is usually drawn from a distribution that is independent from the cost function, say N(0, O). In this case, for a given value of w1, there is no chance of eliminating possible w0’s because any corresponding w0 is in the sample space of the initialization distribution.

heading: 2.1 GENERIC CONVEX CASE WITHOUT STRONG CONVEXITY
text: Under Lemma 1, we have the following proposition, which serves as the stepping stone for our main convergence results. Proposition 1. Suppose each 51, . . . , 5= is convex, then at iteration C of SPPA (Algorithm 1) we have
2_C ( (wC+1) − (w★)) ≤ E[‖wC − w★‖2 | wC+1] − ‖wC+1 − w★‖2, (2) where w★ denotes an optimal solution of equation 1.
Proof. We start by the equation
‖wC − w★‖2 = ‖wC − wC+1 + wC+1 − w★‖2
= ‖wC+1 − w★‖2 + 2(wC − wC+1)>(wC+1 − w★) + ‖wC − wC+1‖2. (3)
According to the definition wC+1 = argminw _C 58C (w) + (1/2)‖w − wC ‖2, we know that wC − wC+1 is a subgradient of _C 58C at wC+1. Therefore _C 58C (wC+1) + (wC − wC+1)>(w − wC+1) ≤ _8 58C (w) at any w. Let w = w★ and substitute it in equation 3, we have
‖wC − w★‖2 ≥ ‖wC+1 − w★‖2 + 2_C ( 58C (wC+1) − 58C (w★)) + ‖wC − wC+1‖2
≥ ‖wC+1 − w★‖2 + 2_C ( 58C (wC+1) − 58C (w★)). Finally, taking condition expectations on both sides given wC+1 (but not wC nor 8C )
E[‖wC − w★‖2 | wC+1] ≥ ‖wC+1 − w★‖2 + 2_CE[( 58 (wC+1) − 58 (w★) | wC+1] . According to Lemma 1, the conditional distribution over 8C is uniform, thus
E[ 58C (wC+1) − 58C (w★) | wC+1] = (wC+1) − (w★), and we obtain equation 2.
With the help of Proposition 1, the rest of the results follows straight-forwardly. Theorem 1. If all 51, . . . , 5= are convex, and an initialization w0 is picked from a distribution with the sample space the same as the domain of , then the sequence {wC } generated by SPPA (Algorithm 1) satisfies
lim )→∞ inf C≤)
E[ (wC )] − (w★) ≤ 1 2 ∑∞ C=1 _C E‖w0 − w★‖2, (4)
where w★ denotes an optimal solution of equation 1. The term E‖w0 − w★‖2 is a constant for any approach of initializing w0.
Proof. Taking total expectation over equation 2, we have
2_C (E[ (wC+1)] − (w★)) ≤ E[‖wC − w★‖2] − E[‖wC+1 − w★‖2] . Sum over all C = 1, 2, . . . , ) , we have
2 )∑ C=1 _C (E[ (wC+1)] − (w★)) ≤ E‖w0 − w★‖2 − E‖w) − w★‖2 ≤ E‖w0 − w★‖2.
On the left-hand-side we have (by definition) inf E[ (wC )] ≤ E[ (wg)] for any g. Applying that, dividing both sides by 2 ∑ C _C , and letting ) →∞, we get equation 4.
Remark. The proofs here take conditional expectations backwards, i.e., conditioned on wC+1 and average over the previous iteration. This may be a little counter-intuitive. However, there is nothing wrong mathematically—expectation is a linear operator so E[ + ] = E[ ] + E[ ] under all circumstances; if an inequality 5 (G) ≤ 6(G) holds for all values of G in its sample space, then E[ 5 (G)] ≤ E[6(G)] for any distribution over G since it is just a nonnegative sum on both sides. Theorem 1 states a generic results on the convergence of SPPA. The left-hand-side of equation 4 is, by definition, nonnegative; the right-hand-side, however, goes to zero if the infinite sum ∑ C _C →∞. This implies that the infimum if E[ (wC )] goes to zero for appropriately chosen step sizes. Somewhat surprisingly, the only assumption we made was that the functions are convex. The celebrated SGD, on the other hand, requires at least two more assumptions: Lipschitz smoothness and that the stochastic gradients have bounded variance. What is more, the flexible choice of _C means that we can make the convergence rate arbitrarily fast, by allowing _C to be increasing (rather than decreasing in most gradient-based methods). Of course this may lead to large variance of the sequence, which we may want to avoid in practice.
Here we provide the convergence rate for two commonly used step size rules. Corollary 1. If all 51, . . . , 5= are convex, then the sequence {wC } generated by SPPA (Algorithm 1) with a constant step size _C = _ satisfies
inf C≤)
E[ (wC )] − (w★) ≤ 1
2_) ‖w0 − w★‖2, (5)
where w★ denotes an optimal solution of equation 1.
The proof is straight-forward by substituting _C with _ in equation 4. This shows that a constant step size (regardless of its value) gives a O(1/)) sublinear convergence rate. Corollary 2. If all 51, . . . , 5= are convex, then the sequence {wC } generated by SPPA (Algorithm 1) with an increasing step size rule _C = C_ satisfies
inf C≤)
E[ (wC )] − (w★) ≤ 1
_) () + 1) ‖w0 − w★‖ 2, (6)
where w★ denotes an optimal solution of equation 1.
The proof comes from the arithmetic series 1+2+· · ·+) = ) () +1)/2. This shows that the increasing step size rule _C = C_ gives a O(1/)2) sublinear convergence rate, the same as Nesterov’s optimal gradient algorithm (2013).
While the expected convergence looks excellent, we notice that in practice the performance also depends on the variance of the sequence generated by SPPA, which could be large if the expected convergence rate is too fast. We can even have linear convergence rate by letting _C increase exponentially, but the variance would be so large that it makes little practical sense. As we will see in the experiment section, a constant step size rule still gives the best performance in most cases.

heading: 2.2 STRONGLY CONVEX CASE
text: Similar to most other algorithms, convergence can be significantly improved if the loss functions are strongly convex. Proposition 2. Suppose each 51, . . . , 5= is strongly convex with parameter `, and an initialization w0 is picked from a distribution with the sample space the same as the domain of , then at iteration C of SPPA (Algorithm 1) we have
2_C ( (wC+1) − (w★)) ≤ E[‖wC − w★‖2 | wC+1] − (1 + _C/`)‖wC+1 − w★‖2, (7) where w★ denotes an optimal solution of equation 1.
Proof. Following equation equation 3, we again have that wC − wC+1 is a subgradient of 58C , which satisfies the following inequality if it is strongly convex
_C 58C (wC+1) + (wC − wC+1)>(w − wC+1) + _C
2` ‖w − wC+1‖2 ≤ _8 58C (w)
at any w. Let w = w★ and substitute it in equation 3, we have
‖wC − w★‖2 ≥ ‖wC+1 − w★‖2 + 2_C ( 58C (wC+1) − 58C (w★)) + _C
` ‖w − wC+1‖2 + ‖wC − wC+1‖2
≥ (1 + _C/`)‖wC+1 − w★‖2 + 2_C ( 58C (wC+1) − 58C (w★)). Taking condition expectations on both sides given wC+1 (but not wC nor 8C ) and invoking Lemma 1, we obtain equation 7.
The main convergence of SPPA for strongly convex functions is the following. Theorem 2. If all 51, . . . , 5= are strongly convex with parameter `, then the sequence {wC } generated by SPPA (Algorithm 1) with constant step sizes _C = _ satisfies
inf C≤)
E[ (wC )] − (w★) ≤ 1/(2`)
(1 + _/`)) − 1E‖w0 − w★‖ 2, (8)
where w★ denotes an optimal solution of equation 1. The term E‖w0 − w★‖2 is a constant for any approach of initializing w0. For large enough ) , the denominator on the right-hand-side of equation 7 is almost (1 + _/`)) , which indicates a linear convergence rate.
Proof. Letting _C = _, taking total expectation over equation 7 and multiplying both sides by (1 + _/`)C , we have 2_(1 + _/`)C (E[ (wC+1)] − (w★)) ≤ (1 + _/`)CE[‖wC −w★‖2] − (1 + _/`)C+1E[‖wC+1 −w★‖2] .
Sum over all C = 1, 2, . . . , ) , we have 2_ )∑ C=1 (1+_/`)C (E[ (wC+1)] − (w★)) ≤ E‖w0−w★‖2−(1+_/`)) E‖w) −w★‖2 ≤ E‖w0−w★‖2.
On the left-hand-side we have (by definition) inf E[ (wC )] ≤ E[ (wg)] for any g. Replacing every E[ (wC )] with the infimum and applying the geometric series
)∑ C=0 (1 + _/`)C = (1 + _/`) ) − 1 _/` ,
we get equation 8.
Final remark. In this section, we provided convergence of SPPA in expectation, following the strategies utilized in the tutorial (Bottou et al., 2018) which is easy to understand. In the supplementary we also provide a stronger almost sure convergence using Martingales theorem.

heading: 2.3 EXTENSION TO STOCHASTIC DOUGLAS-RACHFORD SPLITTING
text: In many cases the ERM formulation includes a regularization term 6(w) as
minimize w∈R3
1
= =∑ 8=1 58 (w) + 6(w). (9)
The most widely adopted regularization include the Tikhonov regularization d‖w‖2, the lasso regularization d‖w‖1 to promote sparsity, and the group-lasso regularization to promote groupsparsity. The Tikhonov regularization can easily be incorporated into the ERM parts without affecting the efficiency of the proximal updates. However, for the nonsmooth ones it may not be ideal to spread it into the stochastic terms as 58 (w) + 6(w). We propose to modify the Douglas-Rachford splitting in order to handle situation like this, which keeps track of two sequences {wC } and {w̃C }:
Algorithm 2 Stochastic Douglas-Rachford splitting (SDRS) 1: repeat 2: wC+1 ← Prox_C6 (w̃C ) 3: randomly draw 8C uniformly from {1, . . . , =} 4: w̃C+1 ← w̃C + Prox_C 58C (2wC+1 − w̃C ) − wC+1 5: until convergence
There has been some attempts to study the convergence of SDRS in the name of stochastic ADMM (Ouyang et al., 2013; Zhong & Kwok, 2014; Huang et al., 2019) or online ADMM (Wang & Banerjee, 2012). However, most of them do not consider a full proximal update but rather a ‘linearized’ update similar in the form of a (stochastic) gradient step; Wang & Banerjee (2012) considers the full proximal update but in an online setting using regret bounds to prove convergence, which resulted in a slower convergence rate as is typical comparing incremental verses stochastic methods.

heading: 3 EXPERIMENTS
text: We show some real data experiments to demonstrate effectiveness of the proposed SPPA implementations in linear regression and binary classification problems. We compare our proposed algorithms with SGD and its three variants: SGD with momentum (Sutskever et al., 2013), AdaBelief (Zhuang et al., 2020), and Adam (Diederik P. Kingma, 2014), with the default settings suggested in their original papers. We used the PyTorch (Paszke et al., 2019) implementation of the algorithms and used the same platform for implementing our customized optimization algorithm. We ran SPPA with our proposed implementations with various _ values among {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1}. The performance was consistent with different settings. We ran all the experiments with 100 different seeds, 0-100, to show the robustness of the algorithms to different stochastic settings. We ran all the experiments in Google Colab Pro environment. 1
1https://colab.research.google.com/notebooks/intro.ipynb

heading: 3.1 LINEAR CLASSIFICATION
text: We perform binary classification using SVM and logistic regression. Our aim is to show the performance of our proposed SPPA algorithms compared to other state of the art optimization algorithms using two different datasets, Bank Note Authentication and IMDB. The inner loop in SPPA algorithm is run at most 20 times for the results presented below.
Bank Note Authentication Dataset. The bank note authentication dataset is a publicly available dataset in UCI machine learning repository Dua & Graff (2017). Data consists of genuine and counterfeit banknote images. The dataset has 1372 instances. Each instance has 5 attributes out of which 4 are the features and one is the target attribute. The target attribute contains two values: 0 and 1, where 0 represents genuine note and 1 represents fake note. The ratio of the two classes are balanced 55/45 (genuine/counterfeit). We pre-process the data to have target values to be -1 and 1. We perform classification using SVM and logistic regression. Considering the per iteration complexity of SPPA, we plotted the loss-time, accuracy-time figures, Figure 1 and Figure 2 respectively. The Figure 1 demonstrates how the regularized hinge loss value decreases with respect to time (seconds) using different algorithms. The Figures 1 and 2 show that with our efficient implementation, SPPA takes less amount of time to reach to a smaller loss, and a greater accuracy despite of its per iteration complexity. The range shown in the plots clearly shows that SPPA is more robust to change in stochastic settings.Due to lack of space we present further experiments on the performance of SPPA using logistic regression on Bank Note Authentication dataset in supplementary.
IMDB Dataset. IMDB is large movie review dataset, used for binary sentiment classification Maas et al. (2011). There are 25,000 movie reviews for training, and 25,000 movie reviews for testing. We used the processed bag of words format to run our experiments. The ratio of the two classes are balanced, 50/50 (positive sentiment/negative sentiment). The target attribute contains two values: −1 and 1, where −1 represents negative sentiment whereas 1 represents positive sentiment. Similar to Bank Note Authentication dataset, we perform classification using SVM and logistic regression with IMDB dataset. SPPA outperforms the other algorithms when we use logistic regression to do the classification. Surprisingly, SGD Momentum diverged with this particular data set for logistic regression problem. Even though SGD Momentum performs similar to our proposed SPPA in most cases, this experiment shows that it is not as stable as SPPA. We tried different parameter settings for all algorithms and the presented results are the best performance of each algorithm. The Figure 3 shows the logistic loss, and the Figure 4 shows the accuracy per seconds. Focusing on the Figure 3, we see that the logistic loss follows a decreasing trend with SPPA whereas with other algorithms the decrease is not too visible. Moreover, Figure 4 shows that SPPA reaches to a higher test accuracy compared to the other algorithms.The performance of SPPA is similar, but not particularly better than the state of the art methods when we use SVM. Related experiment results are available in supplementary.

heading: 3.2 DEEP NEURAL NET CLASSIFICATION
text: CIFAR10. CIFAR10 (Krizhevsky et al., 2010) consists of 50,000 training and 10,000 test images of size 32 × 32 in 10 classes. The network we use for this experiment is based on convolutional deep neural network. It consists of 5 layers, each of the first 2 being a combination of 5 × 5 convolutional filters and 2 × 2 max pooling and last 3 being fully connected. The loss function is cross entropy loss, and the activation function is RELU. Batch size is 4. To show the behaviour of the algorithms more in detail we recorded the loss values at every 500 stochastic update, the experiment is run for 5 epochs and in total 125 updates are presented (25 × 500 × 4 giving the number of training samples). The accuracy is calculated on the test data at every 500 stochastic update overall for 5 epochs, using the accuracy calculation in PyTorch tutorial. Convergence plots are shown in Figures 5 and 6. To show that SPPA works not only on the CNN architecture, we also tried training a residual neural network (ResNet) on CIFAR10. For the details about the network, one can refer to 20 layer ResNet architecture on CIFAR10 (He et al., 2016). In our settings the batch size is 32. The results are shown in Figures 7 and 8. As we can see, the performance is indeed consistent as in the CNN case.
MNIST. MNIST Handwritten Digits Data set (LeCun et al., 2010) is a common data set used for showing the effectiveness of many state of art the algorithms. It consists of 60,000 training, 10,000 test images of size 28 × 28 in 10 classes. The network we used for this experiment is based on convolutional deep neural network. It consists of 5 layers, each of the first 3 being combination of 3×3 convolutional filters and 2×2 max pooling with stride of 2 and last 2 being fully connected. The loss function is cross entropy loss, and the activation function is RELU. Batch size is 64. To show the behaviour of the algorithms more in detail we recorded the loss values at every 25 stochastic update, the experiment is run for 1 epoch and in total 38 updates are presented (25×38×64 giving the number of training samples). The accuracy is calculated on the test data at every 500 stochastic update, using the accuracy calculation in PyTorch tutorial. 2 As observed in Figures 9 and 10, SPPA-based methods
2https://github.com/pytorch/tutorials/blob/master/beginner_source/ blitz/cifar10_tutorial.py
outperform all SGD-based algorithms in terms of both the cross entropy loss and prediction accuracy on the test set.

heading: 4 CONCLUSION
text: In this paper we presented efficient implementations of the stochastic proximal point algorithm (SPPA) for general ERM problems and selected deep learning problems. We showed that for each of the ERM problem, one iteration of SPPA can be expressed in closed form or even as a simple optimization algorithm that can be easily solved by bisection method. We also proved that with our proposed efficient implementation, SPPA has same order of complexity as SGD. There do exist some convergence results of SPPA for convex problems different than ours, but our approach require only convexity of loss functions. The new convergence analysis also shows great resemblance to its deterministic counterpart, which is not the case from other works. The resulting algorithm is more robust in different stochastic settings and has cheap per-iteration complexity, while enjoying convergence under milder conditions. We demonstrated the performance of SPPA with our proposed efficient implementation on some well-known classification and regression data sets for ERM problems, and showed that they outperform many existing methods.

heading: A EFFICIENT IMPLEMENTATION
text: In this section we introduce several efficient methods to calculate the proximal operator update in Algorithm 1 line 3. At first glance it may seem as hard as solving the batch problem itself, but we will see that there are some interesting properties when 58 involve the loss evaluation of only one or a few data points in a data fitting scenario. On the other hand, it is not clear how to parallelize SPPA with mini-batches, which leaves room for future work.

heading: A.1 LEAST SQUARES LOSS
text: For the least squares loss 58 = (H8 − x>8w)2, the proximal update is a linear least squares problem with closed-form solution wC+1 = (O + _Cx8x>8 )−1 (wC + x8H8). It can be efficiently calculated by invoking the Sherman-Morrison formula and avoid directly inverting a 3 × 3 matrix for w ∈ R3 as
wC+1 = wC − x> 8 wC − H8
‖x8 ‖2 + 1/_ x8 .
The update rule looks like, but is not exactly the same as, the normalized least mean squares (NLMS) algorithm (Haykin, 2002). It has appeared in (Pătraşcu, 2020), but we include it here for completeness.

heading: A.2 LOGISTIC (CROSS-ENTROPY) LOSS
text: In general, the one-sample loss function can be written as 58 (w) = x>8w = ℓ(x>8w). For example, the one-sample loss for logistic regression is 58 (w) = log(1 + exp(−H8x>8w)). Suppose it is differentiable, then by applying the chain rule, we have ∇ 58 (w) = x8ℓ ′ (x> 8 w), where ℓ′ is the derivative of the scalar function ℓ(·). Plugging this into the optimality condition _C∇ 58 (wC+1) + wC+1 − wC = 0, we see that wC+1 has the form wC+1 = wC − Ux8 , where U is some scalar. The proximal update reduces to finding the scalar U that solves the following convex problem,
minimize U
ℓ(x>8wC − U‖x8 ‖2) + U2
2_C ‖x8 ‖2.
This can be done by setting the derivative equal to zero. Due to convexity, its derivative is a monotonic function, so the root can be obtained via bisection.
As a concrete example, consider the logistic loss 58 (w) = log(1 + exp(−H8x)8 w)) with H8 = ±1. According to the aforementioned arguments, we need to solve the nonlinear equation
U = H8_C exp(x> 8 wC − U‖x8 ‖2)
1 + exp(x> 8 wC − U‖x8 ‖2)
.
Notice that the right-hand-side is a number between 0 and ±_C , which gives the initial upper and lowerbound on U. Furthermore, we see that x>
8 wC and ‖x8 ‖2 only need to be calculated once, with
O(3) flops for w ∈ R3; each bisection step takes constant time and it needs no more than 20 ∼ 30 scalar computations to render an accurate-enough solution.

heading: A.3 HINGE LOSS
text: Regarding nonsmooth optimization problems, it turns out many of the widely used loss functions admit closed form updates. The main idea is to consider the subgradient calculus, and find the point where 0 is in the subdifferential. Take support vector machine (SVM) as an example, in which the loss function is the hinge loss 58 (w) = [1 − H8x>8w]+. Its subdifferential is
m 58 (w) =  {0} 1 − H8x>8w < 0, {−H8x8} 1 − H8x>8w > 0, {−UH8x8 |0 ≤ U ≤ 1} 1 − H8x>8w = 0.
When added with a proximal term, the subdifferential set is added with w − wC . As a result, the update rule is simply
wC+1 =  wC H8x > 8 wC > 1, wC + _H8x8 H8x>8wC < 1 − _ (C) ‖x8 ‖2, wC + _H8x8 1−H8x>8wC H8 ‖x8 ‖2 otherwise.
An interesting observation here is that it looks like the perceptron algorithm with an adaptively chosen step-size.

heading: A.4 ABSOLUTE ERROR LOSS
text: Robust regression using absolute error loss 58 (w) = |H8 − x>8w |, is another example to nonsmooth convex loss function. The derivation is similar to that of SVM: the subdifferential for the one-sample loss is
m 58 (w) =  {x8} H8 − x>8w < 0, {−x8} H8 − x>8w > 0, {Ux8 | − 1 ≤ U ≤ 1} H8 − x>8w = 0.
After adding the proximal term, the update rule is
wC+1 =  wC − _x8 H8 − x>8wC < −_‖x8 ‖2, wC + _x8 H8 − x>8wC > _‖x8 ‖2, wC + H8−x>8wC ‖x8 ‖2 x8 otherwise.

heading: A.5 GENERIC LOSSES
text: Even for generic nonlinear programming problems, it is still possible to efficiently evaluate the proximal update beyond merely a simple gradient step. The idea is to apply the limited-memory BFGS algorithm (Nocedal & Wright, 2006), or L-BFGS for short. L-BFGS is a memory-efficient implementation of the famous quasi-Newton algorithm BFGS. In a nut shell, L-BFGS is an iterative algorithm that makes use of the second-order information from the optimization loss function, but does not require solving matrix inverses and only requires explicitly evaluating first-order gradients. For a prescribed number of iterations, it requires one matrix-vector multiplication and multiple vector multiplications. As a result, the overall complexity is again O(3) if the initial guess of the Hessian matrix is diagonal.
While there are certain limitations for applying L-BFGS to general nonlinear programming problems, we reckon that it fits perfectly in the context of SPPA implementations.
• L-BFGS has to specify a good initial guess of the Hessian approximation matrix. While in many cases people simply use the identity matrix to start, it may result in very poor approximation. Fortunately, thanks to the proximal term, the identity matrix is in fact a very good initial guess for the Hessian matrix for SPPA updates.
• In order to save memory consumption, L-BFGS has to prescribe the number of iterations before running the algorithm. Obviouly, if the prescribed number of iteration is too large, we incur unnecessary computations, while if it is too small we need to invoke another round with a new estimated Hessian matrix. However, again in the context of SPPA, the proximal term naturally provides a good initialization wC . Our experience show that prescribing 10 iterations of L-BFGS updates is more than enough to obtain accurate solutions.
On the other hand, it perhaps makes more sense to simply use some more advanced first-order methods such as Nesterov’s accelerated gradient descent (1983; 2013) to calculate the proximal update. Both L-BFGS and accelerated gradient descent evaluates the gradient of the loss function with O(3) complexity, and the question is how to leverage convergence rate versus sophistication.

heading: B PROOF OF LEMMA 1
text: At every iteration of SPPA, a component function 58 is picked uniformly at random to obtain the update. Therefore ?(8C |wC ) = 1/=, obviously. Lemma 1, however, states that the conditional
probability of 8C given the next update wC+1 is still uniform. Recall that the update rule is
wC+1 = argmin w _C 58C (w) + (1/2)‖w − wC ‖2. (10)
What Lemma 1 implies is that if we know wC+1 is obtained at iteration C, without knowing the previous iterate wC , we gain no additional information regarding which component function 58 is more likely to have rendered the update equation 10.
Consider the sequence of indexes 81, 82, . . .. At every iteration of SPPA, an index is independently drawn from {1, . . . , =} with equal probabilities. We may also treat this as a Markov chain with the transition probability ?(8C |8C−1) = 1/=. This seems unnecessary at the moment, since an independent sequence is a trivial special case of a Markov chain. However, this implies ?(8C |wC ) = ?(8C |8C−1) and subsequently ?(8C |wC+1) = ?(8C |8C+1). The question then boils down to studying the transition probability of the reversed Markov chain.
Kolmogorov’s criterion gives a necessary and sufficient condition for a Markov chain to be identical to its time-reversed version (Kelly, 2011). Theorem 3. A stationary Markov chain is reversible if and only if its transition probability satisfies
?( 91 | 92)?( 92 | 93) · · · ?( 9<−1 | 9<)?( 9< | 91) = ?( 91 | 9<)?( 92 | 91) · · · ?( 9< | 9<−1) for any finite sequence of states 91, . . . , 9<.
It is easy to see that in our case the Kolmogorov’s criterion trivially holds since the transition probability is 1/= in all cases. This proves that ?(8C |8C+1) = ?(8C |8C−1) = 1/= and thus ?(8C |wC+1) = 1/=.

heading: C ALMOST SURE CONVERGENCE
text: In the main paper we presented the result that the infimum of E[ (wC )] goes to zero, following the techniques presented by Bottou et al. (2018). What this implies is that for some iteration B we may have a very small E[ (wB)], but for C > B the expected loss may go up again. To make a stronger claim that the random variable (wC ) converges almost surely, we are going to use martingales convergence results. The difference in this case is that we are going to use the backwards version of martingale convergence. Here we list the key results referenced from (Pollard, 2002). Definition 1. Let {-C : C ∈ N} be a sequence of integrable random variables, adapted to a decreasing filtration {GC : C ∈ N}. We call {(-C ,GC ) : C ∈ N} a reversed supermartingale if E[-C |GC+1] ≥ -C+1 for all C.
In other words, {(-C ,GC ) : C ∈ N} is a reversed supermartingale if and only if {(-−C ,G−C ) : C ∈ −N} is a supermartingale. In particular, for each fixed C, the finite sequence -C , -C−1, . . . , -0 is a supermartingale with respect to the filtration GC ⊆ GC−1 ⊆ · · · ⊆ G0. The theory for reversed positive supermartingale is analogous to the regular supermartingale theory. Theorem 4. For every reversed, positive supermartingale {(-C ,GC ) : C ∈ N}, there exists an -∞ for which -C → -∞ almost surely.
Equipped with the reversed supermartingale theory, we have the following: Theorem 5. If all 51, . . . , 5= are convex and _C is lowerbounded, then the sequence {wC } generated by SPPA (Algorithm 1) satisfies that (wC ) → (w★) almost surely.
Proof. Define the filtration GC = {wC , wC+1, . . .} for all C = 0, 1, . . .. Consider inequality equation 2 from Proposition 1. Notice that
E[‖wC − w★‖2 | wC+1] = E[‖wC − w★‖2 | GC+1] . This means {(‖wC − w★‖2,GC ) : C ∈ N} is a reversed supermartingale. Applying Theorem 4, we have that {‖wC − w★‖2} converges almost surely; in other words ‖wC − w★‖2 − ‖wC+1 − w★‖2 → 0 almost surely.
Notice that Proposition 1 also implies
‖wC − w★‖2 − ‖wC+1 − w★‖2 ≥ 2_C ( (wC+1) − (w★)). Combining this with the previous result shows that
_C ( (wC+1) − (w★)) → 0 almost surely. Since we assume that _C is lowerbounded for all C, this means (wC+1) → (w★) almost surely.

heading: D ADDITIONAL EXPERIMENTS

heading: D.1 CLASSIFICATION
text: In the main paper, we showed the performance of our unique implementation of SVM on Bank Note Authentication dataset and Logistic Regression on IMDB dataset. For the sake of completeness, we will show the performance of SPPA with our implementation for Logistic Regression on Bank Note Authentication dataset. Figure 11 shows the logistic loss value per second using Bank Note Authentication dataset. The performance of SPPA with our implementation is comparable with SGD with momentum. It still performs better than most of the other algorithms and it is more robust to change in seeds. Figure 12 shows that SPPA with our implementation reaches to the highest accuracy value on test data in less than a second, whereas the other algorithms take around 3 seconds to reach to an even smaller accuracy value. The robustness of SPPA with our implementation is consistent in this experiment too.

heading: D.1.1 CLASSIFICATION WITH L1 REGULARIZATION
text: We perform binary classification using SVM and logistic regression with L1 regularization. Our aim is to show the performance of our proposed SPPA algorithms in presence of L1 regularization term compared to other state of the art optimization algorithms using two different datasets, Bank Note Authentication and IMDB. In the results presented below, for each algorithm every stochastic update is followed by soft-thresholding to impose the effect of L1 regularization term. The inner loop in SPPA algorithm is run at most 20 times for the results presented below.
Bank Note Authentication Dataset. We perform classification using SVM and logistic regression with L1 regularization. Considering the per iteration complexity of SPPA, we plotted the loss-time, accuracy-time figures, Figure 13 and Figure 14 respectively. The Figure 13 demonstrates how the regularized hinge loss value decreases with respect to time (seconds) using different algorithms. The Figures 13 and 14 show that with our efficient implementation, SPPA takes less amount of time to reach to a smaller loss, and a greater accuracy despite of its per iteration complexity. The range shown in the plots clearly shows that SPPA is more robust to change in stochastic settings.
We present further experiments on the performance of SPPA using logistic regression with L1 regularization on Bank Note Authentication dataset loss-time, accuracy-time figures, Figure 15 and
Figure 16 respectively. In Figure 15 SPPA reaches to a smaller loss value in shorter amount of time compared to other algorithms, the closest performance is by Stochastic Gradient with Momentum algorithm. In Figure 16 SPPA reaches to the largest accuracy value in shortest amount of time and it has smallest variance compared to the other algorithms.
IMDB Dataset. We perform classification using SVM and logistic regression with L1 regularization on IMDB dataset. SPPA outperforms the other algorithms when we use logistic regression to do the classification. We aim to see the change in performances of the algorithms in presence of L1 regularization. According to the experiments presented there is not a big difference in the performance of algorithms with L1 regularization except for the slightly increase on the loss values. The Figure 17 shows the loss, and the Figure 18 shows the accuracy per seconds when logistic regression model is used with different optimization algorithms to perform classification on IMDB dataset.

heading: D.2 REGRESSION
text: Bike Sharing. Bike Sharing is a publicly available data set containing hourly or daily count of rental bikes as well as environmental and seasonal settings. The data set we are using is from UC Irvine Machine Learning Repository. 3 The core data set is related to the two-year historical log corresponding to years 2011 and 2012 from Capital Bike share system, Washington D.C., USA (Fanaee-T & Gama, 2013). It consists of 17379 samples (12512 training, 3476 test, 1391 validation)
3https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset
with 16 features each to predict hourly count of rental bikes. We perform linear regression with two different loss functions, MSE loss. The Figure 19 shows mean squared error loss per seconds. We used 100 seeds to test the robustness of algorithms. SPPA reaches to a smaller loss value in less number of iterations. Batch size is 16.
Hotel Average Daily Rates (ADR). Hotel ADR is a publicly available data set consisting of hotel demand data for two different types of hotels; resort hotels and city hotels. Each hotel data set consists of 40,060 samples with 31 features each to predict average daily rate values (Antonio et al., 2019). In our experiments, we used resort hotels with 40,060 samples (30045 training, 10015 test) with 8 features. We perform linear regression with mean squared error (MSE). The Figure 20 shows mean squared error loss per seconds. SPPA reaches to a smaller loss value in less number of iterations and in less time. Batch size is 16. In Figure 20, MSE loss goes to very small value only in first few iterations. Our experiments show that overall performance of SPPA is better in linear regression problems using Hotel ADR dataset.

heading: D.3 REGRESSION WITH L1 REGULARIZATION
text: We perform linear regression with L1 regularization and show the results of different algorithms including our efficient implementation on both Bike Sharing 4 and Hotel Average Daily Rates (ADR) datasets (Antonio et al., 2019).
Bike Sharing. We perform linear regression with L1 regularization on Bike sharing dataset and compare the performances of different optimization algorithms. The Figure 21 shows the decrease in the loss values over time for different algorithms. SPPA outperforms the other algorithms and reaches to a smaller value in a shorter time. We tried different settings for the algorithms, in the results presented in Figure 21 the learning rate for SPPA is 1/B@AC (C) where C is the number of iterations.
4https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset
The other advantage, stability of SPPA is also showcased in Figure 21. Compared to the results without the L1 regularizer, the L1 regularizer helps SPPA to get smaller loss values.
Hotel Average Daily Rates (ADR). We perform linear regression with L1 regularization on Hotel ADR dataset and compare the performances of different optimization algorithms. The Figrure 22 shows the decrease in loss values over time for different algorithms. SPPA converges to a very small loss value much faster, in less than a second it converges to the minimum loss. We tried different settings for the algorithms, in the results presented in Figure 22 the learning rate for SPPA is 1/C where C is the number of iterations.

heading: D.4 MEAN-ABSOLUTE-ERROR REGRESSION
text: We showcased the performance of our algorithm with mean squared error loss in regression problem with two real data-sets in the main paper. Here, we will present our results on the same data sets using absolute error loss. Figure 23 shows that, SPPA with our unique implementation reaches to a smaller loss value in less amount of time and significantly outperforms the alternative methods on Bike Sharing dataset. However, the error range of the algorithms are similar. Figure 24 shows the performance on Hotel ADR dataset with absolute error loss, our implementation of SPPA reaches to a smaller loss value in a significantly less amount of time, in the first few iterations. In addition to its great performance, looking at the error bar one can see that it is much more robust to different seed settings compared to other algorithms.

Review summaries:
# === Update ===  I have read the authors reviews as well as the response by the authors.  I'd like to thank the authors for engaging in a discussion around Lemma 1. It is unfortunate that the result does not appear to hold even in a restricted setting.  I will maintain my score. I wish the authors the best in their future work.  # ======  This submission investigates the convergence rate of the stochastic proximal point algorithm (SPPA) for convex, finite-sum functions. The authors derive new convergence rates for SPPA under the assumption that each function in the finite sum is convex (or strongly convex). Key to these derivations is a lemma characterizing the probabilities of sampling a sub-function form a time-reversible Markov chain.  In addition to these results, the authors derive efficient implementations of SPPA for a variety of loss functions, generally by reducing the problem to one-dimensional root-finding. The paper concludes with experiments comparing SPPA to SGD, Adam, and other baseline stochastic optimization algorithms for (i) convex optimization and (ii) training CNNs and ResNets on MNIST and CIFAR-10, respectively. The experiments show that SPPA converge faster than stochastic gradient methods and have a similar computational cost.

This paper provides a novel elegant proof of convergence for the Stochastic Proximal Point Algorithm, under the assumption of convexity. It proceeds by surprisingly conditioning on future iterates, reversing usual convergence proofs. Experiments show that the SPPA is competitive, including in nonconvex settings such as training neural networks for image classification.

This paper shows the arbitrary convergence rate of the stochastic proximal point algorithm(SPPA), under the convex assumption of the objective function. The authors use a novel approach to analyze SPPA. Numerical results show the efficiency of SPPA on some realistic data-sets.