Cocoercivity, Smoothness and Bias in Variance-Reduced Stochastic Gradient Methods

TL;DR

This paper introduces SVAG, a variance-reduced stochastic gradient method with adjustable bias, analyzing its convergence under cocoercivity and smoothness assumptions, and examining the practical effects of bias in classification tasks.

Contribution

It presents SVAG, a novel biased variance-reduced method, and provides convergence analysis highlighting the limitations of cocoercivity assumptions in biased settings.

Findings

SVAG's step-size requirements are less restrictive when operators are gradients.

Bias affects early convergence stages more than later stages.

Theoretical convergence conditions match or improve previous results, but differ from practical observations.

Abstract

With the purpose of examining biased updates in variance-reduced stochastic gradient methods, we introduce SVAG, a SAG/SAGA-like method with adjustable bias. SVAG is analyzed in a cocoercive root-finding setting, a setting which yields the same results as in the usual smooth convex optimization setting for the ordinary proximal-gradient method. We show that the same is not true for SVAG when biased updates are used. The step-size requirements for when the operators are gradients are significantly less restrictive compared to when they are not. This highlights the need to not rely solely on cocoercivity when analyzing variance-reduced methods meant for optimization. Our analysis either match or improve on previously known convergence conditions for SAG and SAGA. However, in the biased cases they still do not correspond well with practical experiences and we therefore examine the effect of bias numerically on a set of classification problems. The choice of bias seem to primarily affect the early stages of convergence and in most cases the differences vanish in the later stages of convergence. However, the effect of the bias choice is still significant in a couple of cases.

Figures

Equations

$$\underset{x\in {\mathbb{R}}^N}{\min }\frac{1}{n}\sum _{i=1}^n{f}_i(x)$$

$$x^{k+1}$$

$$y_{i^k}^{k+1}$$

$$y_j^{k+1}$$

$$x^{k+1}$$

$$0=Rx:=\frac{1}{n}\sum _{i=1}^n{R}_{i}x$$

$$R{x}^k\approx {\tilde{R}}^k=\frac{\theta }{n}(R_{i^k}{x}^k-y_{i^k}^k)+\frac{1}{n}\sum _{j=1}^n{y}_j^k.$$

$$E[{\tilde{R}}^k|x^k,y_1^k,\dots ,y_n^k]=\frac{\theta }{n}R{x}^k+(1-\frac{\theta }{n})\frac{1}{n}\sum _{j=1}^n{y}_j^k$$

$$(A\otimes B{)}^T=A^T\otimes B^T,(A\otimes B)(C\otimes D)=(AC)\otimes (BD).$$

$${\tau }_i{\mu }_j\text{ for all }i\in \{1,\dots ,m\},j\in \{1,\dots ,l\}$$

$$C_1\times C_2=\{(c_1,c_2)\mid c_1\in C_1,c_2\in C_2\}.$$

$$X^{\star }=\{x|0=\frac{1}{n}{\sum }_{i=1}^n{R}_{i}x\}$$

$$Z^{\star }=\{(x,R_{1}x,\dots ,R_{n}x)|0=\frac{1}{n}{\sum }_{i=1}^n{R}_{i}x\}.$$

$$⟨Rx-Ry,x-y⟩\ge \frac{1}{L}\Vert Rx-Ry{\Vert }^2$$

$$f(y)\le f(x)+⟨\nabla f(x),y-x⟩+\frac{L}{2}\Vert y-x{\Vert }^2$$

$$E[z^{k+1}|{\mathcal{F}}_k]\le (1+{\beta }^k)z^k+{\xi }^k-{\zeta }^k$$

$$z^{k+1}=z^k-(U_{i^k}\otimes I)\mathbf{B}{z}^k$$

$$U_i=\left[\begin{array}{cc}\frac{\lambda }{n}\theta e_i^T& -\frac{\lambda }{n}\theta e_i^T+\frac{\lambda }{n}{1}^T\\ -E_i& E_i\end{array}\right]$$

$$\forall z^{\star }\in Z^{\star },\forall i\in \{1,\dots ,n\}:0=(U_i\otimes I)\mathbf{B}{z}^{\star }.$$

$$0=⟨{\sum }_{i=1}^n{R}_{i}x-{\sum }_{i=1}^n{R}_{i}y,x-y⟩={\sum }_{i=1}^n⟨R_{i}x-R_{i}y,x-y⟩$$

$$\ge {\sum }_{i=1}^n\frac{1}{L}\Vert R_{i}x-R_{i}y{\Vert }^2\ge 0,$$

$$z(\omega )=(\bar{x},{\bar{y}}_1,\dots ,{\bar{y}}_n),z^k(\omega )=({\bar{x}}^k,{\bar{y}}_1^k,\dots ,{\bar{y}}_n^k)$$

$$\Vert \frac{1}{n}{\sum }_{i=1}^n{R}_i\bar{x}\Vert \le \Vert \frac{1}{n}{\sum }_{i=1}^n{R}_i{\bar{x}}^{k(l)}-\frac{1}{n}{\sum }_{i=1}^n{R}_i\bar{x}\Vert +\Vert \frac{1}{n}{\sum }_{i=1}^n{R}_i{\bar{x}}^{k(l)}\Vert$$

$$\le L\Vert {\bar{x}}^{k(l)}-\bar{x}\Vert +\Vert \frac{1}{n}{\sum }_{i=1}^n{R}_i{\bar{x}}^{k(l)}\Vert \to \Vert \frac{1}{n}{\sum }_{i=1}^n{R}_i{x}^{\star }\Vert =0$$

$$E[\Vert z^{k+1}-z^{\star }{\Vert }_{H\otimes I}^2|{\mathcal{F}}_k]$$

$$\le \Vert z^k-z^{\star }{\Vert }_{H\otimes I}^2-\Vert \mathbf{B}{z}^k-\mathbf{B}{z}^{\star }{\Vert }_{(2M-E[U_{i^k}^{T}H{U}_{i^k}]-\xi I)\otimes I}^2$$

$$-\xi nL⟨R{x}^k,x^k-x^{\star }⟩$$

$$H=\left[\begin{array}{cc}1& -\frac{\lambda }{n}(n-\theta )1^T\\ -\frac{\lambda }{n}(n-\theta )1& \frac{\lambda }{L}I+\frac{{\lambda }^2}{n^2}(n-\theta {)}^2{11}^T\end{array}\right]$$

$$M=\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right]\otimes \frac{1}{2n}\frac{\lambda }{L}I-\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\otimes \frac{{\lambda }^2}{2n^2}(n-\theta ){11}^T.$$

