Extended Regularized Dual Averaging Methods for Stochastic Optimization

TL;DR

This paper proposes the extended regularized dual averaging (XRDA) algorithm for stochastic optimization, offering flexible backward step size control to improve convergence while maintaining sparsity and theoretical guarantees.

Contribution

It introduces XRDA, a novel extension of RDA that allows bounded backward step sizes, enhancing convergence and practical performance in stochastic optimization.

Findings

XRDA achieves comparable convergence rates to RDA.

Flexible step size control improves convergence speed.

The method preserves sparsity in solutions.

Abstract

We introduce a new algorithm, extended regularized dual averaging (XRDA), for solving regularized stochastic optimization problems, which generalizes the regularized dual averaging (RDA) method. The main novelty of the method is that it allows a flexible control of the backward step size. For instance, the backward step size used in RDA grows without bound, while for XRDA the backward step size can be kept bounded. We demonstrate experimentally that additional control over the backward step size can significantly improve the convergence rate of the algorithm while preserving desired properties of the iterates, such as sparsity. Theoretically, we show that the XRDA method achieves the same convergence rate as RDA for general convex objectives.

Tables

Table 1. Table 2: Demonstrates how particular parameter choices in the XRDA algorithm ( 1 ) correspond to SDA, RDA, and forward-backward gradient descent.

	${\mu }_n=0$	${\mu }_n=1$
${\alpha }_n=1$ $s_n=n^{-\frac{1}{2}}$	Forward-Backward Gradient descent $x_{n+\frac{1}{2}}=x_n-s_n{g}_n$ $x_{n+1}=\arg {\min }_{x\in V}\frac{1}{2}{\Vert x-x_{n+\frac{1}{2}}\Vert }_V^2+s_{n}G\left(x\right)$	Forward-Backward RDA $x_{n+\frac{1}{2}}=x_{n+1}-s_n{g}_n$ $x_{n+1}=\arg {\min }_{x\in V}\frac{1}{2}{\Vert x-x_{n+\frac{1}{2}}\Vert }_V^2+\left({\sum }_{i=1}^n{s}_n\right)G\left(x\right)$
${\alpha }_n=\sqrt{n}$ $s_n=1$	Simple Dual Averaging (SDA) introduced in [18] (for $G=0$)	Regularized Dual Averaging (RDA) introduced in [24]

Equations

$$x_{n+1}=x_n-s_n{g}_n,$$

$$x_{n+1}=\arg \underset{x}{\min }\left(\sum _{i=1}^n⟨s_i{g}_i,x⟩+\frac{{\alpha }_{n+1}}{2}\Vert x-x_1{\Vert }^2\right),$$

$$x_{n+1}=\frac{{\alpha }_n}{{\alpha }_{n+1}}{x}_n+\left(1-\frac{{\alpha }_n}{{\alpha }_{n+1}}\right)x_1-{\tilde{s}}_n{g}_n,$$

$$\arg \underset{x\in A}{\min }[f(x)=F(x)+G(x)],$$

$$\begin{array}{rl}& x_{n+\frac{1}{2}}=x_n-s_n{g}_n\\ & x_{n+1}-x_{n+\frac{1}{2}}\in -s_n\partial G(x_{n+1}).\end{array}$$

