A Simple Practical Accelerated Method for Finite Sums

Introduction

A large body of recent developments in optimization have focused on minimization of convex finite sums of the form:

a very general class of problems including the empirical risk minimization (ERM) framework as a special case. Any function $h$ can be written in this form by setting $f_{1}(x)=h(x)$ and $f_{i}=0$ for $i\neq 1$, however when each $f_{i}$ is sufficiently regular in a way that can be made precise, it is possible to optimize such sums more efficiently than by treating them as black box functions.

In most cases recently developed methods such as SAG (Schmidt et al., 2013) can find an $\epsilon$-minimum faster than either stochastic gradient descent or accelerated black-box approaches, both in theory and in practice. We call this class of methods fast incremental gradient methods (FIG).

FIG methods are randomized methods similar to SGD, however unlike SGD they are able to achieve linear convergence rates under Lipschitz-smooth and strong convexity conditions (Mairal, 2014; Defazio et al., 2014b; Johnson and Zhang, 2013; Konečný and Richtárik, 2013). The linear rate in the first wave of FIG methods directly depended on the condition number ($L/\mu$) of the problem, whereas recently several methods have been developed that depend on the square-root of the condition number (Lan and Zhou, 2015; Lin et al., 2015; Shalev-Shwartz and Zhang, 2013c; Nitanda, 2014). Analogous to the black-box case, these methods are known as accelerated methods.

In this work we develop another accelerated method, which is conceptually simpler and requires less tuning than existing accelerated methods. The method we give is a primal approach, however it makes use of a proximal operator oracle for each $f_{i}$ instead of a gradient oracle, unlike other primal approaches. The proximal operator is also used by dual methods such as some variants of SDCA (Shalev-Shwartz and Zhang, 2013a).

Algorithm

Our algorithm’s main step makes use of the proximal operator for a randomly chosen $f_{i}$. For convenience, we use the following compact notation:

This proximal operator can be computed efficiently or in closed form in many cases, see Section 4 for details. Like SAGA, we also maintain a table of gradients $g_{i}$, one for each function $f_{i}$. We denote the state of $g_{i}$ at the end of step $k$ by $g_{i}^{k}$. The iterate (our guess at the solution) at the end of step $k$ is denoted $x^{k}.$ The starting iterate $x^{0}$ may be chosen arbitrarily.

The full algorithm is given as Algorithm 1. The sum of gradients $\frac{1}{n}\sum_{i=1}^{n}g_{i}^{k}$ can be cached and updated efficiently at each step, and in most cases instead of storing a full vector for each $g_{i}$, only a single real value needs to be stored. This is the case for linear regression or binary classification with logistic loss or hinge loss, in precisely the same way as for standard SAGA. A discussion of further implementation details is given in Section 4.

the expected convergence rate in terms of squared distance to the solution is given by:

as $\epsilon\rightarrow 0$. This rate matches the lower bound known for this problem (Lan and Zhou, 2015) under the gradient oracle. We conjecture that this rate is optimal under the proximal operator oracle as well. Unlike other accelerated approaches though, we have only a single tunable parameter (the step size $\gamma$), and the algorithm doesn’t need knowledge of $L$ or $\mu$ except for their appearance in the step size.

Compared to the $O\left(\left(L/\mu+n\right)\log\left(1/\epsilon\right)\right)$ rate for SAGA and other non-accelerated FIG methods, accelerated FIG methods are significantly faster when $n$ is small compared to $L/\mu$, however for $n\geq L/\mu$ the performance is essentially the same. All known FIG methods hit a kind of wall at $n\approx L/\mu$, where they decrease the error at each step by no more than $1-\frac{1}{n}$. Indeed, when $n\geq L/\mu$ the problem is so well conditioned so as to be easy for any FIG method to solve it efficiently. This is sometimes called the big data setting (Defazio et al., 2014b).

Our convergence rate can also be compared to that of optimal first-order black box methods, which have rates of the form $k=O\left(\left(\sqrt{L/\mu}\right)\log\left(1/\epsilon\right)\right)$ per epoch equivalent. We are able to achieve a $\sqrt{n}$ speedup on a per-epoch basis, for $n$ not too large. Of course, all of the mentioned rates are significantly better than the $O\left(\left(L/\mu\right)\log\left(1/\epsilon\right)\right)$ rate of gradient descent.