$$E[F((K\otimes I)z^{k+1})|{\mathcal{F}}_k]\le F((K\otimes I)z^k)-\Vert \mathbf{B}{z}^k-\mathbf{B}{z}^{\star }{\Vert }_{\frac{1}{2}S\otimes I}^2$$

$$S$$

$$+\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\otimes \frac{\lambda }{n^2}{11}^T.$$

$$E[⟨(U_{i^k}\otimes I)(\mathbf{B}{z}^k-\mathbf{B}{z}^{\star }),z^k-z^{\star }{⟩}_{H\otimes I}|{\mathcal{F}}_k]$$

$$=⟨(HE[U_{i^k}]\otimes I)(\mathbf{B}{z}^k-\mathbf{B}{z}^{\star }),z^k-z^{\star }⟩.$$

Full text

Introduction

Variance-reduced stochastic gradient (VR-SG) methods is a family of iterative optimization algorithms that combine the low per-iteration computational cost of the ordinary stochastic gradient descent and the attractive convergence properties of gradient descent. Just as ordinary stochastic gradient descent, VR-SG methods solve smooth optimization problems on finite sum form,

[TABLE]

where, for all $i\in\{1,\dots,n\}$, $f_{i}:\mathbb{R}^{N}\to\mathbb{R}$ is a convex function that is $L$-smooth, i.e., $f_{i}$ is differentiable with $L$-Lipschitz continuous gradient. These types of problems are common in model fitting, supervised learning, and empirical risk minimization which, together with the nice convergence properties of VR-SG methods, has lead to a great amount of research on VR-SG methods and the development of several different variants, e.g., [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15].

Broadly speaking, VR-SG methods form a stochastic estimate of the objective gradient by combining one or a few newly evaluated terms of the gradient with all previously evaluated terms. Classic examples of this can be seen in the SAG [1, 2] and SAGA [3] algorithms. Given some initial iterates $x^{0},y_{1}^{0},\dots,y_{n}^{0}\in\mathbb{R}^{N}$ and step-size $\lambda>0$, SAGA samples $i^{k}$ uniformly from $\{1,\dots,n\}$ and then updates the iterates as

[TABLE]

for $k\in\{0,1,\dots\}$. The update of $x^{k+1}$ is said to be unbiased since the expected value of $x^{k+1}$ at iteration $k$ is equal to an ordinary gradient descent update. This is in contrast to the biased SAG, which is identical to SAGA except that the update of $x^{k+1}$ is

[TABLE]

and the expected value of $x^{k+1}$ now includes a term containing the old gradients $\frac{1}{n}\sum_{i=1}^{n}y_{i}^{k}$. Although SAG shows that unbiasedness is not essential for the convergence of VR-SG methods, the effects of this bias are unclear. The majority of VR-SG methods are unbiased but existing works have not established any clear advantage of either the biased SAG or the unbiased SAGA. This paper will examine the effect of bias and its interplay with different problem assumptions for SAG/SAGA-like methods.

Problem and Algorithm

Instead of solving (1) directly, we consider a closely related but more general root-finding problem. Throughout the paper, we consider the Euclidean space $\mathbb{R}^{N}$ and the problem of finding $x\in\mathbb{R}^{N}$ such that

[TABLE]

where $R_{i}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $\frac{1}{L}$-cocoercive—see Section 2—for all $i\in\{1,\dots,n\}$. Since $L$-smoothness of a convex function is equivalent to $\frac{1}{L}$-cocoercivity of the gradient [16, Corollary 18.17], the smooth optimization problem in (1) can be recovered by setting $R_{i}=\nabla f_{i}$ for all $i\in\{1,\dots,n\}$ in (2). Problem (2) is also interesting in its own right with it and the closely related fixed point problem of finding $x\in\mathbb{R}^{N}$ such that $x=(\operatorname{Id}-\alpha R)x$ where $\alpha\in(0,2L^{-1})$ both having applications in for instance feasibility and non-linear signal recovery problems, see [17, 18, 19] and the references therein. To solve this problem, we present the Stochastic Variance Adjusted Gradient (SVAG) algorithm.

SVAG is heavily inspired by SAG and SAGA with both being special cases, $\theta=1$ and $\theta=n$ respectively. Just like SAG and SAGA, in each iteration, SVAG evaluates one operator $R_{i^{k}}$ and stores the results in $y_{i^{k}}^{k+1}$. An estimate of the full operator is then formed as

[TABLE]

The scalar $\theta$ determine how much weight should be put on the new information gained from evaluating $R_{i^{k}}x^{k}$. If the innovation, $R_{i^{k}}x^{k}-y_{i^{k}}^{k}$, is highly correlated with the total innovation, $Rx^{k}-\tfrac{1}{n}\sum_{j=1}^{n}y_{j}^{k}$, a large innovation weight $\theta$ can be chosen and vice versa. The innovation weight $\theta$ also determines the bias of SVAG. Taking the expected value $\widetilde{R}^{k}$ given the information at iteration $k$ gives

[TABLE]

which reveals that $\widetilde{R}^{k}$ is an unbiased estimate of $Rx^{k}$ if $\theta=n$, i.e., in the SAGA case. Any other choice, for instance SAG where $\theta=1$, yields a bias towards $\tfrac{1}{n}\sum_{j=1}^{n}y_{j}^{k}$.

Contribution

The theory behind finding roots of monotone operators in general, and cocoercive operators in particular, has been put to good use when analyzing first order optimization methods, examples include [20, 16, 21, 22, 23, 24]. For instance can both the proximal-gradient and ADMM methods be seen as instances of classic root-finding fixed-point iterations and analyzed as such, namely forward-backward and Douglas–Rachford respectively. The resulting analyses can often be simple and intuitive and even though the root-finding formulation is more general—not all cocoercive operators are gradients of convex functions—they are not necessarily more conservative. For example, analyzing proximal-gradient as forward-backward splitting yield the same rates and step-size conditions as analyzing it as a minimization method in the smooth/cocoercive setting, see for instance [25, Theorem 2.1.14] and [16, Example 5.18 and Proposition 4.39]. However, the main contribution of this paper is to show that the same is not true for VR-SG methods, in particular it is not true for SVAG when it is biased.