$$\begin{array}{rl}& x_{n+\frac{1}{2}}=x_n-s_n{\tilde{g}}_i\\ & x_{n+1}=\arg \underset{x}{\min }\left(\lambda s_n\Vert x{\Vert }_1+\frac{1}{2}\Vert x-x_{n+\frac{1}{2}}{\Vert }^2\right)=\{\begin{array}{ll}{x}_{n+\frac{1}{2}}-\lambda s_n& x>\lambda s_n\\ 0& |x|\le \lambda s_n\\ x_{n+\frac{1}{2}}+\lambda s_n& x<-\lambda s_n.\end{array}\end{array}$$

$$x_{n+1}=\arg \underset{x}{\min }\left(\sum _{i=1}^n⟨s_i{g}_i,x⟩+\frac{{\alpha }_{n+1}}{2}\Vert x-x_1{\Vert }^2+{\gamma }_{n}G(x)\right),$$

$$\begin{array}{rl}& x_{n+\frac{1}{2}}=x_{n-\frac{1}{2}}-s_n{\tilde{g}}_i\\ & x_{n+1}=\arg \underset{x}{\min }\left({\gamma }_{n}G(x)+\frac{1}{2}\Vert x-x_{n+\frac{1}{2}}{\Vert }^2\right),\end{array}$$

$$\begin{array}{rl}& x_{n+\frac{1}{2}}=\left(1-{\mu }_n\right)x_{n-\frac{1}{2}}+{\mu }_n{x}_n-s_n{\tilde{g}}_i\\ & x_{n+1}=\arg \underset{x}{\min }\left(\lambda {\gamma }_n|x{|}_1+\frac{1}{2}\Vert x-x_{n+\frac{1}{2}}{\Vert }^2\right)\end{array}$$

$$\nabla f\in V^{\ast },\partial f\subset V^{\ast }.$$

$$\varphi (\alpha x+(1-\alpha )y)\le \alpha \varphi (x)+(1-\alpha )\varphi (y)-\frac{\sigma }{2}\Vert x-y{\Vert }_V^2.$$

$$(\partial \varphi {)}^{-1}(h)=\arg \underset{x\in V}{\min }[\varphi (x)-⟨h,x⟩].$$

$$h_{n+1}=h_n-s_n{g}_n,x_{n+1}=(\partial \varphi {)}^{-1}(h_{n+1}),$$

$$D_{\varphi ,h}(x,y)=\varphi (x)-\varphi (y)-⟨h,x-y⟩,$$

$$x_{n+1}=\arg \underset{x\in V}{\min }\left(s_n⟨g_n,x⟩+D_{\varphi ,h_n}(x,x_n)\right),h_{n+1}=h_n-s_n{g}_n.$$

$$\varphi (x)=\{\begin{array}{ll}\sum _{i=1}^n{x}_i\log (x_i)& x\in A\\ \infty & \text{otherwise}\end{array}$$

$$(\partial \varphi {)}^{-1}(x{)}_i={\left(\sum _{j=1}^n{e}^{x_j}\right)}^{-1}{e}^{x_i},$$

$$D_{\varphi }(x,y)=\varphi (x)-\varphi (y)-⟨\nabla \varphi ,x-y⟩,$$

$$x_{n+1}=\arg \underset{x\in V}{\min }\left(s_n⟨g_n,x⟩+D_{\varphi }(x,x_n)\right),$$

$$\nabla \varphi (x_{n+1})=\nabla \varphi (x_n)-s_n{g}_{n+1},$$

$$x_{n+1}=\arg \underset{x\in V}{\min }(s_{n}G(x)+D_{\varphi }(x,x_n)).$$

$$G(x)=i_C(x)=\{\begin{array}{ll}0& \text{if x∈C}\\ +\infty & \text{if x∈/C},\end{array}$$

$$l_i(z)=f(x_i)+⟨g_i,z-x_i⟩.$$

$${\bar{x}}_n=S_n^{-1}\sum _{i=1}^n{s}_i{x}_i,$$

$$f({\bar{x}}_n)\le S_n^{-1}\sum _{i=1}^n{s}_{i}f(x_i),$$

$$f(z)\ge S_n^{-1}\sum _{i=1}^n{s}_i{l}_i(z)=S_n^{-1}\sum _{i=1}^n{s}_i(f(x_i)+⟨g_i,z-x_i⟩).$$

$$f({\bar{x}}_n)-f(z)\le S_n^{-1}\sum _{i=1}^n{s}_i(f(x_i)-l_i(z))=S_n^{-1}\sum _{i=1}^n{s}_i⟨g_i,x_i-z⟩$$

$$x_{n+1}=\arg \underset{x}{\min }\left(\sum _{i=1}^n{s}_i{l}_i(x)+{\alpha }_{n+1}{D}_{\varphi }(x,x_1)\right)$$

$$\sum _{i=1}^n⟨g_i,x_i-z⟩\le \left({\alpha }_n{D}_{\varphi }(z,x_1)+\frac{M^2}{2\sigma }\sum _{i=1}^n\frac{s_i^2}{{\alpha }_i}\right),$$

$$\arg \underset{x\in A}{\min }[f(x)=F(x)+G(x)].$$

Full text

K****eywords Convex Optimization  $\cdot$ Subgradient Methods  $\cdot$ Structured Optimization  $\cdot$ Non-smooth Optimization

1 Introduction

Optimizing convex objectives with a regularization term which promotes a certain structure of the minimizer, for example the $\ell^{1}$-norm promoting sparsity or the nuclear norm promoting low rank [8], are very common and important in machine learning. Stochastic optimization of such objectives poses a unique challenge since traditional methods such as stochastic gradient descent (SGD) [14] are often not effective at producing the desired structure (i.e. sparsity or low-rank) of the iterates [24]. This motivated the introduction of the well-known RDA method and its variants, which are able to effectively produce the desired structure of the iterates [24, 9]. We introduce a generalization of RDA, called extended regularized dual averaging (XRDA), which we show can significantly improve convergence while still preserving the desired structure of the iterates.

Let us begin by describing these issues in detail and giving an overview of RDA. Consider the subgradient descent and dual averaging methods [18] for minimizing a Lipschitz convex function $F$, given by

[TABLE]

where $g_{n}\in\partial F(x_{n})$. It is well known that with a step size $s_{n}=n^{-\frac{1}{2}}$ this method attains a convergence rate of $O(n^{-\frac{1}{2}}\log{n})$. The simple dual averaging (SDA) method of Nesterov [18], which is given by

[TABLE]

generalizes subgradient descent (note that setting $\alpha_{n}=1$ recovers the iteration (1)). The advantage of the more general SDA method is that by setting $s_{n}=1$ and $\alpha_{n}=\sqrt{n}$, the logarithmic factor in the convergence rate can be removed [18]. To illustrate the relationship between this new method and original subgradient descent, we rewrite this new iteration as

[TABLE]

where the effective step size is $\tilde{s}_{n}=\frac{s_{n}}{\alpha_{n+1}}$. Thus the difference between this new method and the subgradient descent method is an averaging with the initial iterate $x_{1}$. It is remarkable that this provides a significant improvement in the convergence rate. These results hold in more generality with the $\ell^{2}$ distance replaced by the Bregman distance $D_{\phi}(x,y)$ with respect to a convex function $\phi$, which we describe in more detail in Section 2, but for simplicity we stay with the current setting throughout the introduction.

Next, we consider the composite (or regularized) optimization problem

[TABLE]

where $F(x)$ is a convex Lipschitz function and $G(x)$ is a convex function. A standard method for solving this is the forward-backward subgradient method

[TABLE]

where $g_{n}\in\partial F(x_{n})$ [3]. Note that the second step above corresponds to backward Euler and is known as the proximal map for $G$ [6]. With a choice of step size $s_{n}=O(n^{-\frac{1}{2}})$, this method also achieves a convergence rate of $O(n^{-\frac{1}{2}}\log{n})$. As before the logarithm can be removed by introducing a similar averaging with $x_{1}$ as in the SDA method. Note also that the constants in the convergence rates only depend upon the Lipschitz constant of $F$ and not $G$, which is the advantage of using the forward-backward splitting.

In many cases of practical interest, the subgradients $g_{i}\in\partial F(x_{n})$ in the forward step are not computed exactly, but rather replaced by an unbiased sample $\tilde{g}_{i}$. Using this sample in (5) results in forward-backward stochastic gradient descent. Stochastic gradient descent has proven extremely useful for training a variety of machine learning models [14, 15]. However, for problems where the minimizer is expected to have a special structure, forward-backward stochastic gradient descent often has the drawback that the iterates it produces do not have the desired structure. A very common example is sparsity. For instance, consider the forward-backward stochastic gradient descent algorithm applied to an objective $F(x)+G(x)$ with $G(x)=\lambda\|x\|_{1}$ an $l^{1}$ regularization term:

[TABLE]

Here $\tilde{g}_{i}$ denotes a sampled subgradient. Notice that the soft-thresholding parameter in the above iteration is $\lambda s_{n}$, which must be very small since the step size $s_{n}$ must go to [math] if we wish to converge to the exact optimizer [24] (note here that the gradient samples are stochastic). This means that the iterates generated by this methods will not be sparse even if the true optimizer is sparse. Further, due to the slow convergence rate of the algorithm, thresholding after training may introduce significant errors even though the true optimizer is sparse [24]. This property of stochastic gradient descent also applies to other types of structure, such as low rank structure expected when using the nuclear norm as a regularizer [8].

This motivated the introduction of the regularized dual averaging (RDA) method [24], which is obtained by modifying the SDA algorithm in the following way

[TABLE]

where $\gamma_{n}=\sum_{i=1}^{n}s_{i}$. We can clarify the relationship between this method and the forward-backward subgradient descent method (5) by setting $\alpha_{n}=1$ and rewriting it in the following ‘leap-frog’ form

[TABLE]

where $\gamma_{n}=\sum_{i=1}^{n}s_{i}$. This method achieves a convergence rate of $O(n^{-\frac{1}{2}}\log{n})$ and the logarithm here can be removed by introducing an averaging with $x_{1}$ similar to SDA (this corresponds to a different choice of the parameters $s_{n}$ and $\alpha_{n}$). Notice that in the RDA iteration (8), the backward step size grows as the sum of the previous forward step sizes $s_{n}$ and grows without bound as $n\rightarrow\infty$. This is in contrast to forward-backward stochastic gradient descent (5), where the backward step and forward step are the same, and thus the backward step size goes to [math]. As a result, the RDA method produces sparse iterates while forward backward SGD does not. However, it is quite unintuitive that the backward step size should grow without bound, which we attempt to address in this paper.

We introduce the extended regularized dual averaging (XRDA) method, which generalize the RDA method. This method gives more precise control over the backward step size, for instance a special case of XRDA is the ‘averaged leap-frog scheme’

[TABLE]

where $0\leq\mu_{i}\leq 1$ and $\gamma_{n}=(1-\mu_{n})\gamma_{n-1}+s_{n}$. As we show in Section 4, the convergence rate of this algorithm is $O(n^{-\frac{1}{2}}\log{n})$ and a slightly modified version achieves the optimal black-box convergence rate of $O(n^{-\frac{1}{2}})$ [17, 19, 24]. Further, in Section 5 we show experimentally that choosing the parameter $\mu_{n}$ so that the backward step $\gamma_{n}$ remains bounded leads to significantly improved convergence over RDA while preserving the sparsity of the iterates.

Concerning non-convex applications of the XRDA method, we note that the BinaryConnect [10] and BinaryRelax [25] methods, which have been used to train quantized neural networks, are simply the XRDA method applied to non-convex objectives. Specifically, BinaryConnect consists of the leap frog iteration (8) applied to the case where $G=\chi_{Q}$ is the indicator function of the set $Q$ of quantized parameters, while BinaryRelax is just (8) applied to $G(x)=\text{dist}(x,Q)^{2}$. In addition, the RDA method has also been successfully applied to the training of sparse neural networks [12] and recent numerical experiments show that the XRDA method improves upon this for the training of sparse neural networks [22].

The paper is organized as follows. In Section 2, we introduce some basic notation and the basic problem formulation. In Section 3, we recall the basic ideas behind the dual averaging [18] and RDA [24] methods. Then, in Section 4, we introduce the extended regularized dual averaging (XRDA) method and derive its convergence rate. Here we also exhibit a variety of special cases of the XRDA method and compare with RDA and SDA. Then, in Section 5 we present experimental results showing that the XRDA method significantly exhibits significantly improved convergence relative to RDA while still preserving the sparsity of the iterates. Finally, we give concluding remarks.

2 Notation and Setup

In this section, we introduce the basic setting and notation that we will use. This section is essentially a review of the basics of mirror descent [16] and is a common framework for subgradient methods in convex optimization [7]. We refer to [21] for the basic notions of relative interior, lower semi-continuity, and subgradients of convex functions.

We are interested in first-order methods for optimizing convex functions $f$ on a closed convex subset $A\subset V$ of a Banach space $V$. A universal issue, which applies to all first-order methods, including gradient descent, subgradient descent, accelerated gradient descent, and conjugate gradient, is that gradients and subgradients of the objective lie naturally in the dual space $V^{*}$, i.e.

[TABLE]

Consequently, it does not a priori make sense to add gradients of $f$ to primal points in $V$. In order to circumvent this, we will introduce a convex lower semicontinuous function $\phi:V\rightarrow\mathbb{R}\cup\{\infty\}$, called a mirror function (see [7, 16, 1, 17]). We require that $\phi(x)=\infty$ for $x\notin A$, and that $\phi$ is strongly convex with parameter $\sigma$, i.e. that for $x,y\in V$ and $\alpha\in[0,1]$, we have

[TABLE]

Note that this is only non-trivial if $x,y\in A$. The subdifferential of the mirror function $\phi$ gives a map $(\partial\phi)^{-1}:V^{*}\rightarrow A$, defined by

[TABLE]

Our assumptions on $\phi$ imply that the above minimizer will lie in $A$. Further, the strong convexity of $\phi$ implies that the minimizer above is unique. From the objective (12), it is also clear that $x=(\partial\phi)^{-1}(h)$ iff $h\in\partial\phi(x)$.

The idea of mirror descent is to perform gradient descent in the dual space and to use $(\partial\phi)^{-1}$ to move iterates from the dual to the primal space. Specifically, we let $h_{0}\in V^{*}$ and set $x_{0}=(\partial\phi)^{-1}(h_{0})$. The mirror gradient step [16, 1, 17] with step size $s_{n}$ is given by

[TABLE]

where $g_{n}\in\partial f(x_{n})$. This iteration clearly preserves the property that $h_{n}\in\partial\phi(x_{n})$ at every step.

The iteration (13) is commonly written using the Bregman distance [5]

[TABLE]

defined for $x,y\in V$ and $h\in\partial\phi(x)$, as follows

[TABLE]

This formulation will also be important in our later analysis.

Let us consider a few common examples to illustrate this abstract framework. Suppose that $V$ is a Hilbert space and $A=V$. In this case, a simple choice for $\phi$ is the squared norm $\phi(x)=\frac{1}{2}\|x\|_{V}^{2}$. The map $(\partial\phi)^{-1}$ in this case is the identity and the mirror descent iteration (13) reduces to regular gradient descent.

A common variation on this example is to choose the mirror function $\phi(x)=\frac{1}{2}\langle x,Ax\rangle_{V}$ where $A$ is a positive definite operator on the Hilbert space $V$. In this case, the map $(\partial\phi)^{-1}$ is given by $(\partial\phi)^{-1}(x)=A^{-1}x$, and mirror descent (13) corresponds to preconditioned gradient descent with (linear) preconditioner $A^{-1}$.

Finally, another common special case[13, 17, 2, 7, 4] is when $V=\ell^{1}_{n}$ and $A=\{x\geq 0:x^{T}1=1\}$ is the probability simplex. In this case, one chooses

[TABLE]

as the negative entropy, which is strongly convex with respect to the $|\cdot|_{1}$ norm. This fact is known as Pinsker’s inequality [20]. In this case, the map $(\partial\phi)^{-1}$ is given by

[TABLE]

which is called the softmax function. See [7] and the references therein for more examples.

In the following, we will need to consider a special case of this setup where we additionally assume that $\phi$ restricted to the affine hull $\text{aff}(A)$ is differentiable on the relative interior $\text{relint}(A)$ (see [21, 3, 26], for instance, for the definition and basic properties of the affine hull and relative interior). In addition, we will need that $\nabla\phi$ (here we are restricting to the affine hull) satisfies $\nabla\phi(x)\rightarrow\infty$ as $x\rightarrow\partial_{ri}A$ approaches the relative boundary of $A$. Note that all of the examples given here satisfy these stronger assumptions.

In fact, in this case it is convenient to make the additional assumption that $\text{aff}(A)=V$. This can be assumed without loss of generality by restricting to $\text{aff}(A)$. These additional assumptions imply that $\phi$ is differentiable on $A^{o}$ and that $\nabla\phi:A^{o}\rightarrow V^{*}$ and $(\nabla\phi)^{-1}:V^{*}\rightarrow A^{o}$ (here $A^{o}$ denotes the interior of $A$) gives inverse bijections between $V^{*}$ and $A^{o}$. Additionally, the Bregman distance is given by

[TABLE]

for $x,y\in A^{o}$ in the interior of $A$ (see [16, 7] for these general properties and theory). The mirror descent iteration (15) simplifies to

[TABLE]

where $g_{n}\in\partial f(x_{n})$. Equivalently, we can write $\nabla\phi(x_{n+1})=\nabla\phi(x_{n})-s_{n}g_{n}$.

The iteration (17) corresponds to a forward gradient step. Analogously, there is a corresponding backward gradient step for an objective $G(x)$. This corresponds to solving the backward equation

[TABLE]

where $g_{n+1}\in\partial G(x_{n+1})$. This is easily seen to be equivalent to

[TABLE]

For practical implementation, it is important that the optimization problem appearing in (19) can be efficiently solved. We note that if $\phi=\frac{1}{2}\|\cdot\|_{2}^{2}$, then two common choices for $f$ are

[TABLE]

in which case the mirror-prox update (19) is a projection onto the set $C$, and $G(x)=|x|_{1}$, in which case the mirror-prox (19) update is the well-known soft-thresholding step [23].

In the following, we will use combinations of the forward (17) and backward steps (19) to optimize a composite objective $f(x)=F(x)+G(x)$.

3 Dual Averaging and Regularized Dual Averaging

We begin by briefly recalling the dual averaging method [18] for minimizing a convex objective $f$ to help explain the idea behind both RDA [24] and the XRDA method.

The idea behind the dual averaging method is to observe that, due to the convexity of the objective $f$, the sequence of iterates $x_{i}$ and subgradients $g_{i}\in\partial f(x_{i})$ generated by our algorithm leads to a sequence of lower bounds on the objective

[TABLE]

Now let $s_{i}$ be a sequence of weights and consider the primal average

[TABLE]

where $S_{n}=\sum_{i=1}^{n}s_{i}$. From the convexity of $f$, we obtain a bound on the objective at $\bar{x}_{n}$

[TABLE]

while the corresponding dual average (i.e. average of the lower bounds $l_{i}$) is a lower bound on $f$, i.e.

[TABLE]

Combining equations (23) and (24) we see that

[TABLE]

holds for all $z$, in particular it holds for $z=x^{*}$. The challenge now is to choose the sequence of iterates $x_{i}$ and weights $s_{i}$ so that the right hand side of equation (25) goes to [math] as $n\rightarrow\infty$ for $z=x^{*}$. The solution to this problem presented in the seminal paper [18] by Nesterov is to choose

[TABLE]

for a sequence of parameters $s_{i}$ and $\alpha_{i}$, in which case we have

[TABLE]

where $M$ is the Lipschitz constant of $f$ (i.e. $\|g_{i}\|_{V^{*}}\leq M$). The choice $\alpha_{i}=1$ and $s_{i}=i^{-\frac{1}{2}}$, which corresponds to subgradient descent, gives a convergence rate which is $O(n^{-\frac{1}{2}}\log{n})$. Optimal choices of the weights $s_{i}$ and parameters $\alpha_{i}$ allow one to remove the logarithmic factor and obtain an objective convergence rate which is $O(n^{-\frac{1}{2}})$ [18].

Next, we wish to generalize the dual averaging method to composite objectives $f$ which are the sum of two pieces

[TABLE]

Here $F$ is a convex function which is Lipschitz with respect to the norm $\|\cdot\|_{V}$ and $G$ is a convex function for which we can solve the backward step (19). As discussed in the previous section, common choices for $G$ include, for example, the $l_{1}$-norm (if $\phi=\frac{1}{2}\|x\|_{2}^{2}$) or the characteristic function of a convex set.

In his seminal work [24], Xiao extended Nesterov’s analysis to the composite case by using the lower bound

[TABLE]

where $g_{i}\in\partial F(x_{i})$, in (26). So here we are using a linear lower bound on the piece $F$ and using the function $G$ itself in the lower bound. This results in the RDA method [24]

[TABLE]

where $\gamma_{n+1}=\sum_{i=1}^{n}s_{i}$. From the properties of the Bregman distance discussed in the preceding section, the optimizer above satisfies

[TABLE]

where $\tilde{g}_{n+1}\in\partial G(x_{n+1})$. Introducing the intermediate iterate $x_{n+\frac{1}{2}}$ defined by

[TABLE]

we can rewrite (30) in terms of a forward and backward step as

[TABLE]

The advantages of this method are twofold. First, the convergence behavior only depends upon the Lipschitz constant of the function $F$. Second, the backward step $\gamma_{n}$ grows, which allows the method to obtain sparse iterates even with large step sizes and noisy gradients [24].

4 Extended Regularized Dual Averaging

In this section, we slightly generalize the RDA method to obtain the XRDA method. To do this, we consider the more general lower bound

[TABLE]

where $g_{i}\in\partial F(x_{i})$, $h_{i}\in\partial G(x_{i})$, and $q_{i}\geq 0$ is a parameter. This is easily seen to be a lower bound since $h_{i}\in\partial G(x_{i})$ so that $G(x_{i})+\langle h_{i},z-x_{i}\rangle\leq G(z)$. Note that this lower bound is not even necessarily convex if $q_{i}>1$, but remarkably we can still obtain convergence of the resulting method. We plug this lower bound into equation (26), to get the method (setting $t_{i}=s_{i}q_{i}$)

[TABLE]

Here $\gamma_{n+1}=\sum_{i=1}^{n}(1-q_{i})s_{i}=\sum_{i=1}^{n}(s_{i}-t_{i})$ (restrictions introduced later on the parameters will guarantee that $\gamma_{n+1}>0$). If we introduce the new parameter $\mu_{i}=t_{i}/\gamma_{i}$, we obtain the perhaps more enlightening relation $\gamma_{n+1}=(1-\mu_{n})\gamma_{n}+s_{n}$.

Of course, this is only useful if we can determine an $h_{i}\in\partial G(x_{i})$. We note that the optimality condition corresponding to the iteration (34), combined with the simple form of the Bregman distance (16), implies that

[TABLE]

provides an element in $\partial G(x_{n+1})$. Furthermore, if we assume that $x_{1}\in\arg\min_{x}G(x)$, then we may choose $h_{1}=0$. Utilizing this choice of $h_{n}$ in (34), we obtain the XRDA method, which we have summarized in forward form in Table (1).

An illustrative special case of the algorithm can be obtained by setting $\alpha_{n}=1$, which gives

[TABLE]

Note that the last line above is a backward proximal step and can also be written as

[TABLE]

This choice results in a method which is analogous to subgradient descent, and achieves a convergence rate which is $O(n^{-\frac{1}{2}}\log{n})$. Removing the logarithmic factor is possible by taking a more sophisticated choice of the parameter $\alpha_{n}$, analogous to dual averaging [18]. Finally, in the special case where $A=V=\mathbb{R}^{n}$ and $\phi(x)=\frac{1}{2}\|x\|_{2}^{2}$, we get the iteration

[TABLE]

Since $\gamma_{n+1}=(1-\mu_{n})\gamma_{n}+s_{n}$, we can control the backward stepsize $\gamma_{n}$ by choosing the averaging parameter $\mu_{n}$ appropriately. We give a full comparison of the XRDA method with forward-backward gradient descent, the SDA method and the RDA method in table 2. Here we can see how different parameter choices result in different special cases of the method.

Similar to the dual averaging and RDA methods, the XRDA method is robust to noisy gradients. This means that in more generality, we can let $g_{i}$ be independent random variables whose expectation satisfy $\mathbb{E}(g_{i})\in\partial F(x_{i})$.

Specifically, let $(\Sigma,\mathcal{F},\mathcal{P})$ denote a probability space and define the subgradient sample function (dependent on the function $F$) $g:V\times\Sigma\rightarrow V^{*}$. Given a point $x\in V$ and a random sample $\xi\in\Sigma$, $g$ returns an unbiased sample of an element in the subgradient of the objective $F$ at $x$. Mathematically, this means that we require

[TABLE]

Using this in place of $g_{i}$ in algorithm (1), we obtain a stochastic version of the XRDA method. We have the following convergence theorem.

Theorem 1**.**

Assume that $F$ is convex and the gradient samples are bounded by $M$, i.e. that $\|g(x,\xi)\|_{V^{*}}\leq M$. Let $x^{*}\in\arg\min_{x}f(x)$. Let $x_{n}$ be determined by algorithm (1) with $0\leq t_{i}\leq\gamma_{i}$, $\alpha_{i}\leq\alpha_{i+1}$, and $s_{i+1}\leq s_{i}$. Then, if $x_{1}\in\arg\min_{x}G(x)$ and $h_{1}=0$ (i.e. $x_{\frac{1}{2}}=x_{1}$), we have

[TABLE]

where $S_{n}=\sum_{i=1}^{n}s_{i}$ and $\bar{x}_{n}=S_{n}^{-1}\sum_{i=1}^{n}s_{i}x_{i}$ is a weighted average of the iterates with weights $s_{i}$. We also have

[TABLE]

Notice that $\gamma_{1}$ is the empty sum and so is [math]. This means since $0\leq t_{1}\leq\gamma_{1}$ that $t_{1}=0$ as well. In addition, the condition that $0\leq t_{1}\leq\gamma_{1}$ together with the relation $\gamma_{n+1}=\sum_{i=1}^{n}(s_{i}-t_{i})$ imply that $\gamma_{i}\geq s_{i-1}>0$.

Proof.

First consider the case where the subgradients are not sampled stochastically, but are deterministic. We begin by noting that by convexity we have

[TABLE]

and since the minimum is smaller than the average we have

[TABLE]

Thus it suffices to prove that

[TABLE]

For this we first use the sub-gradient property (and the assumption that $t_{i}\geq 0$) to get

[TABLE]

and rewrite the first sum above as

[TABLE]

We proceed by bounding the term

[TABLE]

Recall from (34) that $x_{i}=\arg\min_{x}f_{i}(x)$ where

[TABLE]

The definition of $h_{i}\in\partial G(x_{i})$ given in (35) (note that $i=1$ is not a problem since $x_{1}\in\arg\min_{x}G(x)$ and $h_{1}=0$) means that

[TABLE]

So as long as $t_{i}\leq\gamma_{i}$, $f_{i}^{\prime}$, defined by

[TABLE]

is a strongly convex function with convexity parameter $\alpha_{i}\sigma$ (since $\phi$ is strongly convex with parameter $\sigma$), and $0\in\partial f^{\prime}_{i}(x_{i})$. This implies that (see Lemma 6 in [18])

[TABLE]

Adding $\langle s_{i}g_{i},x_{i}-z\rangle$ to both sides gives

[TABLE]

Hölder’s inequality, combined with the Cauchy-Schwartz inequality implies that

[TABLE]

and we get (since $g_{i}\in\partial F(x_{i})$ and so the Lipschitz assumption on $F$ implies $\|g_{i}\|_{V^{*}}\leq M$)

[TABLE]

Expanding the definition of $f^{\prime}_{i}$, we obtain

[TABLE]

Plugging this bound into (46) with $z=x_{i+1}$ and $z=x^{*}$, (recalling that $\gamma_{n}=\sum_{i=1}^{n-1}s_{i}-t_{i}$ and $\gamma_{1}=t_{1}=0$), we obtain

[TABLE]

Plugging this into (45) and rearranging terms, we obtain

[TABLE]

Finally, we use the assumption that $\alpha_{i-1}-\alpha_{i}\leq 0$, that $s_{i-1}-s_{i}\geq 0$, and that $x_{1}\in\arg\min_{x}G(x)$ to conclude that

[TABLE]

and that

[TABLE]

This yields

[TABLE]

which completes the proof in the deterministic case.

Now suppose that the subgradients are stochastic. Taking the expectation of equations (42) and (43), we see that it suffices to prove that

[TABLE]

Using the assumption that $\mathbb{E}_{\xi_{i}}(g(x_{i},\xi_{i}))\in\partial f(x_{i})$, the fact that $x_{i}$ only depends on $\xi_{1},...,\xi_{i-1}$, and the assumption that $t_{i}\geq 0$, we see that

[TABLE]

Since $h_{i}$ and $x_{i}$ are independent of $\xi_{i}$, we can pull the $\mathbb{E}_{\xi_{i}}$ out of the inner product above and the left hand side simply becomes

[TABLE]

We now bound the term inside the expectation exactly as before to obtain

[TABLE]

Combining this with (60) and (59), we get

[TABLE]

which completes the proof.

∎

5 Experimental Results

In this section, we provide experimental evidence that the XRDA method exhibits improved convergence behavior over RDA on some convex optimization problems. In particular, we consider the problem of calculating a sparse logistic regression on the MNIST dataset [14]. Specifically, the objective we consider is

[TABLE]

where $W\in\mathbb{R}^{k\times n}$, $b\in\mathbb{R}^{k}$ are the parameters of the logistic regression, $x_{i}\in\mathbb{R}^{n}$ is the training data and $y_{i}\in\mathbb{R}^{k}$ are the training labels [11]. For MNIST the specific parameters are $N=60000$ (number of training data points), $n=784$ (dimension of the image data), $k=10$ (number of classes).

We optimize the objective (63) with $\lambda=5\cdot 10^{-4}$ using stochastic forward backward-splitting, which is given by

[TABLE]

Here $\theta_{n}=(W_{n},b_{n})$ are the parameters to be optimized. Next, we optimize the same objective using forward-backward RDA, given by

[TABLE]

Finally, we optimize the objective using XRDA which is set to asymptotically produce a fixed backward step-size of $M$. This iteration is given by

[TABLE]

Note that this corresponds to the XRDA iteration 1 with $\alpha_{n}=1$, $\mu_{n}=1-s_{n}/M$ and $\phi=\frac{1}{2}\|\cdot\|^{2}$. We take $M=500,1000,2500,5000,$ and $10000$. We run all of these methods with an initial iterate $\theta_{1}=\theta_{1/2}=0$, a batch size of $10$, and step size $s_{n}=3.0\cdot n^{-\frac{1}{2}}$ for a total of $50$ epochs. Note that this scaling of the step size results in the optimal theoretical convergence rate of $O(n^{-\frac{1}{2}})$ for all methods. We plot the number of non-zero components, the test accuracy and the training loss of each of these methods as a function of iteration in Figure 1.

From this experiment, we draw the following conclusions. First, we see that forward-backward stochastic gradient descent doesn’t generate sparse solutions, as expected. In contrast, RDA and XRDA all generate sparse solutions. Further, the larger the backward step size $M$ used in XRDA, the sparser the iterates are. This extends to RDA which generates the sparsest iterates. However, the difference between the sparsity generated by these methods is not particularly large.

Next, considering the training loss (i.e. the objective (63) we are trying to minimize), we see that the larger backward step size used in XRDA improves the convergence rate over forward-backward SGD. In particular, the loss converges fastest with asymptotic backward step size $M=500$ and $M=1000$. With an asymptotic backward step size of $M=2500$, XRDA converges at the roughly the same rate as forward-backward SGD initially and then converges faster for larger iterations. However, once the backward step size of XRDA is pushed beyond this point, the convergence rate slows dramatically as can be seen for $M=5000$ and $M=10000$. Finally, RDA exhibits by far the slowest convergence rate. This shows that the large backward step size produced by the RDA iteration can dramatically slow down convergence. The best results are obtained using XRDA with parameters chosen to give an appropriate backward step size.

6 Conclusion

We have derived and analysed the extended regularized dual averaging (XRDA) method in Table 1, a class of methods for solving regularized optimization problems which generalize RDA [24]. We derived and analysed the convergence of these methods for regularized convex optimization problems. The main novelty of these new methods is that we can control the proximal (or backward) step size. Finally, we have demonstrated experimentally that controlling the backward step size can significantly improve the convergence of the resulting methods while still preserving the desired structure (such as sparsity) of the iterates.

We also remark that these methods are useful for non-convex optimization problems appearing in machine learning as well, since BinaryConnect [10] and BinaryRelax [25] are special cases of the method applied to non-convex neural networks trianing, RDA has successfully been applied to non-convex problems such as the training of sparse neural networks [12], and the XRDA method has also been shown to effectively train sparse neural networks [22].

7 Acknowledgements

We would like to thank Dr. Xiao Lin for helpful comments and discussion during the preparation of this work. This work was supported by the Verne M. Willaman Fund, the Center for Computational Mathematics and Applications (CCMA) at the Pennsylvania State University, and the National Science Foundation (Grant No. DMS-1819157).