For non-smooth but strongly convex problems, we prove a $1/\epsilon$-type rate under a standard iterate averaging scheme. This rate does not require the use of decreasing step sizes, so our algorithm requires less tuning than other primal approaches on non-smooth problems.

Relation to other approaches

Our method is most closely related to the SAGA method. To make the relation clear, we may write our method’s main step as:

The difference is the point at which the gradient of $f_{j}$ is evaluated at. The proximal operator has the effect of evaluating the gradient at $x^{k+1}$ instead of $x^{k}$. While a small difference on the surface, this change has profound effects. It allows the method to be applied directly to non-smooth problems using fixed step sizes, a property not shared by SAGA or other primal FIG methods. Additionally, it allows for much larger step sizes to be used, which is why the method is able to achieve an accelerated rate.

It is also illustrative to look at how the methods behave at $n=1$. SAGA degenerates into regular gradient descent, whereas our method becomes the proximal-point method (Rockafellar, 1976):

The proximal point method has quite remarkable properties. For strongly convex problems, it converges for any $\gamma>0$ at a linear rate. The downside being the inherent difficulty of evaluating the proximal operator. For the $n=2$ case, if each term is an indicator function for a convex set, our algorithm matches Dykstra’s projection algorithm if we take $\gamma=2$ and use cyclic instead of random steps.

Several acceleration schemes have been recently developed as extensions of non-accelerated FIG methods. The earliest approach developed was the ASDCA algorithm (Shalev-Shwartz and Zhang, 2013b, c). The general approach of applying the proximal-point method as the outer-loop of a double-loop scheme has been dubbed the Catalyst algorithm Lin et al. (2015). It can be applied to accelerate any FIG method. Recently a very interesting primal-dual approach has been proposed by Lan and Zhou (2015). All of the prior accelerated methods are significantly more complex than the approach we propose, and have more complex proofs.

Theory

In this section we rehash some simple bounds from proximal operator theory that we will use in this work. Define the short-hand $p_{\gamma f}(x)=\text{prox}_{\gamma f}(x)$, and let $g_{\gamma f}(x)=\frac{1}{\gamma}\left(x-p_{\gamma f}(x)\right)$, so that $p_{\gamma f}(x)=x-\gamma g_{\gamma f}(x)$. Note that $g_{\gamma f}(x)$ is a subgradient of $f$ at the point $p_{\gamma f}(x)$. This relation is known as the optimality condition of the proximal operator. Note that proofs for the following two propositions are in the supplementary material.

In operator theory this property is known as $(1+\mu\gamma)$-cocoerciveness of $p_{\gamma f}$.

Recall our definition of $g_{\gamma f}(x)=\frac{1}{\gamma}\left(x-p_{\gamma f}(x)\right)$ also. After combining, the following relation thus holds between the proximal operator of the conjugate $f^{*}$ and $g_{\gamma f}$:

We will apply cocoerciveness of the proximal operator of $f^{*}$ as it appears in the decomposition. Note that L-smoothness of $f$ implies $1/L$-strong convexity of $f^{*}$. In particular we apply it to the points $\frac{1}{\gamma}x$ and $\frac{1}{\gamma}y$:

Pulling $\frac{1}{\gamma}$ from the right side of the inner product out, and plugging in Equation 9, gives the result. ∎

2 Notation

Let $x^{*}$ be the unique minimizer (due to strong convexity) of $f$. In addition to the notation used in the description of the algorithm, we also fix a set of subgradients $g_{j}^{*}$, one for each of $f_{j}$ at $x^{*}$, chosen such that $\sum_{j=1}^{n}g_{j}^{*}=0$. We also define $v_{j}=x^{*}+\gamma g_{j}^{*}.$ Note that at the solution $x^{*}$, we want to apply a proximal step for component $j$ of the form:

(Technical lemma needed by main proof) Under Algorithm 1, taking the expectation over the random choice of $j$, conditioning on $x^{k}$ and each $g_{i}^{k}$, allows us to bound the following inner product at step $k$:

The proof is in the supplementary material.

3 Main result

(single step Lyapunov descent) We define the Lyapunov function $T^{k}$ of our algorithm (Point-SAGA) at step $k$ as:

for $c=1/\mu L$. Then using step size $\gamma=\frac{\sqrt{(n-1)^{2}+4n\frac{L}{\mu}}}{2Ln}-\frac{1-\frac{1}{n}}{2L}$, the expectation of $T^{k+1}$, over the random choice of $j$, conditioning on $x^{k}$ and each $g_{i}^{k}$, is:

Term 1 of $T^{k+1}$ is straight-forward to simplify:

For term $2$ of $T^{k+1}$ we start by applying cocoerciveness (Theorem 11):

where we have pulled out the quadratic term by using $E[z_{j}^{k}-v_{j}]=x^{k}-x^{*}$ (we can take the expectation since the left hand side of the inner product doesn’t depend on $j$). We now expand $E\left\langle x^{k+1}-x^{k}\,,\,z_{j}^{k}-v_{j}\right\rangle$ further:

We further split the left side of the inner product to give two separate inner products:

The first inner product in Equation 4 is the quantity we bounded in Lemma 8 by $\gamma^{2}\frac{1}{n}\sum_{i=1}^{n}\left\|g_{i}^{k}-g_{i}^{*}\right\|^{2}$. The second inner product in Equation 4, can be simplified using Theorem 3 (note the right side of the inner product is equal to $z_{j}^{k}-v_{j}$):

Combing these gives the following bound on $(1+\mu\gamma)E\left\|x^{k+1}-x^{*}\right\|^{2}$:

Define $\alpha=\frac{1}{1+\mu\gamma}=1-\kappa$, where $\kappa=\frac{\mu\gamma}{1+\mu\gamma}$. Now we multiply the above inequality through by $\alpha$ and combine with the rest of the Lyapunov function, giving:

We want an $\alpha$ convergence rate, so we pull out the required terms:

Now to complete the proof we note that $c=1/\mu L$ and $\gamma=\frac{\sqrt{(n-1)^{2}+4n\frac{L}{\mu}}}{2Ln}-\frac{1-\frac{1}{n}}{2L}$ ensure that both terms inside the round brackets are non-positive, giving $ET^{k+1}\leq\alpha T^{k}$. These constants were found by equating the equations in the brackets to zero, and solving with respect to the two unknowns, $\gamma$ and $c$. It is easy to verify that $\gamma$ is always positive, as a consequence of the condition number $L/\mu$ always being at least 1.∎

(Smooth case) Chaining Theorem 5 gives a convergence rate for Point-SAGA at step $k$ under the constants given in Theorem 5 of:

where $\bar{x}^{k}=\frac{1}{k}E\sum_{t=1}^{k}x^{t}.$ The proof of this theorem is included in the supplementary material.

Implementation

Care must be taken for efficient implementation, particularly in the sparse gradient case. We discuss the key points below. A fast Cython implementation is available on the author’s website incorporating these techniques.

For the most common binary classification and regression methods, implementing the proximal operator is straight-forward. We include details of the computation of the proximal operators for the hinge, square and logistic losses in the supplementary material. The logistic loss does not have a closed form proximal operator, however it may be computed very efficiently in practice using Newton’s method on a 1D subproblem. For problems of a non-trivial dimensionality the cost of the dot products in the main step is much greater than the cost of the proximal operator evaluation. We also detail how to handle a quadratic regularizer within each term’s prox operator, which has a closed form in terms of the unregularized prox operator.

Instead of setting $g_{i}^{0}=f_{i}^{\prime}(x^{0})$ before commencing the algorithm, we recommend using $g_{i}^{0}=0$ instead. This avoids the cost of a initial pass over the data. In practical effect this is similar to the SDCA initialization of each dual variable to 0.

Experiments

We tested our algorithm which we call Point-SAGA against SAGA (Defazio et al., 2014a), SDCA (Shalev-Shwartz and Zhang, 2013a), Pegasos/SGD (Shalev-Shwartz et al., 2011) and the catalyst acceleration scheme (Lin et al., 2015). SDCA was chosen as the inner algorithm for the catalyst scheme as it doesn’t require a step-size, making it the most practical of the variants. Catalyst applied to SDCA is essentially the same algorithm as proposed in Shalev-Shwartz and Zhang (2013c). A single inner epoch was used for each SDCA invocation. Accelerated MISO as well as the primal-dual FIG method (Lan and Zhou, 2015) were excluded as we wanted to test on sparse problems and they are not designed to take advantage of sparsity. The step-size parameter for each method ($\kappa$ for catalyst-SDCA) was chosen using a grid search of powers of $2$. The step size that gives the lowest error at the final epoch is used for each method.