The results consist of two main convergence theorems for SVAG: one in the cocoercive operator case and one in the cocoercive gradient case, the later being equivalent to the minimization of a smooth and convex finite sum. Both of these theorems match or improve upon previously known results for the SAG and SAGA special cases. Comparing the two settings reveal that SVAG can use significantly larger step-sizes, with faster convergence as a result, in the cocoercive gradient case compared to the general cocoercive operator case. In the operator case, an upper bound on the step-size that scales as $\mathcal{O}(n^{-1})$ is found where $n$ is the number of terms in (2). However, the restrictions on the step-size loosen with reduced bias and the unfavorable $\mathcal{O}(n^{-1})$ scaling disappears completely when SVAG is unbiased. In the gradient case, this bad scaling never occurs, regardless of bias. We provide examples in which SVAG diverges with step-sizes larger than the theoretical upper bounds in the operator case. Since the gradient case is proven to converge with much larger step-sizes, this verifies the difference between the convergence behavior of cocoercive operators and gradients.

These results indicate that it is inadvisable to only rely on more general monotone operator theory and not explicitly use the gradient property when analyzing VR-SG methods meant for optimization. However, the large impact of bias in the cocoercive operator setting also raises the question regarding its importance in other non-gradient settings as well. One such setting of interest, where the operators are not gradients of convex functions, is the case of saddle-point problems. These problems are of importance in optimization due to their use in primal-dual methods but recently they have also gained a lot of attention due to their applications in the training of GANs in machine learning. Because of this, and due to the attractive properties of VR-SG methods in the convex optimization setting, efforts have gone into applying VR-SG methods to saddle-point problems as well [26, 27, 28, 29, 30]. Most of these efforts have been unbiased, something our analysis suggests is wise. With that said, it is important to note that our analysis is often not directly applicable due the fact that saddle-point problems rarely are cocoercive.

The main reason for the recent rise in popularity of variance-reduced stochastic methods is their use in the optimization setting, but, although bias plays a big role in the cocoercive operator case, our results are not as clear in this setting. For instance, the theoretical results for the SAG and SAGA special cases yield identical rates and step-size conditions with no clear advantage to either special case. Further experiments are therefore performed where several different choices of bias in SVAG are examined on a set of logistic regression and SVM optimization problems. However, the results of these experiments are in line with existing works with no significant advantage of any particular bias choice in SVAG, these choices include both SAG and SAGA. Although the performance difference is significant in some cases, no single choice of bias performs best for all problems and all bias choices eventually converge with the same rate in the majority of the cases. Furthermore, the theoretical maximal step-size can routinely be exceeded in these experiments, indicating that there is room for further theoretical improvements.

Related Work

There is a large array of options for solving (2). For $n\in\{1,2,3,4\}$, several operator splitting methods exist with varying assumptions on the operator properties, see for instance [31, 32, 33, 34, 24, 22] and the references therein. However, while these methods also can be applied for larger $n$ by simply regrouping the terms, they do not utilize the finite sum structure of the problem. Algorithms have therefore been designed to utilize this structure for arbitrary large $n$ with the hopes of reducing the total computational costs, e.g., [35, 36, 37, 19]. In particular the problem and method in [19] is closely related to the root-finding problem and algorithm considered in this paper.

Using the notation of [19], when $T_{0}=\operatorname{Id}$, the fixed point problem of [19] can be mapped to (2) via $R_{i}=\omega_{i}(\operatorname{Id}-T_{i})$ and vice verse. 111If $T_{i}$ is $\alpha_{i}$-averaged, as assumed in [19], $R_{i}$ is $(2\alpha_{i}\omega_{i})^{-1}$-cocoercive. Many applications considered in [19] can therefore, at least in part, be tackled with our algorithm as well. In particular, the problem of finding common fixed points of firmly nonexpansive operators can directly be solved by our algorithm. However, [19] is more general in that it allows for $T_{0}\neq\operatorname{Id}$ and works in general real Hilbert spaces. Comparing with the algorithm of [19] we see that, just as our algorithm is a generalization of SAG/SAGA, it can be seen as a generalization of Finito [8], another classic VR-SG method. It generalize Finito in several way, for instance it allows for an additional proximal/backward step and it replace the stochastic selection with a different selection criteria. However, in the optimization setting it still suffers from the same drawback as Finito when compared to SAG/SAGA-like algorithms. It still needs to store a full copy of the iterate for each term in objective. Since SAG, SAGA, and SVAG only need to store the gradient of each term, they can utilize any potential structure of the gradients to reduce the storage requirements [1]. Although the differences above are interesting in their own right, the notion of bias we examine in this paper is not applicable to Finito-like algorithms.

SAG and SAGA were compared in [3] but with no direct focus on the effects of bias. Other examples of research on SAG and SAGA include acceleration, sampling strategy selection, and ways to reduce the memory requirement [6, 38, 39, 40, 41, 42]. However, none of these works, including [41] that was written by the authors, analyze the biased case we consider in this paper. Even the works considering non-uniform sampling of gradients [38, 39, 40, 41] perform some sort of bias correction in order to remain unbiased. Furthermore, in order to keep the focus on the effects of the bias we have refrained from bringing in such generalizations into this work, making it distinct from the above research. To the authors’ knowledge, the only theoretical convergence result for biased VR-SG methods are the ones for SAG [1, 2]. But, since they only consider SAG, they fail to capture the breadth of SVAG and our proof is the first to simultaneously capture SAG, SAGA, and more.

Since the release of the first preprint of this paper, [43] has also provided a proof covering the gradient case of both SAG and SAGA, and some choices of bias in SVAG. All though [43] does not consider cocoercive operators, it is some sense more general with them considering a general biased stochastic estimator of the gradient. This generality comes at the cost of a more conservative analysis with their step-size scaling with $\mathcal{O}(n^{-1})$ in all cases.

Preliminaries and Notation

Let $\mathbb{R}$ denote the real numbers and let the natural numbers be denoted $\mathbb{N}=\{0,1,2,\dots\}$. Let ${\langle}{\cdot},{\cdot}{\rangle}$ denote the standard Euclidean inner product and ${\|}{\cdot}{\|}=\sqrt{{\langle}{\cdot},{\cdot}{\rangle}}$ the standard 2-norm. The scaled inner product and norm on we denote as ${\langle}{\cdot},{\cdot}{\rangle}_{\Sigma}={\langle}{\Sigma(\cdot)},{\cdot}{\rangle}$ and ${\|}{\cdot}{\|}_{\Sigma}=\sqrt{{\langle}{\cdot},{\cdot}{\rangle}_{\Sigma}}$ where $\Sigma$ is a positive definite matrix. If $\Sigma$ is not positive definite, ${\|}{\cdot}{\|}_{\Sigma}$ is not a norm but we keep the notation for convenience.

Let $n$ be the number of operators in (2). The vector $\mathbf{1}$ is the vector of all ones in $\mathbb{R}^{n}$ and $e_{i}$ be the vector in $\mathbb{R}^{n}$ of all zeros except the $i$:th element that contains a $1$. The matrix $I$ is an identity matrix with the size derived from context and $E_{i}=e_{i}e_{i}^{T}$.

The symbol $\otimes$ denotes the Kronecker product of two matrices. The Kronecker product is linear in both arguments and the following properties hold

[TABLE]

In the last property it is assumed that the dimensions are such that the matrix multiplications are well defined. The eigenvalues of $A\otimes B$ are given by

[TABLE]

where $\tau_{i}$ and $\mu_{j}$ are the eigenvalues of $A$ and $B$ respectively.

The Cartesian product of two sets $C_{1}$ and $C_{2}$ is defined as

[TABLE]

From this definition we see that if $C_{1}$ and $C_{2}$ are closed and convex, so is $C_{1}\times C_{2}$.

Let $X^{\star}$ be the set of all solutions of $\eqref{eq:zerprob}$,

[TABLE]

and define $Z^{\star}$ as the set of primal-dual solutions

[TABLE]

Assuming they exists, $x^{\star}$ denotes a solution to (2) and $z^{\star}$ denotes a primal-dual solution, i.e., $x^{\star}\in X^{\star}$ and $z^{\star}\in Z^{\star}$.

A single valued operator $R:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $\frac{1}{L}$-cocoercive if

[TABLE]

holds for all $x,y\in\mathbb{R}^{N}$. An operator that is $\frac{1}{L}$-cocoercive is $L$-Lipschitz continuous. The set of zeros of a cocoercive operator $R$ is closed and convex.

A differentiable convex function $f:\mathbb{R}^{N}\to\mathbb{R}$ is called $L$-smooth if the gradient is $\frac{1}{L}$-cocoercive. Equivalently, a differentiable convex function is $L$-smooth if

[TABLE]

holds for all $x,y\in\mathbb{R}^{N}$.

If $f_{i}:\mathbb{R}^{N}\to\mathbb{R}$ is a differentiable convex functions for each $i\in\{1,\dots,n\}$, the minimization of $\sum_{i=1}^{n}f_{i}(x)$ is equivalent to (2) with $R_{i}=\nabla f_{i}$.

For more details regarding monotone operators and convex functions see [25, 16].

To establish almost sure sequence convergence of the stochastic algorithm, the following propositions will be used. The first is from [44] and establishes convergence of non-negative almost super-martingales. The second is based on [45] and provides the tool to show almost sure sequence convergence.

Proposition 2.1**.**

Let $(\Omega,\mathcal{F},P)$ be a probability space and $\mathcal{F}_{0}\subset\mathcal{F}_{1}\subset\dots$ be a sequence of sub-$\sigma$-algebras of $\mathcal{F}$. For all $k\in\mathbb{N}$, let $z^{k}$, $\beta^{k}$, $\xi^{k}$ and $\zeta^{k}$ be non-negative $\mathcal{F}_{k}$-measurable random variables. If $\sum_{i=0}^{\infty}\beta^{i}<\infty$, $\sum_{i=0}^{\infty}\xi^{i}<\infty$ and

[TABLE]

hold almost surely for all $k\in\mathbb{N}$, then $z^{k}$ converges a.s. to a finite valued random variable and $\sum_{i=0}^{\infty}\zeta^{i}<\infty$ almost surely.

Proof.

See [44, Theorem 1]. ∎

Proposition 2.2**.**

Let $Z$ be a non-empty closed subset of a finite dimensional Hilbert space $H$, let $\phi:[0,\infty)\to[0,\infty)$ be a strictly increasing function such that $\phi(t)\to\infty$ as $t\to\infty$, and let $(x^{k})_{k\in\mathbb{N}}$ be a sequence of $H$-valued random variables. If $\phi({\|}{x^{k}-z}{\|})$ converges a.s. to a finite valued non-negative random variable for all $z\in Z$, then the following hold:

1. (i)

$(x^{k})_{k\in\mathbb{N}}$* is bounded almost surely.* 2. (ii)

Suppose the cluster points of $(x^{k})_{k\in\mathbb{N}}$ are a.s. in $Z$, then $(x^{k})_{k\in\mathbb{N}}$ converge a.s. to a $Z$-valued random variable.

Proof.

In finite dimensional Hilbert spaces, these two statements are the same as statements (ii) and (iv) of [45, Proposition 2.3]. Hence, consider the proof of [45, Proposition 2.3] restricted to finite dimensional Hilbert spaces. The proof of (ii) in [45, Proposition 2.3] only relies on the a.s. convergence of $\phi({\|}{x^{k}-z}{\|})$ and hence is implied by the assumptions of this proposition. This proves our first statement. The proof of (iv) in [45, Proposition 2.3] only relies on (iii) of [45, Proposition 2.3] which in turn is implied by (ii) of [45, Proposition 2.3], i.e., our first statement. This proves our second statement. ∎

Convergence

Throughout the analysis we will use the following two assumptions on the operators of (2).

Assumption 3.1**.**

For each $i\in\{1,\dots,n\}$, let $R_{i}$ be $\frac{1}{L}$-cocoercive and $X^{\star}\neq\emptyset$, i.e., (2) has at least one solution.

Assumption 3.2**.**

For each $i\in\{1,\dots,n\}$, let $R_{i}=\nabla f_{i}$ for some differentiable function $f_{i}$ and define $F=\tfrac{1}{n}\sum_{i=1}^{n}f_{i}$. Furthermore, let Assumption 3.1 hold, i.e., $f_{i}$ is $L$-smooth and convex and $\operatorname*{arg\,min}F(x)$ exists.

Reformulation

We begin by formalizing and reformulating Algorithm 1 into a more convenient form. Let $(\Omega,\mathcal{F},P)$ be the underlying probability space of Algorithm 1. The index selected at iteration $k$ is then a uniformly distributed random variable $i^{k}:\Omega\to\{1,\dots,n\}$. For each $k\in\mathbb{N}$, define the random variable $z^{k}:\Omega\to\mathbb{R}^{N(n+1)}$ as $z^{k}=(x^{k},y_{1}^{k},\dots,y_{n}^{k})$ where $x^{k}$ and $y_{i}^{k}$ for $i\in\{1,\dots,n\}$ are the iterates of Algorithm 1. Let $\mathcal{F}_{0}\subset\mathcal{F}_{1}\subset\dots$ be a sequence of sub-$\sigma$-algebras of $\mathcal{F}$ such that $z^{k}$ are $\mathcal{F}_{k}$-measurable and $i^{k}$ is independent of $\mathcal{F}_{k}$. With the operator $\mathbf{B}:\mathbb{R}^{N(n+1)}\to\mathbb{R}^{2Nn}$ defined as $\mathbf{B}(x,y_{1},\dots,y_{n})=(R_{1}x,\dots,R_{n}x,y_{1},\dots,y_{n})$, one iteration of Algorithm 1 can be written as