We selected a set of commonly used datasets from the LIBSVM repository (Chang and Lin, 2011). The pre-scaled versions were used when available. Logistic regression with $L_{2}$ regularization was applied to each problem. The $L_{2}$ regularization constant for each problem was set by hand to ensure $f$ was not in the big data regime $n\geq L/\mu$; as noted above, all the methods perform essentially the same when $n\geq L/\mu$. The constant used is noted beneath each plot. Open source code to exactly replicate the experimental results is available at https://github.com/adefazio/point-saga.

The key property that distinguishes accelerated FIG methods from their non-accelerated counterparts is their performance scaling with respect to the dataset size. For large datasets on well-conditioned problems we expect from the theory to see little difference between the methods. To this end, we ran experiments including versions of the datasets subsampled randomly without replacement in 10% and 5% increments, in order to show the scaling with $n$ empirically. The same amount of regularization was used for each subset.

Figure 1 shows the function value sub-optimality for each dataset-subset combination. We see that in general accelerated methods dominate the performance of their non-accelerated counter-parts. Both SDCA and SAGA are much slower on some datasets comparatively than others. For example, SDCA is very slow on the 5 and 10% COVTYPE datasets, whereas both SAGA and SDCA are much slower than the accelerated methods on the AUSTRALIAN dataset. These differences reflect known properties of the two methods. SAGA is able to adapt to inherent strong convexity while SDCA can be faster on very well-conditioned problems.

There is no clear winner between the two accelerated methods, each gives excellent results on each problem. The Pegasos (stochastic gradient descent) algorithm with its slower than linear rate is a clear loser on each problem, almost appearing as an almost horizontal line on the log scale of these plots.

Non-smooth problems

We also tested the RCV1 dataset on the hinge loss. In general we did not expect an accelerated rate for this problem, and indeed we observe that Point-SAGA is roughly as fast as SDCA across the different dataset sizes.

References

Appendix A Proximal operators

For the most common binary classification and regression methods, implementing the proximal operator is straight-forward. In this section let $y_{j}$ be the label or target for regression, and $X_{j}$ the data instance vector. We assume for binary classification that $y_{j}\in\{-1,1\}$.

The proximal operator has a closed form expression:

Logistic loss:

There is no closed form expression, however it can be computed very efficiently using Newton iteration, since it can be reduced to a 1D minimization problem. In particular, let $c_{0}=0$, $\gamma^{\prime}=\gamma\left\|X_{j}\right\|^{2}$, and $a=\left\langle z,X_{j}\right\rangle$. Then iterate until convergence:

The prox operator is then $\text{prox}_{\gamma f_{j}}(z)=z-\left(a-c^{k}\right)X_{j}/\left\|X_{j}\right\|^{2}$. Three iterations are generally enough, but ill-conditioned problems or large step sizes may require up to 12. Correct initialization is important, as it will diverge when initialized with a point on the opposite side of 0 from the solution.

Squared loss:

Let $\gamma^{\prime}=\gamma\left\|X_{j}\right\|^{2}$ and $a=\left\langle z,X_{j}\right\rangle$. Define:

Then $\text{prox}_{\gamma f_{j}}(z)=z-\left(a-c\right)X_{j}/\left\|X_{j}\right\|^{2}.$

L2 regularization

Including a regularizer within each $f_{i}$, i.e. $F_{i}(x)=f_{i}(x)+\frac{\mu}{2}\left\|x\right\|^{2},$ can be done using the proximal operator of $f_{i}$. Define the scaling factor:

Then $\text{prox}_{\gamma F_{i}}(z)=\text{prox}_{\rho\gamma f_{i}}(\rho z)$.

Appendix B Proofs

Under Algorithm 1, taking the expectation over the random choice of $j$, conditioning on $x^{k}$ and each $g_{i}^{k}$, allows us to bound the following inner product at step $k$:

We start by splitting on the right hand side of the inner product:

The first inner product has expectation on the left hand side (Recall that $E[g_{j}^{*}]=0$), so it’s simply 0 in expectation (we may take expectation on the left since the right doesn’t depend on $j$). The second inner product is the same on both sides, so we may convert it to a norm-squared term. So we have:

The inequality used is just an application of the variance formula $E[\left(X-E[X]\right)^{2}]=E[X^{2}]-E[X]^{2}\leq E[X^{2}].$∎

Chaining the main theorem gives a convergence rate for point-saga at step $k$ under the constants given in of:

First we simplify $T^{0}$ using $c=1/\mu L$ and use Lipschitz smoothness:

Now recall that the main theorem gives a bound $E\left[T^{k+1}\right]\leq\left(1-\kappa\right)T^{k}$ where the expectation is conditional on $x^{k}$ and each $g_{i}^{k}$ from step $k$, taking expectation over the randomness in the choice of $j$. We can further take expectation with respect to $x^{k}$ and each $g_{i}^{k}$, giving the unconditional bound:

where $\bar{x}^{k}=\frac{1}{k}E\sum_{t=1}^{k}x^{t}.$

Recall the bound on the Lyapunov function established in the main theorem:

In the non-smooth case this holds with $L=\infty$. In particular, if we take $c=\alpha\gamma^{2}n$, then:

Recall that this expectation is (implicitly) conditional on $x^{k}$ and each $g_{i}^{k}$ from step $k$, Taking expectation over the randomness in the choice of $j$. We can further take expectation with respect to $x^{k}$ and each $g_{i}^{k}$, and negate the inequality, giving the unconditional bound:

We can drop the $-E\left[T^{k}\right]$ since it is always negative. Dividing through by $k$:

Now using Jensen’s inequality on the left gives:

where $\bar{x}^{k}=\frac{1}{k}E\sum_{t=1}^{k}x^{t}.$ Now we plug in $T^{0}=\frac{c}{n}\sum_{i}\left\|g_{i}^{0}-g_{i}^{*}\right\|^{2}+\left\|x^{0}-x^{*}\right\|^{2}$ with $c=\alpha\gamma^{2}n\leq\gamma^{2}n$:

Now we plug in the bounds in terms of $B$ and $R$:

In order to balance the terms on the right, we need:

So we can take $\gamma=R/B\sqrt{n}$, giving a rate of:

Appendix C Proximal operator bounds with proofs

In this section we prove some simple bounds from proximal operator theory that we will use in this work. Define the short-hand $p_{\gamma f}(x)=\text{prox}_{\gamma f}(x)$, and let $g_{\gamma f}(x)=\frac{1}{\gamma}\left(x-p_{\gamma f}(x)\right)$, so that $p_{\gamma f}(x)=x-\gamma g_{\gamma f}(x)$. Note that $g_{\gamma f}(x)$ is a subgradient of $f$ at the point $p_{\gamma f}(x)$. This relation is known as the optimality condition of the proximal operator.

Using strong convexity of $f,$ we apply Equation 6 at the (sub-)gradients $g_{\gamma f}(x)$ and $g_{\gamma f}(y)$, and their corresponding points $p_{\gamma f}(x)$ and $p_{\gamma f}(y)$:

We now multiply both sides by $\gamma$, then add $\left\|p_{\gamma f}(x)-p_{\gamma f}(y)\right\|^{2}$ to both sides:

leading to the bound by using the optimality condition: $p_{\gamma f}(x)+\gamma g_{\gamma f}(x)=x.$∎

Recall our definition of $g_{\gamma f}(x)=\frac{1}{\gamma}\left(x-p_{\gamma f}(x)\right)$ also. After combining, the following relation thus holds between the proximal operator of the conjugate $f^{*}$ and $g_{\gamma f}$:

Let $u=p_{\gamma f}(x)$, and $v=\frac{1}{\gamma}\left(x-u\right)$. Then $v\in\partial f(u)$ by the optimality condition of the proximal operator of $f$ (namely if $u=p_{\gamma f}(x)$ then $u=x-\gamma v\Leftrightarrow v\in\partial f(u)$). It follows by conjugacy of $f$ that $u\in\partial f^{*}(v).$ Thus we may interpret $v=\frac{1}{\gamma}\left(x-u\right)$ as the optimality condition of a proximal operator of $f^{*}$ :

Plugging in the definition of $v$ then gives:

Further plugging in $u=p_{\gamma f}(x)$ and rearranging gives the result. ∎