[TABLE]

where $z^{0}\in\mathbb{R}^{N(n+1)}$ is given and

[TABLE]

for all $i\in\{1,\dots,n\}$. The vector $e_{i}$ and the matrix $E_{i}$ are defined in Section 2.

The following lemma characterizes the zeros of $(U_{i}\otimes I)\mathbf{B}$ and hence the fixed points of (6) and Algorithm 1.

Lemma 3.1**.**

Let Assumption 3.1 hold, each $z^{\star}$ in $Z^{\star}$ is then a zero of $(U_{i}\otimes I)\mathbf{B}$ for all $i\in\{1,\dots,n\}$, i.e.

[TABLE]

Furthermore, the set $Z^{\star}$ is closed and convex and $R_{i}x^{\star}=R_{i}\bar{x}^{\star}$ for all $x^{\star},\bar{x}^{\star}\in X^{\star}$ and for all $i\in\{1,\dots,n\}$.

Proof of Lemma 3.1.

The zero statement, $0=(U_{i}\otimes I)\mathbf{B}z^{\star}$, follows from definition of $z^{\star}$. For closedness and convexity of $Z^{\star}$, we first prove that $R_{i}x^{\star}$ is unique for each $i\in\{1,\dots,n\}$. Taking $x,y\in X^{\star}$, which implies ${\textstyle\sum_{i=1}^{n}}R_{i}x={\textstyle\sum_{i=1}^{n}}R_{i}y=0$, and using cocoercivity (4) of each $R_{i}$ gives

[TABLE]

hence must $R_{i}x=R_{i}y$ for all $i\in\{1,\dots,n\}$. The set $Z^{\star}$ is a Cartesian product of $X^{\star}$ and the points $r_{i}=R_{i}x^{\star}$ for $i\in\{1,\dots,n\}$ and any $x^{\star}\in X^{\star}$. A set consisting of only one point is closed and convex and $X^{\star}$ is closed and convex since $\frac{1}{n}\sum_{i=1}^{n}R_{i}$ is cocoercive [16, Proposition 23.39], hence is $Z^{\star}$ closed and convex. ∎

The operator $\mathbf{B}$ in the reformulated algorithm can be used to enforce the following property on the sequence $(z^{k})_{k\in\mathbb{N}}$.

Lemma 3.2**.**

Let $(\Omega,\mathcal{F},P)$ be a probability space and $(z^{k})_{k\in\mathbb{N}}$ be a sequence of random variables $z^{k}:\Omega\to\mathbb{R}^{N(n+1)}$. If $\mathbf{B}z^{k}\to\mathbf{B}z^{\star}$ a.s. where $z^{\star}\in Z^{\star}$, then any cluster point of $(z^{k})_{k\in\mathbb{N}}$ will almost surely be in $Z^{\star}$.

Proof of Lemma 3.2.

Let $z$ be a cluster point of $(z^{k})_{k\in\mathbb{N}}$. Take an $\omega\in\Omega$ such that $\mathbf{B}z^{k}(\omega)\to\mathbf{B}z^{\star}$. For this $\omega$ and for all $k\in\mathbb{N}$, we define the realizations of $z$ and $z^{k}$ as

[TABLE]

where $\bar{x},\bar{y}_{1},\dots,\bar{y}_{n}\in\mathbb{R}^{N}$ and $\bar{x}^{k},\bar{y}_{1}^{k},\dots,\bar{y}_{n}^{k}\in\mathbb{R}^{N}$ for all $k\in\mathbb{N}$.

Since $\mathbf{B}\bar{z}^{k}\to\mathbf{B}z^{\star}$ we directly have $\bar{y}_{i}^{k}\to R_{i}x^{\star}$ for $x^{\star}\in X^{\star}$ and hence must $\bar{y}_{i}=R_{i}x^{\star}$ for all $i\in\{1,\dots,n\}$. Note, $R_{i}x^{\star}$ is independent of which $x^{\star}\in X^{\star}$ was chosen, see Lemma 3.1. Furthermore, $\mathbf{B}\bar{z}^{k}\to\mathbf{B}\bar{z}^{\star}$ implies that $R_{i}\bar{x}^{k}\to R_{i}x^{\star}$ for all $i\in\{1,\dots,n\}$. Let $(\bar{x}^{k(l)})_{l\in\mathbb{N}}$ be a subsequence converging to $\bar{x}$, then

[TABLE]

as $l\to\infty$ where $L$-Lipschitz continuity of $\frac{1}{n}\sum_{i=1}^{n}R_{i}$ was used. This concludes that $\bar{x}\in X^{\star}$ and since $\bar{y}_{i}=R_{i}x^{\star}=R_{i}\bar{x}$ for all $i\in\{1,\dots,n\}$ by Lemma 3.1, we have that $z(\omega)\in Z^{\star}$. Since this hold for any $\omega$ such that $\mathbf{B}z^{k}(\omega)\to\mathbf{B}z^{\star}$ and the set in $\mathcal{F}$ of all such $\omega$ have probability one due to the almost sure convergence of $\mathbf{B}z^{k}\to\mathbf{B}z^{\star}$, we have $z\in Z^{\star}$ almost surely. ∎

The reformulation (6) further allows us to concisely formulate two Lyapunov inequalities.

Lemma 3.3**.**

Let Assumption 3.1 hold, the update (6) then satisfies

[TABLE]

for all $k\in\mathbb{N}$ and $\xi\in[0,\frac{2\lambda}{nL}]$, where matrices $H$ and $M$ are given by

[TABLE]

and

[TABLE]

Lemma 3.4**.**

Let Assumption 3.2 hold, the update (6) then satisfies

[TABLE]

for all $k\in\mathbb{N}$, where $K=\begin{bmatrix}1&\frac{\lambda}{n}\mathbf{1}^{T}\end{bmatrix}$ and

[TABLE]

Proof of Lemma 3.3.

Take $k\in\mathbb{N}$, note that since $U_{i^{k}}$ is independent of $\mathcal{F}_{k}$ and $z_{k}$ is $\mathcal{F}_{k}$-measurable we have

[TABLE]

The matrix $H\operatorname{\mathbb{E}}[U_{i^{k}}]$ is given by

[TABLE]

see the supplementary material for verification of this and other matrix identities. We also note that

[TABLE]

Taking $\xi\in[0,\frac{2\lambda}{nL}]$ and putting these two expression together yield

[TABLE]

Using $\frac{1}{L}$-cocoercivity of $R_{i}$ for each $i\in\{1,\dots,n\}$ gives

[TABLE]

Setting

[TABLE]

gives

[TABLE]

where $M=\frac{1}{2}(\bar{M}+\bar{M}^{T})$ is the matrix in the Lemma. Finally, using this inequality and $0=(U_{i^{k}}\otimes I)\mathbf{B}z^{\star}$ from Lemma 3.1 gives

[TABLE]

∎

Proof of Lemma 3.4.

Take $k\in\mathbb{N}$ and note that

[TABLE]

where $Q_{i^{k}}=\frac{\lambda}{n}\begin{bmatrix}(\theta-1)e_{i^{k}}^{T}&-(\theta-1)e_{i^{k}}^{T}\end{bmatrix}$. Furthermore, with $G=\tfrac{1}{n}\begin{bmatrix}\mathbf{1}^{T}&0\end{bmatrix}$, we have $\nabla F(x^{k})=(G\otimes I)\mathbf{B}z^{k}$. From the definition of $z^{\star}$ we have $0=(G\otimes I)\mathbf{B}z^{\star}=(Q_{i^{k}}\otimes I)\mathbf{B}z^{\star}$. Using $L$-smoothness, (5), of $F$ yields

[TABLE]

where $S_{L}=\operatorname{\mathbb{E}}[Q_{i^{k}}^{T}G+G^{T}Q_{i^{k}}-LQ_{i^{k}}^{T}Q_{i^{k}}]$.

With $D=\begin{bmatrix}0&\mathbf{1}^{T}\end{bmatrix}$ we have $(K\otimes I)z^{k}=x^{k}+\frac{\lambda}{n}(D\otimes I)\mathbf{B}z^{k}$. Using the first order convexity condition on $F$ and $0=(D\otimes I)\mathbf{B}z^{\star}=(G\otimes I)\mathbf{B}z^{\star}$ yields

[TABLE]

where $S_{C}=\frac{\lambda}{n}(D^{T}G+G^{T}D)$. Combining these two inequalities gives

[TABLE]

where $S=S_{L}+S_{C}$. ∎

Convergence Theorems

We are now ready to state the main convergence theorems for SVAG. They are stated with the notation from Algorithm 1 but are proved at the end of this section with the help of the reformulation in (6) and the lemmas above.

Theorem 3.1**.**

For all $i\in\{1,\dots,n\}$, let $(x^{k})_{k\in\mathbb{N}}$ and $(y_{i}^{k})_{k\in\mathbb{N}}$ be the sequences generated by Algorithm 1. If Assumption 3.1 hold and the step-size, $\lambda>0$, and innovation weight, $\theta\in\mathbb{R}$, satisfy

[TABLE]

then $x^{k}\rightarrow x^{\star}$ and $y_{i}^{k}\rightarrow R_{i}x^{\star}$ almost surely for all $i\in\{1,\dots,n\}$, where $x^{\star}$ is a solution to (2). For all $i\in\{1,\dots,n\}$, the residuals converge a.s. as

[TABLE]

where $c=2+|n-\theta|$ and

[TABLE]

for any $x^{\star}\in X^{\star}$.

Theorem 3.2**.**

For all $i\in\{1,\dots,n\}$, let $(x^{k})_{k\in\mathbb{N}}$ and $(y_{i}^{k})_{k\in\mathbb{N}}$ be the sequences generated by Algorithm 1. If Assumption 3.2 hold and the step-size, $\lambda>0$, and innovation weight, $\theta\in[0,n]$, satisfy

[TABLE]

then $x^{k}\rightarrow x^{\star}$ and $y_{i}^{k}\rightarrow\nabla f_{i}(x^{\star})$ almost surely, where $x^{\star}$ is a solution to (2). For all $i\in\{1,\dots,n\}$, the residuals converge a.s. as

[TABLE]

where

[TABLE]

for any $x^{\star}\in X^{\star}$.

Both Theorem 3.1 and 3.2 give the step-size condition $\lambda\in(0,\tfrac{1}{2L})$ for the SAGA special case, i.e., $\theta=n$. This is the same as the largest upper bound found in the literature [3] and appears to be tight [41]. Theorem 3.2 also give this step-size condition when $\theta=1$, i.e., SAG in the optimization case. This bound improves on upper bound of $\frac{1}{16L}\leq\lambda$ presented in [2].

In the cocoercive operator setting with $\theta\neq n$, Theorem 3.1 gives a step-size condition that scales with $n^{-1}$. This step-size scaling is significantly worse compared to the gradient case in Theorem 3.2 in which the step-size’s dependence on $n$ is $\mathcal{O}(1)$ for all $\theta$. This difference is indeed real and not an artifact of the analysis since we in Section 4 present a problem for which the cocoercivity result appears to be tight. A consequence of this unfavorable step-size scaling in the operator setting is slow convergence. There is therefore little reason to use anything else than $\theta=n$ in SVAG when $R_{i}$ is not a gradient of a smooth function for all $i\in\{1,\dots,n\}$.

The rates of Theorem 3.1 and 3.2 are of $\mathcal{O}(\frac{1}{t+1})$ type with two sets of multiplicative factors. One factor which only depend on the algorithm parameters, $\frac{n}{\lambda(L^{-1}-\lambda c)}$, and one set which depend on how the algorithm initialization relates to the solution set, $C_{R}$ and $C_{R}+C_{F}$. The initialization dependent factors also depend on the algorithm parameters, but, since knowing the exact dependency requires knowing the solution set, we will not attempt to tune the parameters to decrease this factor. Only considering the first factor, the rate becomes better if $c$ is decreased and, since $c$ is independent of $\lambda$, the best choice of step-size is $\lambda=(2Lc)^{-1}$. This means that $\lambda=(4L)^{-1}$ and $\theta=n$ are the best parameter choices in the cocoercive operator setting. In the optimization case the best step-size is also $\lambda=(4L)^{-1}$ but the innovation weight can be selected as either $\theta=n$ or $\theta=1$.

However, in the optimization case we do not believe that these theoretical rates reflects real world performance and parameter choices based on them might therefore not perform particularly well. We base this belief on our experience with numerical experiments. For $\theta\neq n$ and $\theta\neq 1$, we have not found any optimization problem where the step-size condition in Theorem 3.2 appears to be tight. Also, using $\lambda=(2Lc)^{-1}$ as suggested by the Theorem 3.2 can in some cases lead to impractically small step-sizes. For instance, if $\lambda=(2Lc)^{-1}$ was used in the experiments in Section 4, a couple of the experiments would have step-sizes over 1000 times smaller than the ones used now. One can of course not disprove a worst case analysis with experiments but we still feel they indicate a conservative analysis, even though the analysis improves on the previous best results.

Proof of Theorem 3.1.

Apply Lemma 3.3 with $\xi=0$, the iterates given by (6) then satisfy the following for all $z^{\star}\in Z^{\star}$,

[TABLE]

Assuming $H\succ 0$ and $2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]\succ 0$ can Proposition 2.1 be applied. We will later prove that this assumption indeed does hold. Proposition 2.1 gives a.s. summability of ${\|}{\mathbf{B}z^{k}-\mathbf{B}z^{\star}}{\|}_{(2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}])\otimes I}^{2}$ and hence will $\mathbf{B}z^{k}\to\mathbf{B}z^{\star}$ almost surely. Lemma 3.2 then gives that all cluster points of $(z^{k})_{k\in\mathbb{N}}$ are in $Z^{\star}$ almost surely. Finally, since Proposition 2.1 ensures the a.s. convergence of ${\|}{z^{k}-z^{\star}}{\|}_{H\otimes I}^{2}$ and since $\mathbb{R}^{N(n+1)}$ with the inner product ${\langle}{(H\otimes I)\cdot},{\cdot}{\rangle}$ is a finite dimensional Hilbert space, Proposition 2.2 gives the almost sure convergence of $z^{k}\to z^{\star}\in Z^{\star}$.

There always exists a $\lambda$ such that $2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]$ and $H$ are positive definite. First we show that $H\succ 0$ always holds for $\lambda>0$. Taking the Schur complement of $1$ in $H$ gives

[TABLE]

Hence is $H\succ 0$ since the Schur complement is positive definite.

We now show $2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]\succ 0$. Straightforward algebra, see the supplementary material, yields

[TABLE]

Positive definiteness of this matrix is established by ensuring positivity of the smallest eigenvalue $\sigma_{\min}$. The smallest eigenvalue $\sigma_{\min}$ is greater than the sum of the smallest eigenvalue of each term. For the eigenvalues of the Kronecker products, see (3). This gives that

[TABLE]

Since $\lambda>0$ by assumption, if

[TABLE]

we have that $\sigma_{\min}>0$ and that $2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]$ is positiv definite.

Rates are gotten by taking the total expectation of (8) and adding together the inequalities from $k=0$ to $k=t$, yielding

[TABLE]

Putting in the lower bound on $\sigma_{\min}$ and rearranging yields

[TABLE]

From the definition of $H$ in Lemma 3.3 we have

[TABLE]

where $z^{\star}=(x^{\star},R_{1}x^{\star},\dots,R_{n}x^{\star})$. Since this hold for any $z^{\star}\in Z^{\star}$ and hence any $x^{\star}\in X^{\star}$, the results of theorems follows by minimizing the RHS over $x^{\star}\in X^{\star}$. Note, since $R_{i}x^{\star}$ constant for all $x^{\star}\in X^{\star}$, the objective is convex and, since $X^{\star}$ is closed and convex, the minimum is then attained. ∎

Proof of Theorem 3.2.

Combining Lemma 3.3 and 3.4 yield

[TABLE]

which holds for all $k\in\mathbb{N}$, $\xi\in[0,\frac{2\lambda}{nL}]$, and $z^{\star}\in Z^{\star}$. Since $H\succ 0$ for $\lambda>0$, see the proof of Theorem 3.1, the first term is non-negative while the second term is non-negative if $\theta\leq n$. From cocoercivity of $\nabla F$, the last term is non-positive and we assume, for now, that there exists $\lambda>0$ and $\frac{2\lambda}{nL}\geq\xi>0$ such that $2M-\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]+\lambda(n-\theta)S-\xi I\succ 0$, making the third term non-positive.

Applying Proposition 2.1 gives the a.s. summability of

[TABLE]

Since both term are positive, both terms are a.s. summable. From the first term we have the a.s. convergence of $\mathbf{B}z^{k}\to\mathbf{B}z^{\star}$ and Lemma 3.2 then gives that all cluster points of $(z^{k})_{k\in\mathbb{N}}$ are almost surely in $Z^{\star}$. For the second term we note that by convexity we have

[TABLE]

and $F(x^{k})-F(x^{\star})$ then is summable a.s. since $\xi nL>0$. Using smoothness of $F$, (5) and the notation from (7) gives

[TABLE]

since $(G\otimes I)\mathbf{B}z^{k}\to(G\otimes I)\mathbf{B}z^{\star}=0$ and $(D\otimes I)\mathbf{B}z^{k}\to(D\otimes I)\mathbf{B}z^{\star}=0$ almost surely. Therefore we have the a.s. convergence of $F((K\otimes I)z^{k})-F(x^{\star})\to 0$.

From Proposition 2.1 we can also conclude that ${\|}{z^{k}-z^{\star}}{\|}_{H\otimes I}^{2}+2\lambda(n-\theta)(F((K\otimes I)z^{k})-F(x^{\star}))$ a.s. converge to a non-negative random variable. Since $F((K\otimes I)z^{k})-F(x^{\star})\to 0$ a.s. we have that ${\|}{z^{k}-z^{\star}}{\|}_{H\otimes I}^{2}$ also must a.s. converge to a non-negative random variable. Proposition 2.2 then give the almost sure convergence of $(z^{k})_{k\in\mathbb{N}}$ to $Z^{\star}$.

We now show that there exists $\lambda>0$ and $\xi>0$ such that

[TABLE]

We show positive definiteness by ensuring that the smallest eigenvalue is positive. The smallest eigenvalue $\sigma_{\min}$ is greater than the sum of the smallest eigenvalues of each term,

[TABLE]

Assuming $\lambda\leq\frac{1}{2L}$ yields the following lower bound on the smallest eigenvalue

[TABLE]

Selecting

[TABLE]

which satisfy the assumption $\frac{2\lambda}{nL}\geq\xi>0$, yield $\sigma_{\min}\geq\xi$. Since $\lambda>0$ by assumption, if

[TABLE]

we have that $\sigma_{\min}\geq\xi>0$ and hence that the examined matrix is positive definite. Furthermore, if $\lambda$ satisfies the above inequality it also satisfies the assumption $\lambda\leq\frac{1}{2L}$.

Rates are gotten in the same way as for Theorem 3.1, the total expectation is taken of the Lyapunov inequality at the beginning of the proof and the inequalities are summed from $k=0$ to $k=t$.

[TABLE]

Inserting the lower bound on $\sigma_{\min}$, rearranging and minimizing over $x^{\star}\in X^{\star}$ yield the results of the theorem. ∎

Numerical Experiments

A number of experiments, outlined below, were performed to verify the tightness of the theory in the cocoercive operator case and examine the effect of bias in the cocoercive gradient case. The experiments were implemented in Julia [46] and, together with several other VR-SG methods, can be found at https://github.com/mvmorin/VarianceReducedSG.jl.

Cocoercive Operators Case

In order for the difference between cocoercive operators and cocoercive gradients to not be an artifact of our analysis, the results in the operator case can not be overly conservative. We therefore construct a cocoercive operator problem for which the results appear to be tight, thereby verifying the difference. Consider problem (2) where the operator $R_{i}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ is an averaged rotation

[TABLE]

for all $i\in\{1,\dots,n\}$ and some $\tau\in[0,2\pi)$. The operators are $1$-cocoercive and the zero vector is the only solution to (2) if $\tau\neq\pi$. The step-size condition from Theorem 3.1 appears to be tight for $\theta\in[0,n]$ when the angle of rotation $\tau$ approaches $\pi$. We therefore let $\tau=\frac{179}{180}\pi$ and solve the problem with different configurations of step-size $\lambda$ and innovation weight $\theta$.

Figure 1 displays the relative distance to the solution after $100n$ iterations of SVAG together with the upper bound on the step-size. When $\theta\in[0,n]$ and $\lambda$ exceeds the upper bound, the distance to the solution increases for both $n=100$ and $n=10000$, i.e., the method does not converge. Hence, for $\theta\in[0,n]$, the step-size bound in Theorem 3.1 appears to be tight. However, it is noteworthy that for this particular problem it seems beneficial to exceed the step-size bound when $\theta>n$.

Cocoercive Gradients Case

Since, as we stated in Section 3.2, we do not believe that the theoretical rates are particularly tight in the optimization case, we examine the effects of the bias numerically. These experiments can of course not be exhaustive and we choose to focus on only the bias parameter $\theta$ and perform all experiments with the same step-size. This also demonstrate why we believe the analysis to be conservative since the chosen step-size in some cases are a 1000 times larger than upper bound from Theorem 3.2. Convergence with this large of a step-size have also been seen elsewhere with both [2] and [43] disregarding their own the theoretical step-size conditions.

The experiments are done by performing a rough parameter sweep over the innovation weight $\theta$ on two different binary classification problems and we will look for patterns in how the convergence is affected. The first problem is logistic regression,

[TABLE]

The second is SVM with a square hinge loss,

[TABLE]

where $\gamma>0$ is a regularization parameter. In both problems are $y_{i}\in\{-1,1\}$ the label and $a_{i}\in\mathbb{R}^{N}$ the features of the $i$th training data point. Note, although not initially obvious, $\max(0,\cdot)^{2}$ is convex and differentiable with Lipschitz continuous derivative and the second problem is therefore indeed smooth. The logistic regression problem does not necessarily have a unique solution and the distance to the solution set is therefore hard to estimate. For this reason, we examine the convergence of $\|\nabla F(x^{k})\|\to 0$ instead of the distance to the solution set.

The datasets for both these classification problems are taken from the LibSVM[47] collection of datasets. The number of examples in the datasets varies between $n=683$ and $n=60,000$ while the number of features is between $N=10$ and $N=5,000$. Two of the datasets, mnist.scale and protein, consist of more than 2 classes. These are converted to binary classification problems by grouping the different classes into two groups. For the digit classification dataset mnist.scale, the digits are divided into the groups 0-4 and 5-9. For the protein data set, the classes are grouped as 0 and 1-2. The results of solving the classification problems above can be found in Figures 2 and 3.

From Figures 2 and 3 it appears like the biggest difference between the innovation weights are in the early stages of the convergence. Most innovation weight choices appear to eventually converge with the same rate. In the cases where this does not happen, the fastest converging choice of innovation weight actually reaches machine precision. It is therefore not possible to say whether these cases would eventually reach the same rate as well. Since none of the choices of $\theta$ appears to consistently be at a significant disadvantage, even though the step-size used exceeds the upper bound in Theorem 3.2 when $\theta=0.1n$ and $\theta=0.01n$, we conjecture that the asymptotic rates for a given step-size is independent of $\theta$.

The initial phase can clearly have a large impact on the convergence and it can therefore still be a benefit to tuning the bias. However, comparing the different choices of innovation weight yields no clear conclusion since no single choice of innovation weight consistently outperforms another. In most cases do the lower bias choices—$\theta=n$ (SAGA) or $\theta=0.1n$—seem perform best but, when they do not, the high bias choices—$\theta=1$ (SAG) and $\theta=0.01n$—perform significantly better. Another observation is that lowering $\theta$ increases any oscillations if they are present. We speculate that it is due to the increased inertia and we also believe that this inertia is what allows the lower innovation weights to sometimes perform better.

Conclusion

We presented SVAG, a variance-reduced stochastic gradient method with adjustable bias and with SAG and SAGA as special cases. It was analyzed in two scenarios, one being the minimization problem of a finite sum of functions with cocoercive gradients and the other being finding a root of a finite sum of cocoercive operators. The analysis improves on the previously best known analyses in both settings and, more significantly, the two different scenarios gave significantly different convergence conditions for the step-size. In the cocoercive operator setting a much more restrictive condition was found and it was verified numerically. This difference is not present in ordinary gradient descent and can therefore easily be overlooked, however, these results suggest that is inadvisable in the variance-reduced stochastic gradient setting.

The theoretical results in the minimization case was further examined with numerical experiments. Several choices of bias was examined but we did not find the same dependence on the bias that the theory suggests. In fact, the asymptotic convergence behavior was similar for the different choices of bias, indicating that further improvements of the theory is still needed. The bias mainly impacted the early stages of the convergence and in a couple of cases this impact was significant. There might therefore still be benefits to tuning the bias to the particular problem but further work is needed to efficiently do so.

Funding

This work is funded by the Swedish Research Council via grant number 2016-04646.

Appendix A Matrix Identities

Here we verify the matrix identities used in the proofs.

[TABLE]

[TABLE]

For $\operatorname{\mathbb{E}}[U_{i^{k}}^{T}HU_{i^{k}}]$, first note that $U_{i^{k}}$ can be written as

[TABLE]

which gives

[TABLE]

The first term is

[TABLE]

The second term is

[TABLE]

The third and forth term are

[TABLE]

This results in

[TABLE]

and

[TABLE]

[TABLE]

We have $S=S_{L}+S_{C}$ where

[TABLE]

and

[TABLE]

[TABLE]