Convergence Analysis of Inexact Randomized Iterative Methods

TL;DR

This paper analyzes the convergence rates of inexact randomized iterative methods, extending classical algorithms by allowing inexact sub-problem solutions and providing complexity results and numerical validation.

Contribution

It introduces and analyzes inexact variants of several randomized methods, broadening their applicability and understanding of convergence behavior.

Findings

Inexact methods converge under various assumptions on inexactness.

Numerical experiments show inexactness can improve computational efficiency.

Many popular randomized algorithms are special cases of the proposed framework.

Abstract

In this paper we present a convergence rate analysis of inexact variants of several randomized iterative methods. Among the methods studied are: stochastic gradient descent, stochastic Newton, stochastic proximal point and stochastic subspace ascent. A common feature of these methods is that in their update rule a certain sub-problem needs to be solved exactly. We relax this requirement by allowing for the sub-problem to be solved inexactly. In particular, we propose and analyze inexact randomized iterative methods for solving three closely related problems: a convex stochastic quadratic optimization problem, a best approximation problem and its dual, a concave quadratic maximization problem. We provide iteration complexity results under several assumptions on the inexactness error. Inexact variants of many popular and some more exotic methods, including randomized block Kaczmarz, randomized Gaussian Kaczmarz and randomized block coordinate descent, can be cast as special cases. Numerical experiments demonstrate the benefits of allowing inexactness.

Figures

Tables

Table 1. Table 1 : Summary of the iteration complexity results obtained in this paper. ω 𝜔 \omega denotes the stepsize (relaxation parameter) of the method. In all cases, x ∗ = Π ℒ , 𝐁 ​ ( x 0 ) subscript 𝑥 subscript Π ℒ 𝐁 subscript 𝑥 0 x_{*}=\Pi_{{\cal L},\mathbf{B}}(x_{0}) and ρ = 1 − ω ​ ( 2 − ω ) ​ λ min + ∈ ( 0 , 1 ) 𝜌 1 𝜔 2 𝜔 superscript subscript 𝜆 0 1 \rho=1-\omega(2-\omega)\lambda_{\min}^{+}\in(0,1) are the quantities appear in the convergence results (here λ min + superscript subscript 𝜆 \lambda_{\min}^{+} denotes the minimum non zero eigenvalue of matrix 𝐖 𝐖 \mathbf{W} , see equation ( 19 )). Inexactness parameter q 𝑞 q is chosen always in such a way to obtain linear convergence and it can be seen as the quantity that controls the inexactness. In all theorems the quantity of convergence is 𝔼 ​ [ ‖ x k − x ∗ ‖ 𝐁 2 ] 𝔼 delimited-[] subscript superscript norm subscript 𝑥 𝑘 subscript 𝑥 2 𝐁 \mathbb{E}[\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}] (except in Theorem  1 where we analyze 𝔼 [ ∥ x k − x ∗ ∥ 𝐁 ] ) \mathbb{E}[\|x_{k}-x_{*}\|_{\mathbf{B}}]) . As we show in Section  5 , under similar assumptions, iSDSA has exactly the same convergence with iBasic but the upper bounds of the third column are related to the dual function values 𝔼 ​ [ D ​ ( y ∗ ) − D ​ ( y 0 ) ] 𝔼 delimited-[] 𝐷 subscript 𝑦 𝐷 subscript 𝑦 0 {\mathbb{E}\left[D(y_{*})-D(y_{0})\right]} .

Assumption on the Inexactness error ${ϵ}_k$	$\omega$	Upper Bounds	Theorem
Assumption Assumption 1$a$	$(0,2)$	${\rho }^{k/2}{\Vert x_0-x_{\ast }\Vert }_{𝐁}+{\sum }_{i=0}^{k-1}{\rho }^{\frac{k-1-i}{2}}{\sigma }_i$	1
Assumption Assumption 1$b$	$(0,2)$	${\left(\sqrt{\rho }+q\right)}^{2k}{\Vert x_0-x_{\ast }\Vert }_{𝐁}^2$	2
Assumptions Assumption 1,Assumption 2	$(0,2)$	${\rho }^k{\Vert x_0-x_{\ast }\Vert }_{𝐁}^2+{\sum }_{i=0}^{k-1}{\rho }^{k-1-i}{\overline{\sigma }}_i^2$	3(i)
Assumptions Assumption 1$b$,Assumption 2	$(0,2)$	${\left(\rho +q^2\right)}^k{\Vert x_0-x_{\ast }\Vert }_{𝐁}^2$	3(ii)
Assumptions Assumption 1$c$,Assumption 2	$(0,2)$	${\left(\rho +q^2{\lambda }_{\min }^{+}\right)}^k{\Vert x_0-x_{\ast }\Vert }_{𝐁}^2$	3(iii)

Table 2. Table 2 : The exact algorithms under study with the potentially expensive to compute key sub-problems of their update rule. The inexact update rules are presented in the last column for both notions of inexactness (abstract and more structured). We use ⋅ ∗ superscript ⋅ {\cdot}^{*} to define the important quantity that needs to be computed exactly in the update rule of each method and ⋅ ≈ superscript ⋅ {\cdot}^{\approx} to indicate the proposed inexact variant.

Exact Algorithms	Key Subproblem (problem that we solve inexactly)	Inexact Update Rules (abstract and structured inexactness error)
SGD (4)	Exact computation of ${\lambda }_k^{\ast }$, where ${\lambda }_k^{\ast }=\arg {\min }_{\lambda :{𝐌}_k\lambda =d_k}\Vert \lambda \Vert$. Appears in the computation of $\nabla f_{{𝐒}_k}\left(x_k\right)=-{𝐁}^{-1}{𝐀}^{\top }{𝐒}_k{\lambda }_k^{\ast }$	$x_{k+1}=x_k+\omega {𝐁}^{-1}{𝐀}^{\top }{𝐒}_k{\lambda }_k^{\approx }$ $=x_k-\omega \nabla f_{{𝐒}_k}\left(x_k\right)+\underset{{ϵ}_k}{\underbrace{\omega {𝐁}^{-1}{𝐀}^{\top }{𝐒}_k\left({\lambda }_k^{\approx }-{\lambda }_k^{\ast }\right)}}.$
SPM (8)	Exact computation of the projection ${\Pi }_{{ℒ}_{{𝐒}_k},𝐁}^{\ast }\left(x_k\right)=\arg {\min }_{x^{\prime }\in {ℒ}_{{𝐒}_k}}{\Vert x^{\prime }-x_k\Vert }_{𝐁}$	$x_{k+1}=\omega {\Pi }_{{ℒ}_{{𝐒}_k},𝐁}^{\approx }\left(x_k\right)+\left(1-\omega \right)x_k$ $=\omega {\Pi }_{{ℒ}_{{𝐒}_k}}^{𝐁}\left(x_k\right)+\left(1-\omega \right)x_k+\underset{{ϵ}_k}{\underbrace{\omega \left({\Pi }_{{ℒ}_{{𝐒}_k},𝐁}^{\approx }\left(x_k\right)-{\Pi }_{{ℒ}_{{𝐒}_k},𝐁}^{\ast }\left(x_k\right)\right)}}$
SDSA (11)	Exact computation of ${\lambda }_k^{\ast }$, where ${\lambda }_k^{\ast }\in \arg {\max }_{\lambda }D\left(y_k+{𝐒}_k\lambda \right)$.	$y_{k+1}=y_k+\omega {𝐒}_k{\lambda }_k^{\approx }=y_k+\omega {𝐒}_k{\lambda }_k^{\ast }+\underset{{ϵ}_k^d}{\underbrace{\omega {𝐒}_k\left({\lambda }_k^{\approx }-{\lambda }_k^{\ast }\right)}}$

Table 3. Table 3 : Popular inequalities with abbreviations and formulas.

Useful inequalities
Inequalities (Full names)	Abbreviations	Formula	Assumptions
Jensen Inequality	Jensen	$f\left[𝔼\left(x\right)\right]\le 𝔼\left[f\left(x\right)\right]$	f is convex
Conditioned Jensen	cond. Jensen	$f\left(𝔼\left[x|s\right]\right)\le 𝔼\left[f\left(x\right)|s\right]$	f is convex
Cauchy-Swartz (B-norm)	C.S	$\left|{⟨a,b⟩}_{𝐁}\right|\le {\Vert a\Vert }_{𝐁}{\Vert b\Vert }_{𝐁}$	$a,b\in {ℝ}^n$
Variance Inequality	V.I	${\left(𝔼\left[X\right]\right)}^2\le 𝔼\left[X^2\right]$	$X$ random variable

Table 4. Table 4 : Frequently used notation.

The Basics
$𝐀,b$	$m\times n$ matrix and $m\times 1$ vector defining the system $𝐀x=b$
$ℒ$	$\{x:𝐀x=b\}$ (solution set of the linear system)
$𝐁$	$n\times n$ symmetric positive definite matrix
${⟨x,y⟩}_{𝐁}$	$x^{\top }𝐁y$ ($𝐁$-inner product)
${\Vert x\Vert }_{𝐁}$	$\sqrt{{⟨x,x⟩}_{𝐁}}$ ($𝐁$-norm)
${𝐌}^{†}$	Moore-Penrose pseudoinverse of matrix $𝐌$
$𝐒$	a random real matrix with $m$ rows
$𝒟$	distribution from which matrix $𝐒$ is drawn ($𝐒\sim 𝒟$)
$𝐇$	$𝐒{\left({𝐒}^{\top }{\mathrm{𝐀𝐁}}^{-1}{𝐀}^{\top }𝐒\right)}^{†}{𝐒}^{\top }$
$𝐙$	${𝐀}^{\top }\mathrm{𝐇𝐀}$
$\mathrm{Range}\left(𝐌\right)$	range space of matrix $𝐌$
$\mathrm{Null}\left(𝐌\right)$	null space of matrix $𝐌$
$\mathbb{P}(\cdot )$	probability of an event
$𝔼[\cdot ]$	expectation
Projections
${\Pi }_{ℒ,𝐁}\left(x\right)$	projection of $x$ onto $ℒ$ in the $𝐁$-norm
${𝐁}^{-1}𝐙$	projection matrix, in the $𝐁$-norm, onto $\mathrm{Range}\left({𝐁}^{-1}{𝐀}^{\top }𝐒\right)$
Optimization
$𝒳$	set of minimizers of $f$
$x_{\ast }$	a point in $ℒ$
$f_{𝐒}$, $\nabla f_{𝐒}$, ${\nabla }^2{f}_{𝐒}$	stochastic function, its gradient and Hessian
${ℒ}_{𝐒}$	$\{x:{𝐒}^{\top }𝐀x={𝐒}^{\top }b\}$ (set of minimizers of $f_{𝐒}$)
$f$	$𝔼\left[f_{𝐒}\right]$
$\nabla f$	gradient of $f$ with respect to the $𝐁$-inner product
${\nabla }^{2}f$	${𝐁}^{-1}𝔼\left[𝐙\right]$ (Hessian of $f$ in the $𝐁$-inner product)
Eigenvalues
$𝐖$	${𝐁}^{-1/2}𝔼\left[𝐙\right]{𝐁}^{-1/2}$ (psd matrix with the same spectrum as ${\nabla }^{2}f$)
${\lambda }_1,\dots ,{\lambda }_n$	eigenvalues of $𝐖$
${\lambda }_{\max },{\lambda }_{\min }^{+}$	largest and smallest nonzero eigenvalues of $𝐖$
Algorithms
$\omega$	relaxation parameter / stepsize
${ϵ}_k$	Inexactness error
$q$	Inexactness parameter
$\rho$	$1-\omega \left(2-\omega \right){\lambda }_{\min }^{+}$

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x):={\mathbb{E}}_{\mathbf{S}\sim \mathcal{D}}[f_{\mathbf{S}}(x)],$$

$$\mathbf{A}x=b.$$

$$f_{\mathbf{S}}(x):=\frac{1}{2}\Vert \mathbf{A}x-b{\Vert }_{\mathbf{H}}^2=\frac{1}{2}(\mathbf{A}x-b{)}^{\top }\mathbf{H}(\mathbf{A}x-b),$$

$$x_{k+1}=x_k-\omega \nabla f_{{\mathbf{S}}_k}(x_k),$$

$$\text{emStochastic Newton Method (SNM)}:x_{k+1}=x_k-\omega ({\nabla }^2{f}_{{\mathbf{S}}_k}(x_k){)}^{{†}_{\mathbf{B}}}\nabla f_{{\mathbf{S}}_k}(x_k),$$

$$\text{emStochastic Proximal Point Method (SPPM)}:x_{k+1}=\arg \underset{x\in {\mathbb{R}}^n}{\min }\left\{f_{{\mathbf{S}}_k}(x)+\frac{1-\omega }{2\omega }\Vert x-x_k{\Vert }_{\mathbf{B}}^2\right\}.$$

$$\underset{x\in {\mathbb{R}}^n}{\min }P(x):=\frac{1}{2}\Vert x-x_0{\Vert }_{\mathbf{B}}^2\text{subject to}\mathbf{A}x=b.$$

$$x_{k+1}=\omega {\Pi }_{{\mathcal{L}}_{{\mathbf{S}}_k},\mathbf{B}}(x_k)+(1-\omega )x_k,$$

$$x_{k+1}={\Pi }_{{\mathcal{L}}_{\mathbf{S}},\mathbf{B}}(x_k),$$

$$\underset{y\in {\mathbb{R}}^m}{\max }D(y):=(b-\mathbf{A}{x}_0{)}^{\top }y-\frac{1}{2}\Vert {\mathbf{A}}^{\top }y{\Vert }_{{\mathbf{B}}^{-1}}^2.$$

$$y_{k+1}=y_k+\omega {\mathbf{S}}_k{\lambda }_k,$$

$${\lambda }_k={\left({\mathbf{S}}_k^{\top }\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{\mathbf{S}}_k\right)}^{†}{\mathbf{S}}_k^{\top }\left(b-\mathbf{A}(x_0+{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{y}_k)\right).$$

$$x_k=x(y_k)=x_0+{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{y}_k.$$

$$D(y_{\ast })-D(y)=\frac{1}{2}\Vert x(y_{\ast })-x(y){\Vert }_{\mathbf{B}}^2.$$

$${\Pi }_{\mathcal{L},\mathbf{B}}(x):=\arg \underset{x^{\prime }\in \mathcal{L}}{\min }\Vert x^{\prime }-x{\Vert }_{\mathbf{B}}=x-{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }(\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{)}^{†}(\mathbf{A}x-b).$$

$$\mathbf{Z}:={\mathbf{A}}^{\top }\mathbf{H}\mathbf{A}={\mathbf{A}}^{\top }\mathbf{S}({\mathbf{S}}^{\top }\mathbf{A}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }\mathbf{S}{)}^{†}{\mathbf{S}}^{\top }\mathbf{A}.$$

$$f_{\mathbf{S}}(x)=\frac{1}{2}(\mathbf{A}x-b{)}^{\top }\mathbf{H}(\mathbf{A}x-b)=\frac{1}{2}(x-x_{\ast }{)}^{\top }{\mathbf{Z}}_k(x-x_{\ast }),$$

$$\nabla f_{\mathbf{S}}(x)\stackrel{\text{eqref}eq:f_s}{=}{\mathbf{B}}^{-1}{\mathbf{A}}^{\top }\mathbf{H}(\mathbf{A}x-b)={\mathbf{B}}^{-1}{\mathbf{A}}^{\top }\mathbf{H}\mathbf{A}(x-x_{\ast })={\mathbf{B}}^{-1}\mathbf{Z}(x-x_{\ast }),$$

$$\mathbf{W}:={\mathbf{B}}^{-\frac{1}{2}}\mathbb{E}[\mathbf{Z}]{\mathbf{B}}^{-\frac{1}{2}},$$

$$x_{k+1}=x_k+\omega {\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{\mathbf{S}}_k{\lambda }_k^{\approx }=x_k-\omega \nabla f_{{\mathbf{S}}_k}(x_k)+\underset{{ϵ}_k}{\underbrace{\omega {\mathbf{B}}^{-1}{\mathbf{A}}^{\top }{\mathbf{S}}_k({\lambda }_k^{\approx }-{\lambda }_k^{\ast })}}.$$

$$x_{k+1}$$

$$x_{k+1}=x_k-\omega {\mathbf{A}}_{C:}^{\top }({\mathbf{A}}_{C:}{\mathbf{A}}_{C:}^{\top }{)}^{†}({\mathbf{A}}_{C:}{x}_k-b_C)+{ϵ}_k.$$

$$x_{k+1}=x_k-\omega {\mathbf{I}}_{:C}({\mathbf{I}}_{:C}^{\top }\mathbf{A}{\mathbf{I}}_{:C}{)}^{†}{\mathbf{I}}_{:C}^{\top }(\mathbf{A}{x}_k-b)+{ϵ}_k.$$

$$\mathbb{E}[\Vert {ϵ}_k{\Vert }_{\mathbf{B}}^2|x_k,{\mathbf{S}}_k]\le {\sigma }_k^2,$$

$$\mathbb{E}[\Vert {ϵ}_k{\Vert }_{\mathbf{B}}^2|x_k,{\mathbf{S}}_k]\le {\sigma }_k^2=q^2\Vert x_k-x_{\ast }{\Vert }_{\mathbf{B}}^2,$$

$$\mathbb{E}[\Vert {ϵ}_k{\Vert }_{\mathbf{B}}^2|x_k,{\mathbf{S}}_k]\le {\sigma }_k^2=2q^2{f}_{{\mathbf{S}}_k}(x_k),$$

$$\mathbb{E}[{⟨(\mathbf{I}-\omega {\mathbf{B}}^{-1}{\mathbf{Z}}_k)(x_k-x_{\ast }),{ϵ}_k⟩}_{\mathbf{B}}]=0.$$

$$\rho =1-\omega (2-\omega ){\lambda }_{\min }^{+}.$$

$$\mathbb{E}[\Vert x_k-x_{\ast }{\Vert }_{\mathbf{B}}]\le {\rho }^{k/2}\Vert x_0-x_{\ast }{\Vert }_{\mathbf{B}}+\sum _{i=0}^{k-1}{\rho }^{\frac{k-1-i}{2}}{\sigma }_i.$$

Full text

Keywords Inexact methods $\cdot$ Iteration complexity $\cdot$ Linear systems $\cdot$ Randomized block coordinate descent $\cdot$ Randomized block Kaczmarz $\cdot$ Stochastic gradient descent $\cdot$ Stochastic Newton method$\cdot$ Quadratic optimization $\cdot$ Convex optimization

**Mathematical Subject Classifications ** 68Q25 $\cdot$ 68W20 $\cdot$ 68W40 $\cdot$ 65Y20 $\cdot$ 90C15 $\cdot$ 90C20 $\cdot$ 90C25 $\cdot$ 15A06 $\cdot$ 15B52 $\cdot$ 65F10

1 Introduction

In the era of big data where data sets become continuously larger, randomized iterative methods become very popular and they are now playing major role in areas like numerical linear algebra, scientific computing and optimization. They are preferred mainly because of their cheap per iteration cost which leads to the improvement in terms of complexity upon classical results by orders of magnitude and to the fact that they can easily scale to extreme dimensions. However, a common feature of these methods is that in their update rule a particular subproblem needs to be solved exactly. In the case that the size of this problem is large, this step can be computationally very expensive. The purpose of this work is to reduce the cost of this step by incorporating inexact updates in the stochastic methods under study.

1.1 The Setting

In this paper we are interested to solve three closely related problems:

• •

Stochastic Quadratic Optimization Problem

• •

Best Approximation Problem

• •

Concave Quadratic Maximization Problem

We start by presenting the main connections and key relationships between these problems as well as popular randomized iterative methods (with exact updates) for solving each one of them.

Stochastic Optimization Problem:

We study the stochastic quadratic optimization problem

[TABLE]

first proposed in [46] for reformulating consistent linear systems

[TABLE]

In particular, problem (1) is defined by setting:

[TABLE]

where ${\bf H}$ is a random symmetric positive semi-definite matrix ${\bf H}:={\bf S}({\bf S}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S})^{\dagger}{\bf S}^{\top}$ that depends on three different matrices: the data matrix ${\bf A}\in\mathbb{R}^{m\times n}$ of the linear system (2), a random matrix $\mathbf{S}\in\mathbb{R}^{m\times q}\sim{\cal D}$ and on an $n\times n$ positive definite matrix ${\bf B}$ which defines the geometry of the space. Throughout the paper, ${\bf B}$ is used to define a $\mathbf{B}-$inner product in $\mathbb{R}^{n}$ via $\langle x,z\rangle_{\bf B}:=\langle{\bf B}x,z\rangle$ and an induced $\mathbf{B}-$norm, $\|x\|_{\bf B}:=(x^{\top}{\bf B}x)^{1/2}$. By $\dagger$ we denote the Moore-Penrose pseudoinverse.

The expectation in (1) is over random matrices ${\bf S}$ with $m$ rows (and arbitrary number of columns $q$, e.g., $q=1$) drawn from an arbitrary (user defined) distribution ${\cal D}$. The authors of [46] give necessary and sufficient conditions that distribution $D$ needs to be satisfied for the set of solutions of (1) to be equal to the set of solutions of the linear system (2); a property for which the term exactness was coined in (see Section 3 for more details on exactness).

In [46], problem (1) was solved via Stochastic Gradient Descent (SGD)111The gradient is computed with respect to the inner product $\langle{\bf B}x,y\rangle$.:

[TABLE]

and a linear rate of convergence was proved despite the fact that $f$ is not necessarily strongly convex, (1) is not a finite-sum problem and a fixed stepsize $\omega>0$ is used.

The stochastic optimization problem (1) has many unique characteristics mainly because it has constructed in a particular way in order to capture all the information of the linear system (2). For example it holds that $f_{{\bf S}}(x)=\tfrac{1}{2}\|\nabla f_{{\bf S}}(x)\|_{\bf B}^{2},$ and it can be proved that all eigenvalues of its Hessian matrix $\nabla^{2}f(x)$ are upper bounded by 1. Due to these specific characteristics, the update rules of seemingly different randomized iterative methods are identical. In particular the following methods for solving (1) have exactly the same behavior with SGD [46]:

[TABLE]

[TABLE]

In all methods $\omega>0$ is a fixed stepsize and ${\bf S}_{k}$ is sampled afresh in each iteration from distribution ${\cal D}$. See [46] for more insights into the reformulation (1), its properties and other equivalent reformulations (e.g., stochastic fixed point problem, probabilistic intersection problem, and stochastic linear system).

Best Approximation Problem and Sketch and Project Method:

In [46, 29], it has been shown that for the case of consistent linear systems with multiple solutions, SGD (and as a result SNM (5) and SPPM (6)) converges linearly to one particular minimizer of function $f$, the projection of the initial iterate $x_{0}$ onto the solution set of the linear system (2). This naturally leads to the best approximation problem:

[TABLE]

Unlike, the linear system (2) which is allowed to have multiple solutions, the best approximation problem has always (from its construction) a unique solution. For solving problem (7), the Sketch and Project Method (SPM):

[TABLE]

was analyzed in [18, 46]. Here, $\Pi_{{\cal L}_{{\bf S}_{k}},\mathbf{B}}(x_{k})$ denotes the projection of point $x_{k}$ onto ${\cal L}_{{\bf S}_{k}}=\{x\in\mathbb{R}^{n}\;:\;{\bf S}_{k}^{\top}{\bf A}x={\bf S}_{k}^{\top}b\}$ in the ${\bf B}$-norm. In the special case of unit stepsize ($\omega=1$) algorithm (8) simplifies to

[TABLE]

first proposed in [18]. The name Sketch and Project method is justified by the iteration structure which follows two steps: (i) Choose the sketched system ${\cal L}_{{\bf S}_{k}}:=\{x\;:\;{\bf S}^{\top}{\bf A}x={\bf S}^{\top}b\}$, (ii) Project the last iterate $x_{k}$ onto ${\cal L}_{{\bf S}_{k}}$. The Sketch and Project viewpoint will be useful later in explaining the natural interpretation of the proposed inexact update rules. (see Section 4.2).

Dual Problem and SDSA:

The Fenchel dual of (7) is the (bounded) unconstrained concave quadratic maximization problem

[TABLE]

Boundedness follows from consistency. It turns out that by varying ${\bf A},{\bf B}$ and $b$ (but keeping consistency of the linear system), the dual problem in fact captures all bounded unconstrained concave quadratic maximization problems [29].

A direct dual method for solving problem (10) was first proposed in [19]. The dual method—Stochastic Dual Subspace Ascent (SDSA)— updates the dual vectors $y_{k}$ as follows:

[TABLE]

where the random matrix ${\bf S}_{k}$ is sampled afresh in each iteration from distribution ${\cal D}$, and $\lambda_{k}$ is chosen in such a way to maximize the dual objective $D$: $\lambda_{k}\in\arg\max_{\lambda}D(y_{k}+{\bf S}_{k}\lambda)$. More specifically, SDSA is defined by picking the $\lambda_{k}$ with the smallest (standard Euclidean) norm. This leads to the formula:

[TABLE]

It can be proved, [19, 29], that the iterates $\{x_{k}\}_{k\geq 0}$ of the sketch and project method (8) arise as affine images of the iterates $\{y_{k}\}_{k\geq 0}$ of the dual method (11) as follows:

[TABLE]

In [19] the dual method was analyzed for the case of unit stepsize ($\omega=1$). Later in [29] the analysis extended to capture the cases of $\omega\in(0,2)$. Momentum variants of the dual method that provide further speed up have been also studied in [29].

An interesting property that holds between the suboptimalities of the Sketch and Project method and SDSA is that the dual suboptimality of $y$ in terms of the dual function values is equal to the primal suboptimality of $x(y)$ in terms of distance [19, 29]. That is,

[TABLE]

This simple to derive result (by combining the expression of the dual function $D(y)$ (10) and the equation (13)) gives for free the convergence analysis of SDSA, in terms of dual function suboptimality once the analysis of Sketch and Project is available (see Section 5).

1.2 Contributions

In this work we propose and analyze inexact variants of all previously mentioned randomized iterative algorithms for solving the stochastic optimization problem, the best approximation problem and the dual problem. In all of these methods, a certain potentially expensive calculation/operation needs to be performed in each step; it is this operation that we propose to be performed inexactly. For instance, in the case of SGD, it is the computation of the stochastic gradient $\nabla f_{{\bf S}_{k}}(x_{k})$, in the case of SPM is the computation of the projection $\Pi_{{\cal L}_{{\bf S}},\mathbf{B}}(x_{k})$, and in the case of SDSA it is the computation of the dual update ${\bf S}_{k}\lambda_{k}$.

We perform an iteration complexity analysis under an abstract notion of inexactness and also under a more structured form of inexactness appearing in practical scenarios. An inexact solution of these subproblems can be obtained much more quickly than the exact solution. Since in practical applications the savings thus obtained are larger than the increase in the number of iterations needed for convergence, our inexact methods can be dramatically faster.

Let us now briefly outline the rest of the paper:

In Section 2 we describe the subproblems and introduce two notions of inexactness (abstract and structured) that will be used in the rest of the paper. The Inexact Basic Method (iBasic) is also presented. iBasic is a method that simultaneously captures inexact variants of the algorithms (4), (5), (6) for solving the stochastic optimization problem (1) and algorithm (8) for solving the best approximation problem (7). It is an inexact variant of the Basic Method, first presented in [46], where the inexactness is introduced by the addition of an inexactness error $\epsilon_{k}$ in the original update rule. We illustrate the generality of iBasic by presenting popular algorithms that can be cast as special cases.

In Section 3 we establish convergence results of iBasic under general assumptions on the inexactness error $\epsilon_{k}$ of its update rule (see Algorithm 1). In this part we do not focus on any specific mechanisms which lead to inexactness; we treat the problem abstractly. However, such errors appear often in practical scenarios and can be associated with inaccurate numerical solvers, quantization, sparsification and compression mechanisms. In particular, we introduce several abstract assumptions on the inexactness level and describe our generic convergence results. For all assumptions we establish linear rate of decay of the quantity $\mathbb{E}[\|x_{k}-x_{*}\|_{{\bf B}}^{2}]$ (i.e. L2 convergence)444As we explain later, a convergence of the expected function values of problem 1 can be easily obtained as a corollary of L2 convergence..

Subsequently, in Section 4 we apply our general convergence results to a more structured notion of inexactness error and propose a concrete mechanisms leading to such errors. We provide theoretical guarantees for this method in situations when a linearly convergent iterative method (e.g., Conjugate Gradient) is used to solve the subproblem inexactly. We also highlight the importance of the dual viewpoint through a sketch-and-project interpretation.

In Section 5 we study an inexact variant of SDSA, which we called iSDSA, for directly solving the dual problem (10). We provide a correspondence between iBasic and iSDSA and we show that the random iterates of iBasic arise as affine images of iSDSA. We consider both abstract and structured inexactness errors and provide linearly convergent rates in terms of the dual function suboptimality ${\mathbb{E}\left[D(y_{*})-D(y_{0})\right]}$.

Finally, in Section 6 we evaluate the performance of the proposed inexact methods through numerical experiments and show the benefits of our approach on both synthetic and real datasets. Concluding remarks are given in Section 7.

A summary of the convergence results of iBasic under several assumptions on the inexactness error with pointers to the relevant theorems is available in Table 1. We highlight that similar convergence results can be also obtained for iSDSA in terms of the dual function suboptimality ${\mathbb{E}\left[D(y_{*})-D(y_{0})\right]}$ (check Section 5 for more details on iSDSA).

1.3 Notation

For convenience, a table of the most frequently used notation is included in the Appendix C. In particular, with boldface upper-case letters we denote matrices and $\mathbf{I}$ is the identity matrix. By ${\cal L}$ we denote the solution set of the linear system $\mathbf{A}x=b$. By ${\cal L}_{\mathbf{S}}$, where ${\bf S}$ is a random matrix, we denote the solution set of the sketched linear system $\mathbf{S}^{\top}\mathbf{A}x=\mathbf{S}^{\top}b$. In general, we use $\cdot^{*}$ to express the exact solution of a sub-problem and $\cdot^{\approx}$ to indicate its inexact variant. Unless stated otherwise, throughout the paper, $x_{*}$ is the projection of $x_{0}$ onto ${\cal L}$ in the ${\bf B}$-norm: $x_{*}=\Pi_{{\cal L},\mathbf{B}}(x_{0})$. An explicit formula for the projection of point $x$ onto set ${\cal L}$ is given by

[TABLE]

A formula for the projection onto ${\cal L}_{{\bf S}}=\{x\in\mathbb{R}^{n}\;:\;{\bf S}^{\top}{\bf A}x={\bf S}^{\top}b\}$ is obtained by replacing ${\bf A}$ and $b$ with ${\bf S}^{\top}{\bf A}$ and ${\bf S}^{\top}b$ respectively into the above equation. We denote this projection by $\Pi_{{\cal L}_{{\bf S}},\mathbf{B}}(x)$. We also write $[n]:=\{1,2,\dots,n\}$.

In order to keep the expression brief throughout the paper we define555In the $k^{th}$ iterate the expression becomes ${\bf Z}_{k}:={\bf A}^{\top}{\bf S}_{k}({\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k})^{\dagger}{\bf S}_{k}^{\top}{\bf A}$.:

[TABLE]

Using this matrix we can easily express important quantities related to the problems under study. For example the stochastic functions $f_{\bf S}$ of problem (1) can be expressed as

[TABLE]

In addition the gradient and the Hessian of $f_{\bf S}$ with respect to the ${\bf B}$ inner product are equal to

[TABLE]

and $\nabla^{2}f_{{\bf S}}(x)=\mathbf{B}^{-1}\mathbf{Z}$ [46]. Similarly the gradient and Hessian of the objective function $f$ of problem (1) are $\nabla f(x)={\bf B}^{-1}{\mathbb{E}\left[{\bf Z}\right]}(x-x_{*})$ and $\nabla^{2}f(x)={\bf B}^{-1}{\mathbb{E}\left[{\bf Z}\right]},$ respectively.

A key matrix in our analysis is

[TABLE]

which has the same spectrum with the matrix $\nabla^{2}f(x)$ but at the same time is symmetric and positive semi-definite666Note that matrix $\nabla^{2}f(x)$ is not symmetric but it is self-adjoint with respect to the $\mathbf{B}$-inner product.. We denote with $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{n}$ the $n$ eigenvalues of $\mathbf{W}$. With $\lambda_{\min}^{+}$ we indicate the smallest nonzero eigenvalue, and with $\lambda_{\max}=\lambda_{n}$ the largest eigenvalue. It was shown in [46] that $0\leq\lambda_{i}\leq 1$ for all $i\in[n]$.

2 Inexact update rules

In this section we start by explaining the key sub-problems that need to be solved exactly in the update rules of the previously described methods. We present iBasic, a method that solves problems (1) and (7) and we show how by varying the main parameters of the method we recover inexact variants of popular algorithms as special cases. Finally closely related work on inexact algorithms for solving different problems is also presented.

2.1 Expensive Sub-problems in Update Rules

Let us devote this subsection on explaining how the inexactness can be introduced in the current exact update rules of SGD777Note that SGD has identical updates to the Stochastic Newton and Stochastic proximal point method. Thus the inexactness can be added to these updates in similar way. (4), Sketch and Project (8) and SDSA (11) for solving the stochastic optimization, best approximation and the dual problem respectively. As we have shown these methods solve closely related problems and the key subproblems in their update rule are similar. However the introduction of inexactness in the update rule of each one of them can have different interpretation.

For example for the case of SGD for solving the stochastic optimization problem (1) (see also Section 4.1 and 4.2 for more details), if we define $\lambda_{k}^{*}=({\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k})^{\dagger}{\bf S}_{k}^{\top}(b-{\bf A}x_{k})$ then the stochastic gradient of function $f$ becomes \nabla f_{{\bf S}_{k}}(x_{k})\overset{\eqref{eq:grad_f_S}}{=}-\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k}\lambda_{k}^{*} and the update rule of SGD takes the form: $x_{k+1}=x_{k}+\omega{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda_{k}^{*}$. Clearly in this update the expensive part is the computation of the quantity $\lambda_{k}^{*}$ that can be equivalently computed to be the least norm solution of the smaller (in comparison to $\mathbf{A}x=b$) linear system ${\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda={\bf S}_{k}^{\top}(b-{\bf A}x_{k})$. In our work we are suggesting to use an approximation $\lambda_{k}^{\approx}$ of the exact solution and with this way avoid executing the possibly expensive step of the update rule. Thus the inexact update is taking the following form:

[TABLE]

Here $\epsilon_{k}$ denotes a more abstract notion of inexactness and it is not necessary to be always equivalent to the quantity $\omega\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k}(\lambda_{k}^{\approx}-\lambda_{k}^{*})$. It can be interpreted as an expression that acts as an perturbation of the exact update. In the case that $\epsilon_{k}$ has the above form we say that the notion of inexactness is structured. In our work we are interested in both the abstract and more structured notions of inexactness. We first present general convergence results where we require the error $\epsilon_{k}$ to satisfy general assumptions (without caring how this error is generated) and later we analyze the concept of structured inexactness by presenting algorithms where $\epsilon_{k}=\omega\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k}(\lambda_{k}^{\approx}-\lambda_{k}^{*})$.

In similar way, the expensive operation of SPM (8) is the exact computation of the projection $\Pi_{{\cal L}_{{\bf S}_{k}},{\bf B}}^{*}(x_{k})$. Thus we are suggesting to replace this step with an inexact variant and compute an approximation of this projection. The inexactness here can be also interpreted using both, the abstract $\epsilon_{k}$ error and its more structured version $\epsilon_{k}=\omega\left(\Pi_{{\cal L}_{{\bf S}_{k}},{\bf B}}^{\approx}(x_{k})-\Pi_{{\cal L}_{{\bf S}_{k}},{\bf B}}^{*}(x_{k})\right)$. At this point, observe that, by using the expression (15) the structure of the $\epsilon_{k}$ in SPM and SGD has the same form.

In the SDSA the expensive subproblem in the update rule is the computation of the $\lambda_{k}^{*}$ that satisfy $\lambda_{k}^{*}\in\arg\max_{\lambda}D(y_{k}+{\bf S}_{k}\lambda)$. Using the definition of the dual function (10) this value can be also computed by evaluating the least norm solution of the linear system ${\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda={\bf S}_{k}^{\top}\left(b-\mathbf{A}(x_{0}+\mathbf{B}^{-1}\mathbf{A}^{\top}y_{k}\right))$. Later in Section 5 we analyze both notions of inexactness (abstract and more structured) for inexact variants of SDSA.

Table 2 presents the key sub-problem that needs to be solved in each algorithm as well as the part where the inexact error is appeared in the update rule.

2.2 The Inexact Basic Method

In each iteration of the all aforementioned exact methods a sketch matrix $\mathbf{S}\sim{\cal D}$ is drawn from a given distribution and then a certain subproblem is solved exactly to obtain the next iterate. The sketch matrix $\mathbf{S}\in\mathbb{R}^{m\times q}$ requires to have $m$ rows but no assumption on the number of columns is made which means that the number of columns $q$ allows to vary through the iterations and it can be very large. The setting that we are interested in is precisely that of having such large random matrices $\mathbf{S}$. In these cases we expect that having approximate solutions of the subproblems will be beneficial.

Recently randomized iterative algorithms that requires to solve large subproblems in each iteration have been extensively studied and it was shown that are really beneficial when they compared to their single coordinates variants ($\mathbf{S}\in\mathbb{R}^{m\times 1}$) [34, 35, 44, 27]. However, in theses cases the evaluation of an exact solution for the suproblem in the update rule can be computationally very expensive. In this work we propose and analyze inexact variants by allowing to solve the subproblem that appear in the update rules of the stochastic methods, inexactly. In particular, following the convention established in [46] of naming the main algorithm of the paper Basic method we propose the inexact Basic method (iBasic) (Algorithm 1).

The $\epsilon_{k}$ in the update rule of the method represents the abstract inexactness error described in Subsection 2.1. Note that, iBasic can have several equivalent interpretations. This allow as to study the methods (4),(5),(6) for solving the stochastic optimization problem and the sketch and project method (8) for the best approximation problem in a single algorithm only. In particular iBasic can be seen as inexact stochastic gradient descent (iSGD) with fixed stepsize applied to (1). From (17), $\nabla f_{{\bf S}_{k}}(x_{k})={\bf B}^{-1}{\bf A}^{\top}{\bf H}_{k}({\bf A}x_{k}-b)$ and as a result the update rule of iBasic can be equivalently written as: $x_{k+1}=x_{k}-\omega\nabla f_{{\bf S}_{k}}(x_{k})+\epsilon_{k}.$ In the case of the best approximation problem (7), iBasic can be interpreted as inexact Sketch and Project method (iSPM) as follows:

[TABLE]

For the dual problem (10) we devote Section 5 for presenting an inexact variant of the SDSA (iSDSA) and analyze its convergence using the rates obtained for the iBasic in Sections 3 and 4.

2.3 General Framework and Further Special Cases

The proposed inexact methods, iBasic (Algorithm 1) and iSDSA (Section 5), belong in the general sketch and project framework, first proposed from Gower and Richtarik in [18] for solving consistent linear systems and where a unified analysis of several randomized methods was studied. This interpretation of the algorithms allow us to recover a comprehensive array of well-known methods as special cases by choosing carefully the combination of the main parameters of the algorithms.

In particular, the iBasic has two main parameters (besides the stepsize $\omega>0$ of the update rule). These are the distribution ${\cal D}$ from which we draw random matrices $\mathbf{S}$ and the positive definite matrix $\mathbf{B}\in\mathbb{R}^{n\times n}$. By choosing carefully combinations of the parameters ${\cal D}$ and $\mathbf{B}$ we can recover several existing popular algorithms as special cases of the general method. For example, special cases of the exact Basic method are the Randomized Kaczmarz, Randomized Gaussian Kaczmarz888Special case of the iBasic, when the random matrix $\mathbf{S}$ is chosen to be a Gaussian vector with mean $0\in R^{m}$ and a positive definite covariance matrix $\Sigma\in\mathbb{R}^{m\times m}$. That is $\mathbf{S}\sim N(0,\Sigma)$ [18, 29]., Randomized Coordinate Descent and their block variants. For more details about the generality of the sketch and project framework and further algorithms that can be cast as special cases of the analysis we refer the interested reader to Section 3 of [18] and Section 7 of [29]. Here we present only the inexact update rules of two special cases that we will later use in the numerical evaluation.

*Special Cases: * Let us define with $\mathbf{I}_{:C}$ the column concatenation of the $m\times m$ identity matrix indexed by a random subset $C$ of $[m]$.

• •

Inexact Randomized Block Kaczmarz (iRBK): Let $\mathbf{B}=\mathbf{I}$ and let pick in each iteration the random matrix $\mathbf{S}=\mathbf{I}_{:C}\sim{\cal D}$. In this setup the update rule of the iBasic simplifies to

[TABLE]

• •

*Inexact Randomized Block Coordinate Descent (iRBCD)*999In the setting of solving linear systems Randomized Coordinate Descent is known also as Gauss-Seidel method. Its block variant can be also interpret as randomized coordinate Newton method (see [42]).: If the matrix $\mathbf{A}$ of the linear system is positive definite then we can choose $\mathbf{B}=\mathbf{A}$. Let also pick in each iteration the random matrix $\mathbf{S}=\mathbf{I}_{:C}\sim{\cal D}$. In this setup the update rule of the iBasic simplifies to

[TABLE]

For more papers related to Kaczmarz method (randomized, greedy, cyclic update rules) we refer the interested reader to [23, 28, 38, 5, 37, 39, 9, 33, 34, 13, 31, 59, 35, 50]. For the coordinate descent method (a.k.a Gauss-Seidel for linear systems) and its block variant, Randomized Block Coordinate Descent we suggest [25, 36, 44, 45, 40, 41, 43, 7, 24, 14, 1, 54].

2.4 Other Related Work on Inexact Methods

One of the current trends in the large scale optimization problems is the introduction of inexactness in the update rules of popular deterministic and stochastic methods. The rational behind this is that an approximate/inexact step can often computed very efficiently and can have significant computational gains compare to its exact variants.

In the area of deterministic algorithms, the inexact variant of the full gradient descent method, $x_{k+1}=x_{k}-\omega_{k}[\nabla f(x_{k})+\epsilon_{k}]$, has received a lot of attention [49, 11, 51, 16, 32]. It has been analyzed for the cases of convex and strongly convex functions under several meaningful assumptions on the inexactness error $\epsilon_{k}$ and its practical benefit compared to the exact gradient descent is apparent. For further deterministic inexact methods check [10] for Inexact Newton methods, [52, 47] for Inexact Proximal Point methods and [3] for Inexact Fixed point methods.

In the recent years, with the explosion that happens in areas like machine learning and data science inexactness enters also the updating rules of several stochastic optimization algorithms and many new methods have been proposed and analyzed.

In the large scale setting, stochastic optimization methods are preferred mainly because of their cheap per iteration cost (compared to their deterministic variants), their property to scale to extreme dimensions and their improved theoretical complexity bounds. In areas like machine learning and data science, where the datasets become larger rapidly, the development of faster and efficient stochastic algorithms is crucial. For this reason, inexactness has recently introduced to the update rules of several stochastic optimization algorithms and new methods have been proposed and analyzed. One of the most interesting work on inexact stochastic algorithms appears in the area of second order methods. In particular on inexact variants of the Sketch-Newton method and subsampled Newton Method for minimize convex and non-convex functions [48, 2, 4, 56, 57, 58]. Note that our results are related also with this literature since our algorithm can be seen as inexact stochastic Newton method (see equation (5)). To the best or our knowledge our work is the first that provide convergence analysis of inexact stochastic proximal point methods (equation (6)) in any setting. From numerical linear algebra viewpoint inexact sketch and project methods for solving the best approximation problem and its dual problem where also never analyzed before.

As we already mentioned our framework is quite general and many algorithms, like iRBK (21) and iRBCD (22) can be cast as special cases. As a result, our general convergence analysis includes the analysis of inexact variants of all of these more specific algorithms as special cases. In [34] an analysis of the exact randomized block Kacmzarz method has been proposed and in the experiments an inexact variant was used to speedup the method. However, no iteration complexity results were presented for the inexact variant and both the analysis and numerical evaluation have been made for linear systems with full rank matrices that come with natural partition of the rows (this is a much more restricted case than the one analyzed in our setting). For inexact variants of the randomized block coordinate descent algorithm in different settings than ours we suggest [53, 15, 6, 12].

Finally an analysis of approximate stochastic gradient descent for solving the empirical risk minimization problem using quadratic constraints and sequential semi-definite programs has been presented in [22].

3 Convergence Results Under General Assumptions

In this section we consider scenarios in which the inexactness error $\epsilon_{k}$ can be controlled, by specifying a per iteration bound $\sigma_{k}$ on the norm of the error. In particular, by making different assumptions on the bound $\sigma_{k}$ we derive general convergence rate results. Our focus is on the abstract notion of inexactness described in Section 2.1 and we make no assumptions on how this error is generated.

An important assumption that needs to be hold in all of our results is exactness. A formal presentation is presented below. We state it here and we highlight that is a requirement for all of our convergence results (exactness is also required in the analysis of the exact variants [46]).

Exactness.

Note that $f_{\bf S}$ is a convex quadratic, and that $f_{\bf S}(x)=0$ whenever $x\in{\cal L}:=\{x\;:\;{\bf A}x=b\}$. However, $f_{\bf S}$ can be zero also for points $x$ outside of ${\cal L}$. Clearly, $f(x)$ is nonnegative, and $f(x)=0$ for $x\in{\cal L}$. However, without further assumptions, the set of minimizers of $f$ can be larger than ${\cal L}$. The exactness assumption ensures that this does not happen. For necessary and sufficient conditions for exactness, we refer the reader to [46]. Here it suffices to remark that a sufficient condition for exactness is to require ${\mathbb{E}\left[{\bf H}\right]}$ to be positive definite. This is easy to see by observing that $f(x)={\mathbb{E}\left[f_{\bf S}(x)\right]}=\tfrac{1}{2}\|{\bf A}x-b\|^{2}_{{\mathbb{E}\left[{\bf H}\right]}}$. In other words, if ${\cal X}=\text{argmin}f(x)$ is the solution set of the stochastic optimization problem (1) and ${\cal L}=\{x:\mathbf{A}x=b\}$ the solution set of the linear system (2) then the notion of exactness is captured by: ${\cal X}={\cal L}$

3.1 Assumptions on Inexactness Error

In the convergence analysis of iBasic the following assumptions on the inexactness error are used. We note that Assumptions Assumption 1${a}$, Assumption 1${b}$ and Assumption 1${c}$ are special cases of Assumption Assumption 1. Moreover Assumption Assumption 2 is algorithmic dependent and can hold in addition of any of the other four assumptions. In our analysis, depending on the result we aim at, we will require either one of the first four Assumptions to hold by itself, or to hold together with Assumption Assumption 2. We will always assume exactness.

In all assumptions the expectation on the norm of error ($\|\epsilon_{k}\|^{2}$) is conditioned on the value of the current iterate $x_{k}$ and the random matrix $\mathbf{S}_{k}$. Moreover it is worth to mention that for the convergence analysis we never assume that the inexactness error has zero mean, that is $\mathbb{E}[\epsilon_{k}]=0$.

Assumption 1*.*

[TABLE]

where the upper bound $\sigma_{k}$ is a sequence of random variables (that can possibly depends on both the value of the current iterate $x_{k}$ and the choice of the random $\mathbf{S}_{k}$ at the $k^{th}$ iteration).

The following three assumptions on the sequence of upper bounds are more restricted however as we will later see allow us to obtain stronger and more controlled results.

Assumption 1${a}$.

[TABLE]

where the upper bound $\sigma_{k}\in\mathbb{R}$ is a sequence of real numbers.

Assumption 1${b}$.

[TABLE]

where the upper bound is a special sequence that depends on a non-negative inexactness parameter $q$ and the distance to the optimal value $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$.

Assumption 1${c}$.

[TABLE]

where the upper bound is a special sequence that depends on a non-negative inexactness parameter $q$ and the value of the stochastic function $f_{\mathbf{S}_{k}}$ computed at the iterate $x_{k}$.

Finally the next assumption is more algorithmic oriented. It holds in cases where the inexactness error $\epsilon_{k}$ in the update rule is chosen to be orthogonal with respect to the $\mathbf{B}$-inner product to the vector $\Pi_{{\cal L}_{\mathbf{S}_{k}},{\mathbf{B}}}(x_{k})-x_{*}=({\bf I}-\omega{\bf B}^{-1}{\bf Z}_{k})(x_{k}-x_{*})$. This statement may seem odd at this point but its usefulness will become more apparent in the next section where inexact algorithms with structured inexactness error will be analyzed. As it turns out, in the case of structured inexactness error (Algorithm 2) this assumption is satisfied.

Assumption 2*.*

[TABLE]

3.2 Convergence Results

In this section we present the analysis of the convergence rates of iBasic by assuming several combination of the previous presented assumptions.

All convergence results are described only in terms of convergence of the iterates $x_{k}$, that is $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$, and not the objective function values $f(x_{k})$. This is sufficient, because by $f(x)\leq\frac{\lambda_{\rm max}}{2}\|x-x_{*}\|^{2}_{\mathbf{B}}$ (see Lemma 10) we can directly deduce a convergence rate for the function values.

The exact Basic method (Algorithm 1 with $\epsilon_{k}=0$), has been analyzed in [46] and it was shown to converge with $\mathbb{E}[\|x_{k}-x_{*}\|^{2}_{{\bf B}}]\leq\rho^{k}\|x_{0}-x_{*}\|^{2}_{{\bf B}}$ where $\rho=1-\omega(2-\omega)\lambda_{\min}^{+}$. Our analysis of iBasic is more general and includes the convergence of the exact Basic method as special case when we assume that the upper bound is $\sigma_{k}=0,\quad\forall k\geq 0$. For brevity, in he convergence analysis results of this manuscript we also use

[TABLE]

Let us start by presenting the convergence of iBasic when only Assumption Assumption 1${a}$ holds for the inexactness error.

Theorem 1**.**

Let assume exactness and let $\{x_{k}\}_{k=0}^{\infty}$ be the iterates produced by iBasic with $\omega\in(0,2)$. Set $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$ and consider the error $\epsilon_{k}$ be such that it satisfies Assumption Assumption 1${a}$. Then,

[TABLE]

Proof.

See Appendix B.1. ∎

Corollary 1**.**

In the special case that the upper bound $\sigma_{k}$ in Assumption Assumption 1${a}$ is fixed, that is $\sigma_{k}=\sigma$ for all $k>0$ then inequality (28) of Theorem 1 takes the following form:

[TABLE]

This means that we obtain a linear convergence rate up to a solution level that is proportional to the upper bound $\sigma$101010Several similar more specific assumptions can be made for the upper bound $\sigma_{k}$. For example if the upper bound satisfies $\sigma_{k}=\sigma^{k}$ with $\sigma\in(0,1)$ for all $k>0$ then it can be shown that $C\in(0,1)$ exist such that inequality (28) of Theorem 1 takes the form: $\mathbb{E}[\|x_{k}-x_{*}\|_{{\bf B}}]\leq O(C^{k})$ (see [51, 16] for similar results). .

Proof.

See Appendix B.2. ∎

Inspired from [16], let us now analyze iBasic using the sequence of upper bounds that described in Assumption Assumption 1${b}$. This construction of the upper bounds allows us to obtain stronger and more controlled results. In particular using the upper bound of Assumption Assumption 1${b}$ the sequence of expected errors converge linearly to the exact $x_{*}$ (not in a potential neighborhood like the previous result). In addition Assumption Assumption 1${b}$ guarantees that the distance to the optimal solution reduces with the increasing of the number of iterations. However for this stronger convergence a bound for $\lambda_{\min}^{+}$ is required, a quantity that in many problems is unknown to the user or intractable to compute. Nevertheless, there are cases that this value has a close form expression and can be computed before hand without any further cost. See for example [27, 30, 26, 21] where methods for solving the average consensus were presented and the value of $\lambda_{\min}^{+}$ corresponds to the algebraic connectivity of the network under study.

Theorem 2**.**

Assume exactness. Let $\{x_{k}\}_{k=0}^{\infty}$ be the iterates produced by iBasic with $\omega\in(0,2)$. Set $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$ and consider the inexactness error $\epsilon_{k}$ be such that it satisfies Assumption Assumption 1${b}$, with $0\leq q<1-\sqrt{\rho}$. Then

[TABLE]

Proof.

See Appendix B.3. ∎

At Theorem 2, to guarantee linear convergence the inexact parameter $q$ should live in the interval $\left[0,1-\sqrt{\rho}\right)$. In particular, $q$ is the parameter that controls the level of inexactness of Algorithm 1. Not surprisingly the fastest convergence rate is obtained when $q=0$; in such case the method becomes equivalent with its exact variant and the convergence rate simplifies to $\rho=1-\omega(2-\omega)\lambda_{\min}^{+}$. Note also that similar to the exact case the optimal convergence rate is obtained for $\omega=1$ [46].

Moreover, the upper bound $\sigma_{k}$ of Assumption Assumption 1${b}$ depends on two important quantities, the $\lambda_{\min}^{+}$ (through the upper bound of the inexactness parameter $q$) and the distance to the optimal solution $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$. Thus, it can have natural interpretation. In particular the inexactness error is allowed to be large either when the current iterate is far from the optimal solution ($\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$ large) or when the problem is well conditioned and $\lambda_{\min}^{+}$ is large. In the opposite scenario, when we have ill conditioned problem or we are already close enough to the optimum $x_{*}$ we should be more careful and allow less errors to the updates of the method.

In the next theorem we provide the complexity results of iBasic in the case that the Assumption Assumption 2 is satisfied combined with one of the previous assumptions.

Theorem 3**.**

Let assume exactness and let $\{x_{k}\}_{k=0}^{\infty}$ be the iterates produced by iBasic with $\omega\in(0,2)$. Set $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$. Let also assume that the inexactness error $\epsilon_{k}$ be such that it satisfies Assumption Assumption 2. Then:

1. (i)

If Assumption Assumption 1 holds:

[TABLE]

where $\bar{\sigma}_{i}^{2}=\mathbb{E}[\sigma_{i}^{2}],\forall i\in[k-1].$ 2. (ii)

If Assumption Assumption 1${b}$ holds with $q\in\left(0,\sqrt{\rho}\right)$:

[TABLE] 3. (iii)

If Assumption Assumption 1${c}$ holds with $q\in\left(0,\sqrt{\omega(2-\omega)}\right)$:

[TABLE]

Proof.

See Appendix B.4. ∎

Remark 1**.**

In the case that Assumptions Assumption 1${a}$ and Assumption 2 hold simultaneously, the convergence of iBasic is similar to (31) but in this case ${\bar{\sigma}}^{2}_{i}=\sigma_{i}^{2},\,\forall i\in[k-1]$ (due to Assumption Assumption 1${a}$, $\sigma_{k}\in\mathbb{R}$ is a sequence of real numbers). In addition, note that for $q\in(0,\min\{\sqrt{\rho},1-\sqrt{\rho}\})$ having Assumption Assumption 2 on top of Assumption Assumption 1${b}$ leads to improvement of the convergence rate. In particular, from Theorem 2, iBasic converges with rate $(\sqrt{\rho}+q)^{2}=\rho+q^{2}+2\sqrt{\rho}q$ while having both assumptions this is simplified to the faster $\rho+q^{2}$ (32).

4 iBasic with Structured Inexactness Error

Up to this point, the analysis of iBasic was focused in more general abstract cases where the inexactness error $\epsilon_{k}$ of the update rule satisfies several general assumptions. In this section we are focusing on a more structured form of inexactness error and we provide convergence analysis in the case that a linearly convergent algorithm is used for the computation of the expensive key subproblem of the method.

4.1 Linear System in the Update Rule

As we already mentioned in Section 2.1 the update rule of the exact Basic method (Algorithm 1 with $\epsilon_{k}=0$) can be expressed as $x_{k+1}=x_{k}+\omega{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda_{k}^{*}$, where $\lambda_{k}^{*}=({\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k})^{\dagger}{\bf S}_{k}^{\top}(b-{\bf A}x_{k})$.

Using this expression the exact Basic method can be equivalently interpreted as the following two step procedure:

Find the least norm solution111111We are precisely looking for the least norm solution of the linear system $\mathbf{M}_{k}\lambda=d_{k}$ because this solution can be written down in a compact way using the Moore-Penrose pseudoinverse. This is equivalent with the expression that appears in our update: $\lambda_{k}^{*}=({\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k})^{\dagger}{\bf S}_{k}^{\top}(b-{\bf A}x_{k})=\mathbf{M}_{k}^{\dagger}d_{k}$. However it can be easily shown that the method will still converge with the same rate of convergence even if we choose any other solution of the linear system $\mathbf{M}_{k}\lambda=d_{k}$. of $\underbrace{{\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}}_{\mathbf{M}_{k}}\lambda=\underbrace{{\bf S}_{k}^{\top}(b-{\bf A}x_{k})}_{d_{k}}$. That is find $\lambda_{k}^{*}=\arg\min_{\lambda\in{\cal Q}_{k}}\|\lambda\|$ where ${\cal Q}_{k}=\left\{\lambda\in\mathbb{R}^{q}:\mathbf{M}_{k}\lambda=d_{k}\right\}$. 2. 2.

Compute the next iterate: $x_{k+1}=x_{k}+\omega{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda_{k}^{*}.$

In the case that the random matrix $\mathbf{S}_{k}$ is large (this is the case that we are interested in), solving exactly the linear system $\mathbf{M}_{k}\lambda=d_{k}$ in each step can be prohibitively expensive. To reduce this cost we allow the inner linear system $\mathbf{M}_{k}\lambda=d_{k}$ to be solved inexactly using an iterative method. In particular we propose and analyze the following inexact algorithm:

For the computation of the inexact solution of the linear system (34) any known iterative method for solving general linear systems can be used. In our analysis we focus on linearly convergent methods. For example based on the properties of the linear system (34), conjugate gradient (CG) or sketch and project method (SPM) can be used for the execution of step 3. In these cases, we name Algorithm 2, InexactCG and *InexactSP * respectively.

It is known that the classical CG can solve linear systems with positive definite matrices. In our approach matrix $\mathbf{M}_{k}$ is positive definite only when the original linear system $\mathbf{A}x=b$ has full rank matrix $\mathbf{A}$. On the other side SPM can solve any consistent linear system and as a result can solve the inner linear system $\mathbf{M}_{k}\lambda_{k}=d_{k}$ without any further assumption on the original linear system. In this case, one should be careful because the system has no unique solution. We are interested to find the least norm solution of $\mathbf{M}_{k}\lambda_{k}=d_{k}$ which means that the starting point of the sketch and project at the $k^{th}$ iteration should be always $\lambda_{k}^{0}=0$. Recall that any special case of the sketch and project method (Section 2.3) solves the best approximation problem.

Let us now define $\lambda_{k}^{r}$ to be the approximate solution $\lambda_{k}^{\approx}$ of the $q\times q$ linear system (34) obtained after $r$ steps of the linearly convergent iterative method. Using this, the update rule of Algorithm 2, takes the form:

[TABLE]

Remark 2**.**

The update rule (35) of Algorithm 2 is equivalent to the update rule of iBasic (Algorithm 1) when the error $\epsilon_{k}$ is chosen to be,

[TABLE]

This is precisely the connection between the abstract and more concrete/structured notion of inexactness that first presented in Table 2.

Let us now define a Lemma that is useful for the analysis of this section and it verifies that Algorithm 2 with unit stepsize satisfies the general Assumption Assumption 2 presented in Section 3.1.

Lemma 4**.**

Let us denote $x_{k}^{*}=\Pi_{{\cal L}_{\mathbf{S}_{k}},\mathbf{B}}(x_{k})$ the projection of $x_{k}$ onto ${\cal L}_{\mathbf{S}_{k}}$ in the ${\bf B}$-norm and $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$. Let also assume that $\omega=1$ (unit stepsize). Then for the updates of Algorithm 2 it holds that:

[TABLE]

Proof.

Note that $x_{k}^{*}-x_{*}=x_{k}-\nabla f_{{\bf S}_{k}}(x_{k})-x_{*}\in Null(\mathbf{S}_{k}^{\top}\mathbf{A})$ . Moreover $\epsilon_{k}\overset{\eqref{specEpsilon}}{=}\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k}(\lambda_{k}^{r}-\lambda_{k}^{*})\in Range(\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k})$. From the knowledge that the null space of an arbitrary matrix is the orthogonal complement of the range space of its transpose we have that $Null(\mathbf{S}_{k}^{\top}\mathbf{A})$ is orthogonal with respect to the $\mathbf{B}$-inner product to $Range(\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k})$. This completes the proof (see Figure 1 for the graphical interpretation). ∎

4.2 Sketch and Project Interpretation

Let us now give a different interpretation of the inexact update rule of Algorithm 2 using the sketch and project approach. That will make us appreciate more the importance of the dual viewpoint and make clear the connection between the primal and dual methods.

Recall that in the special case of unit stepsize (see equation (9)) the exact sketch and project method perform updates of the form:

[TABLE]

That is, a sketched system ${\bf S}^{\top}{\bf A}x={\bf S}^{\top}b$ is first chosen and then a the next iterate is computed by making a projection of the current iterate $x_{k}$ onto this system.

In general, execute a projection step is one of the most common task in numerical linear algebra/optimization literature. However in the large scale setting even this task can be prohibitively expensive and it can be difficult to execute inexactly. For this reason we suggest to move to the dual space where the inexactness can be easily controlled.

Observe that the update rule of equation (38) has the same structure as the best approximation problem (7) where the linear system under study is the sketched system $\mathbf{S}_{k}^{\top}{\bf A}x=\mathbf{S}_{k}^{\top}b$ and the starting point is the current iterate $x_{k}$. Hence we can easily compute its dual:

[TABLE]

where $\lambda\in\mathbb{R}^{q}$ is the dual variable. The $\lambda_{k}^{*}$ (possibly more than one) that solves the dual problem in each iteration $k$, is the one that satisfies $\nabla D_{k}(\lambda_{k}^{*})=0$. By computing the derivative this is equivalent with finding the $\lambda$ that satisfies the linear system ${\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda={\bf S}_{k}^{\top}(b-{\bf A}x_{k})$. This is the same linear system we desire to solve inexactly in Algorithm 2. Thus, computing an inexact solution $\lambda_{k}^{\approx}$ of the linear system is equivalent with computing an inexact solution of the dual problem (39). Then by using the affine mapping (13) that connects the primal and the dual spaces we can also evaluate an inexact solution of the original primal problem (38).

The following result relates the inexact levels of these quantities. In particular it shows that dual suboptimality of $\lambda_{k}$ in terms of dual function values is equal to the distance of the dual values $\lambda_{k}$ in the $\mathbf{M}_{k}$-norm.

Lemma 5**.**

Let us define $\lambda_{k}^{*}\in\mathbb{R}^{q}$ be the exact solution of the linear system ${\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}\lambda={\bf S}_{k}^{\top}(b-{\bf A}x_{k})$ or equivalently of dual problem (39). Let us also denote with $\lambda_{k}^{\approx}\in\mathbb{R}^{q}$ the inexact solution. Then:

[TABLE]

Proof.

[TABLE]

where in the second equality we use equation (13) to connect the optimal solutions of (38) and (39) and obtain $[\mathbf{S}_{k}^{\top}b-\mathbf{S}_{k}^{\top}\mathbf{A}x_{k}]^{\top}=(\lambda_{k}^{*})^{\top}{\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k}.$ ∎

4.3 Complexity Results

In this part we analyze the performance of Algorithm 2 when a linearly convergent iterative method is used for solving inexactly the linear system (34) in step 3 of Algorithm 2 . We denote with $\lambda_{k}^{r}$ the approximate solution of the linear system after we run the iterative method for $r$ steps.

Before state the main convergence result let us present a lemma that summarize some observations that are true in our setting.

Lemma 6**.**

Let $\lambda_{k}^{*}=({\bf S}_{k}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S}_{k})^{\dagger}{\bf S}_{k}^{\top}(b-{\bf A}x_{k})$ be the exact solution and $\lambda_{k}^{r}$ be approximate solution of the linear system (34). Then, $\|\lambda_{k}^{*}\|^{2}_{\mathbf{M}_{k}}=2f_{\mathbf{S}_{k}}(x_{k})$ and $\|\epsilon_{k}\|_{{\bf B}}^{2}=\|\lambda_{k}^{r}-\lambda_{k}^{*}\|_{\mathbf{M}_{k}}^{2}$.

Proof.

[TABLE]

Moreover,

[TABLE]

∎

Theorem 7**.**

Let us assume that for the computation of the inexact solution of the linear system (34) in step 3 of Algorithm 2, a linearly convergent iterative method is chosen such that 121212In the case that deterministic iterative method is used, like CG, we have that $\|\lambda_{k}^{r}-\lambda_{k}^{*}\|_{\mathbf{M}_{k}}^{2}\leq\rho_{\mathbf{S}_{k}}^{r}\|\lambda_{k}^{0}-\lambda_{k}^{*}\|^{2}_{\mathbf{M}_{k}}$ which is also true in expectation:

[TABLE]

where $\lambda_{k}^{0}=0$ for any $k>0$ and $\rho_{\mathbf{S}_{k}}\in(0,1)$ for every choice of $\mathbf{S}_{k}\sim{\cal D}$. Let exactness hold and let $\{x_{k}\}_{k=0}^{\infty}$ be the iterates produced by Algorithm 2 with unit stepsize ($\omega=1$). Set $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$. Suppose further that there exists a scalar $\theta<1$ such that with probability 1, $\rho_{\mathbf{S}_{k}}\leq\theta$. Then, Algorithm 2 converges linearly with:

[TABLE]

Proof.

Theorem 7 can be interpreted as corollary of the general Theorem 3(iii). Thus, it is sufficient to show that Algorithm 2 satisfies the two Assumptions Assumption 1${c}$ and Assumption 2. Firstly, note that from Lemma 4, Assumption Assumption 2 is true. Moreover,

[TABLE]

which means that Assumption Assumption 1${c}$ also holds with $q=\theta^{r/2}\in(0,1)$. This completes the proof. ∎

Having present the main result of this section let us now state some remarks that will help understand the convergence rate of the last Theorem.

Remark 3**.**

From its definition $\theta^{r}\in(0,1)$ and as a result $\left(1-\theta^{r}\right)\lambda_{\min}^{+}\leq\lambda_{\min}^{+}$. This means that the method converges linearly but always with worst rate than its exact variant.

Remark 4**.**

Let as assume that $\theta$ is fixed. Then as the number of iterations in step 3 of the algorithm ($r\rightarrow\infty$) increasing $(1-\theta^{r})\rightarrow 1$ and as a result the method behaves similar to the exact case.

Remark 5**.**

The $\lambda_{\min}^{+}$ depends only on the random matrices $\mathbf{S}\sim{\cal D}$ and to the positive definite matrix $\mathbf{B}$ and is independent to the iterative process used in step 3. The iterative process of step 3 controls only the parameter $\theta$ of the convergence rate.

Remark 6**.**

Let us assume that we run Algorithm 2 two separate times for two different choices of the linearly convergence iterative method of step 3. Let also assume that the distribution ${\cal D}$ of the random matrices and the positive definite matrix $\mathbf{B}$ are the same for both instances and that for step 3 the iterative method run for $r$ steps for both algorithms. Let assume that $\theta_{1}<\theta_{2}$ then we have that $\rho_{1}=1-\left(1-\theta_{1}^{r}\right)\lambda_{\min}^{+}<1-\left(1-\theta_{2}^{r}\right)\lambda_{\min}^{+}=\rho_{2}$. This means in the case that $\theta$ is easily computable, we should always prefer the inexact method with smaller $\theta$.

The convergence of Theorem 7 is quite general and it holds for any linearly convergent methods that can inexactly solve (34). However, in case that the iterative method is known we can have more concrete results. See below the more specified results for the cases of Conjugate gradient (CG) and Sketch and project method (SPM).

Convergence of InexactCG:

CG is deterministic iterative method for solving linear systems $\mathbf{A}x=b$ with symmetric and positive definite matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ in finite number of iterations. In particular, it can be shown that converges to the unique solution in at most $n$ steps. The worst case behavior of CG is given by [55, 17] 131313A sharper convergence rate of CG [55] for solving $\mathbf{A}x=b$ can be also used

$\|x_{k}-x_{*}\|^{2}_{\mathbf{A}}\leq\left(\frac{\lambda_{n-k}-\lambda_{1}}{\lambda_{n-k}+\lambda_{1}}\right)^{2}\|x_{0}-x_{*}\|^{2}_{\mathbf{A}},$

where matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ has $\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{n}$ eigenvalues.:

[TABLE]

where $x_{k}$ is the $k^{th}$ iteration of the method and $\kappa(\mathbf{A})$ the condition number of matrix $\mathbf{A}$.

Having present the convergence of CG for general linear systems, let us now return back to our setting. We denote $\lambda_{k}^{r}\in\mathbb{R}^{q}$ to be the approximate solution of the inner linear system (34) after $r$ conjugate gradient steps. Thus using (43) we know that $\|\lambda_{k}^{r}-\lambda_{k}^{*}\|_{\mathbf{M}_{k}}^{2}\leq\rho_{\mathbf{S}_{k}}^{4r}\|\lambda_{k}^{0}-\lambda_{k}^{*}\|^{2}_{\mathbf{M}_{k}},$ where $\rho_{\mathbf{S}_{k}}=\left(\frac{\sqrt{\kappa(\mathbf{M}_{k})}-1}{\sqrt{\kappa(\mathbf{M}_{k})}+1}\right)$. Now by making the same assumption as the general Theorem 7 the InexactCG converges with $\mathbb{E}[\|x_{k}-x_{*}\|_{{\bf B}}^{2}]\leq\left[1-\left(1-\theta_{CG}^{r}\right)\lambda_{\min}^{+}\right]^{k}\|x_{0}-x_{*}\|_{{\bf B}}^{2},$ where $\theta_{CG}<1$ such that $\rho_{\mathbf{S}_{k}}=\left(\frac{\sqrt{\kappa(\mathbf{M}_{k})}-1}{\sqrt{\kappa(\mathbf{M}_{k})}+1}\right)^{4}\leq\theta_{CG}$ with probability 1.

Convergence of InexactSP:

In this setting we suggest to run the sketch and project method (SPM) for solving inexactly the linear system (34). This allow us to have no assumptions on the structure of the original system $\mathbf{A}x=b$ and as a result we are able to solve more general problems compared to what problems InexactCG can solve141414Recall that InexactCG requires the matrix $\mathbf{M}_{k}$ to be positive definite (this is true when matrix $\mathbf{A}$ is a full rank matrix). Like before, by making the same assumptions as in Theorem 7 the more specific convergence $\mathbb{E}[\|x_{k}-x_{*}\|_{{\bf B}}^{2}]\leq\left[1-\left(1-\theta_{SP}^{r}\right)\lambda_{\min}^{+}\right]^{k}\|x_{0}-x_{*}\|_{{\bf B}}^{2},$ for the InexactSP can be obtained. Now the quantity $\rho_{\mathbf{S}_{k}}$ denotes the convergence rate of the exact Basic method151515Recall that iBasic and its exact variant ($\epsilon_{k}=0$) can be expressed as sketch and project methods (20). when this applied to solve linear system (34) and $\theta_{SP}<1$ is a scalar such that $\rho_{\mathbf{S}_{k}}\leq\theta_{SP}$ with probability 1.

5 Inexact Dual Method

In the previous sections we focused on the analysis of inexact stochastic methods for solving the stochastic optimization problem (1) and the best approximation (7). In this section we turn into the dual of the best approximation (10) and we propose and analyze an inexact variant of the SDSA (11). We call the new method iSDSA and is formalized as Algorithm 3. In the update rule $\epsilon^{d}_{k}$ indicates the dual inexactness error that appears in the $k^{th}$ iteration of iSDSA.

5.1 Correspondence Between the Primal and Dual Methods

With the sequence of the dual iterates $\{y_{k}\}_{k=0}^{\infty}$ produced by the iSDSA we can associate a sequence of primal iterates $\{x_{k}\}_{k=0}^{\infty}$ using the affine mapping (13). In our first result we show that the random iterates produced by iBasic arise as an affine image of iSDSA under this affine mapping.

Theorem 8**.**

*(Correspondence between the primal and dual methods)

Let $\{x_{k}\}_{k=0}^{\infty}$ be the iterates produced by iBasic (Algorithm 1). Let $y_{0}=0$, and $\{y_{k}\}_{k=0}^{\infty}$ the iterates of the iSDSA. Assume that the two methods use the same stepsize $\omega>0$ and the same sequence of random matrices $\mathbf{S}_{k}$. Assume also that $\epsilon_{k}=\mathbf{B}^{-1}\mathbf{A}^{\top}\epsilon^{d}_{k}$ where $\epsilon_{k}$ and $\epsilon^{d}_{k}$ are the inexactness errors appear in the update rules of iBasic and iSDSA respectively. Then*

[TABLE]

for all $k\geq 0$. That is, the primal iterates arise as affine images of the dual iterates.

Proof.

[TABLE]

Thus by choosing the inexactness error of the primal method to be $\epsilon_{k}=\mathbf{B}^{-1}\mathbf{A}^{\top}\epsilon^{d}_{k}$ the sequence of vectors $\{\phi(y_{k})\}$ satisfies the same recursion as the sequence $\{x_{k}\}$ defined by iBasic. It remains to check that the first element of both recursions coincide. Indeed, since $y_{0}=0$, we have $x_{0}=\phi(0)=\phi(y_{0})$. ∎

5.2 iSDSA with Structured Inexactness Error

In this subsection we present Algorithm 4. It can be seen as a special case of iSDSA but with a more structured inexactness error.

Similar to their primal variants, it can be easily checked that Algorithm 4 is a special case of the iSDSA ( Algorithm 3) when the dual inexactness error is chosen to be $\epsilon_{k}^{d}=\mathbf{S}_{k}(\lambda_{k}^{r}-\lambda_{k}^{*})$. Note that, using the observation of Remark 2 that $\epsilon_{k}=\omega\mathbf{B}^{-1}\mathbf{A}^{\top}\mathbf{S}_{k}(\lambda_{k}^{r}-\lambda_{k}^{*})$ and the above expression of $\epsilon_{k}^{d}$ we can easily verify that the expression $\epsilon_{k}=\mathbf{B}^{-1}\mathbf{A}^{\top}\epsilon^{d}_{k}$ holds. This is precisely the connection between the primal and dual inexactness errors that have already been used in the proof of Theorem 8.

5.3 Convergence of Dual Function Values

We are now ready to state a linear convergence result describing the behavior of the inexact dual method in terms of the function values $D(y_{k})$. The following result is focused on the convergence of iSDSA by making similar assumption to Assumption Assumption 1${b}$. Similar convergence results can be obtained using any other assumption of Section 3.1. The convergence of Algorithm 4, can be also easily derived using similar arguments with the one presented in Section 4 and the convergence guarantees of Theorem 7.

Theorem 9**.**

(Convergence of dual objective). Assume exactness. Let $y_{0}=0$ and let $\{y_{k}\}_{k=0}^{\infty}$ to be the dual iterates of iSDSA (Algorithm 3) with $\omega\in(0,2)$. Set $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$ and let $y_{*}$ be any dual optimal solution. Consider the inexactness error $\epsilon_{k}^{d}$ be such that it satisfies $\mathbb{E}[\|\mathbf{B}^{-1}\mathbf{A}^{\top}\epsilon_{k}^{d}\|^{2}_{\mathbf{B}}\;|\;y_{k},\mathbf{S}_{k}]\leq\sigma_{k}^{2}=q^{2}2\left[D(y_{*})-D(y_{k})\right]$ where $0\leq q<1-\sqrt{\rho}$. Then

[TABLE]

Proof.

The proof follows by applying Theorem 2 together with Theorem 8 and the identity $\tfrac{1}{2}\|x_{k}-x_{*}\|^{2}_{\bf B}=D(y_{*})-D(y_{k})$ (14). ∎

Note that in the case that $q=0$, iSDSA simplifies to its exact variant SDSA and the convergence rate coincide with the one presented in [29, 19]. Following similar arguments to those in [19], the same rate can be proved for the duality gap $\mathbb{E}[P(x_{k})-D(y_{k})]$.

6 Numerical Evaluation

In this section we perform preliminary numerical tests for studying the computational behavior of iBasic with structured inexactness error when is used to solve the best approximation problem (7) or equivalently the stochastic optimization problem (1)161616Note that from Section 5 and the correspondence between the primal and dual methods, iSDSA will have similar behavior when is applied to the dual problem (10).. As we have already mentioned, iBasic can be interpreted as sketch-and-project method, and as a result a comprehensive array of well-known algorithms can be recovered as special cases by varying the main parameters of the methods (Section 2.3). In particular, in our experiments we focus on the evaluation of two popular special cases, the inexact Randomized Block Kaczmarz (iRBK) (equation (21)) and inexact randomized block coordinate descent method (iRBCD) (equation (22))We implement Algorithm 2 presented in Section 4 using CG 171717Recall that in order to use CG, the matrix $\mathbf{M}_{k}$ that appears in linear system (34) should be positive definite. This is true in the case that the matrix $\mathbf{A}$ of the original system has full column rank matrix. Note however that the analysis of Section 4 holds for any consistent linear system $\mathbf{A}x=b$ and without making any further assumption on its structure or the linearly convergence methods. to inexactly solve the linear system of the update rule (equation (34)). Recall that in this case we named the method InexactCG.

The convergence analysis of previous sections is quite general and holds for several combinations of the two main parameters of the method, the positive definite matrix $\mathbf{B}$ and the distribution ${\cal D}$ of the random matrices $\mathbf{S}$. For obtaining iRBK as special case we have to choose $\mathbf{B}=\mathbf{I}\in\mathbb{R}^{n\times n}$ (Identity matrix) and for the iRBCD the given matrix $\mathbf{A}$ should be positive definite and choose $\mathbf{B}=\mathbf{A}$. For both methods the distribution ${\cal D}$ should be over random matrices $\mathbf{S}=\mathbf{I}_{:C}$ where $\mathbf{I}_{:C}$ is the column concatenation of the $m\times m$ identity matrix indexed by a random subset $C$ of $[m]$. In our experiments we choose to have one specific distribution over these matrices. In particular, we assume that the random matrix in each iteration is chosen uniformly at random to be $\mathbf{S}=\mathbf{I}_{:d}$ with the subset $d$ of $[m]$ to have fixed pre-specified cardinality.

The code for all experiments is written in the Julia 0.6.3 programming language and run on a Mac laptop computer (OS X El Capitan), 2.7 GHz Intel Core i5 with 8 GB of RAM.

To coincide with the theoretical convergence results of Algorithm 2 the relaxation parameter (stepsize) of the methods study in our experiments is chosen to be $\omega=1$ (no relaxation). In all implementations, we use $x_{0}=0\in\mathbb{R}^{n}$ as an initial point and in comparing the methods with their inexact variants we use the relative error measure $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}/\|x_{0}-x_{*}\|^{2}_{\mathbf{B}}\overset{x_{0}=0}{=}\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}/\|x_{*}\|^{2}_{\mathbf{B}}$. We run each method (exact and inexact) until the relative error is below $10^{-5}$. For the horizontal axis we use either the number of iterations or the wall-clock time measured using the tic-toc Julia function. In the exact variants, the linear system (34) in Algorithm 2 needs to be solved exactly. In our experiments we follow the implementation of [18] for both exact RBCD and exact RBK where the built-in direct solver (sometimes referred to as ”backslash”) is used.

Experimental setup:

For the construction of consistent linear systems $\mathbf{A}x=b$ we use the following setup:

• •

For iRBK: Let matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ being given (it can be either synthetic or real data). Then a vector $z\in\mathbb{R}^{n}$ is chosen to be i.i.d $\mathcal{N}(0,1)$ and the right hand side of the linear system is set to $b=\mathbf{A}z$. With this way the consistency of the linear system with matrix $\mathbf{A}$ and right hand side $b$ is ensured.

• •

For iRBCD: A Gaussian matrix $\mathbf{P}\in\mathbb{R}^{m\times n}$ is generated and then matrix $\mathbf{A}=\mathbf{P}^{\top}\mathbf{P}\in\mathbb{R}^{n\times n}$ is used in the linear system (with this way matrix $\mathbf{A}$ is positive definite with probability 1). The vector $z\in\mathbb{R}^{n}$ is chosen to be i.i.d $\mathcal{N}(0,1)$ and again to ensure consistency of the linear system, the right hand side is set to $b=\mathbf{A}z$.

6.1 Importance of Large Block Size

Many recent works have shown that using larger block sizes can be very beneficial for the performance of randomized iterative algorithms [18, 44, 34, 27]. In Figure 2 we numerically verify this statement. We show that both RBK and RBCD (no inexact updates) outperform in number of iterations and wall clock time their serial variants where only one coordinate is chosen (block of size $d=1$) per iteration. This justify the necessity of choosing methods with large block sizes. Recall that this is precisely the class of algorithms that could have an expensive subproblem in their update rule which is required to be solved exactly and as a result can benefit the most from the introduction of inexactness.

6.2 Inexactness and Block Size (iRBCD)

In this experiment, we first construct a positive definite linear system following the previously described procedure for iRBCD. We first generate a Gaussian matrix $\mathbf{P}\in\mathbb{R}^{10000\times 7000}$ and then the positive definite matrix $\mathbf{A}=\mathbf{P}^{\top}\mathbf{P}\in\mathbb{R}^{7000\times 7000}$ is used to define a consistent liner system. We run iRBCD in this specific linear system and compare its performance with its exact variance for several block sizes $d$ (numbers of column of matrix $\mathbf{S}$). For evaluating the inexact solution of the linear system in the update rule we run CG for either 2, 5 or 10 iterations. In Figure 3, we plot the evolution of the relative error in terms of both the number of iterations and the wall-clock time.

We observe that for any block size the inexact methods are always faster in terms of wall clock time than their exact variants even if they require (as is expected) equal or larger number of iterations. Moreover it is obvious that the performance of the inexact method becomes much better than the exact variant as the size $d$ increases and as a results the sub-problem that needs to be solved in each step becomes more expensive. It is worth to highlight that for the chosen systems, the exact RBCD behaves better in terms of wall clock time as the size of block increases (this coincides with the findings of the previous experiment).

6.3 Evaluation of iRBK

In the last experiment we evaluate the performance of iRBK in both synthetic and real datasets. For computing the inexact solution of the linear system in the update rule we run CG for pre-specified number of iterations that can vary depending the datasets. In particular, we compare iRBK and RBK on synthetic linear systems generated with the Julia Gaussian matrix functions “randn(m,n)” and “sprandn(m,n,r)” (input $r$ of sprandn function indicates the density of the matrix). For the real datasets, we test the performance of iRBK and RBK using real matrices from the library of support vector machine problems LIBSVM [8]. Each dataset of the LIBSVM consists of a matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ ($m$ features and $n$ characteristics) and a vector of labels $b\in\mathbb{R}^{m}$. In our experiments we choose to use only the matrices of the datasets and ignore the label vectors 181818Note that the real matrices of the Splice and Madelon datasets are full rank matrices.. As before, to ensure consistency of the linear system, we choose a Gaussian vector $z\in\mathbb{R}^{n}$ and the right hand side of the linear system is set to $b={\bf A}z$ (for both the synthetic and the real matrices). By observing Figure 4 it is clear that for all problems under study the performance of iRBK in terms of wall clock time is much better than its exact variant RBK.

7 Conclusion

In this work we propose and analyze inexact variants of several stochastic algorithms for solving quadratic optimization problems and linear systems. We provide linear convergence rate under several assumptions on the inexactness error. The proposed methods require more iterations than their exact variants to achieve the same accuracy. However, as we show through our numerical evaluations, the inexact algorithms require significantly less time to converge.

With the continuously increasing size of datasets, inexactness should definitely be a tool that practitioners should use in their implementations even in the case of stochastic methods that have much cheaper-to-compute iteration complexity than their deterministic variants. Recently, accelerated and parallel stochastic optimization methods [29, 46, 54] have been proposed for solving linear systems. We speculate that the addition of inexactness to these update rules will lead to methods faster in practice. We also believe that our approach and complexity results can be extended to the more general case of minimization of convex and non-convex functions in the stochastic setting. Finally, sketch-and-project algorithms have been used for solving the average consensus problem [27, 20] popular in distributed optimization literature. Our results could also be useful in this area and lead to the development of novel randomized gossip algorithms that use inexactness in their update rule.

8 Acknowledgements

The first author would like to acknowledge Robert Mansel Gower, Georgios Loizou and Rachael Tappenden for useful discussions.

Appendix A Technical Preliminaries

Lemma 10**.**

*(Lemma 4.2 [46]: Quadratic bounds)

For all $x\in\mathbb{R}^{n}$ and $x_{*}\in{\cal L}$ the following hold: $\lambda_{\rm min}^{+}f(x)\leq\frac{1}{2}\|\nabla f(x)\|_{\mathbf{B}}^{2}\leq\lambda_{\rm max}f(x)$ and $f(x)\leq\frac{\lambda_{\rm max}}{2}\|x-x_{*}\|^{2}_{\mathbf{B}}$. Furthermore, if exactness is satisfied and we let $x_{*}=\Pi_{{\cal L},{\bf B}}(x_{0})$ we have*

[TABLE]

Lemma 11**.**

[46*]**

Let $x_{*}\in{\cal L}$ and $\{x_{k}\}_{k\geq 0}$ be the random iterates produced by the exact Basic method (Algorithm 1 with $\epsilon_{k}=0$) with an arbitrary stepsize $\omega\in\mathbb{R}$. Then:*

[TABLE]

By taking expectation condition on $x_{k}$ (that is, the expectation is with respect to $\mathbf{S}_{k}$) and assuming $\omega\in(0,2)$ we can further obtain:

[TABLE]

Remark 7**.**

Let $x$ and $y$ be random vectors and let $\sigma$ positive constant. If we assume $\mathbb{E}[\|x\|_{{\bf B}}^{2}\;|\;y]\leq\sigma^{2}$ then by using the variance inequality (check Table 3) we obtain $\mathbb{E}[\|x\|_{{\bf B}}\;|\;y]\leq\sigma$. In our setting if we assume $\mathbb{E}[\|\epsilon_{k}\|_{{\bf B}}^{2}\;|\;x_{k},\mathbf{S}_{k}]\leq\sigma_{k}^{2}$ where $\epsilon_{k}$ is the inexactness error and $x_{k}$ is the current iterate then by the variance inequality it holds that $\mathbb{E}[\|\epsilon_{k}\|_{{\bf B}}\;|\;x_{k},\mathbf{S}_{k}]\leq\sigma_{k}$.

Appendix B Proofs of Main Results

In our convergence analysis we use several popular inequalities. Look Table 3 in Appendix C for the abbreviations and the relevant formulas.

A key step in the proofs of the theorems is to use the tower property of the expectation. We use it in the form

[TABLE]

where $X$ is some random variable. In all proofs we perform the three expectations in order, from the innermost to the outermost. Similar to the main part of the paper we use $\rho=1-\omega(2-\omega)\lambda_{\rm min}^{+}$.

B.1 Proof of Theorem 1

Proof.

First we decompose:

[TABLE]

Applying the innermost expectation of (49) to (50), we get:

[TABLE]

We now analyze the three expression T1,T2,T3 separately.

Note that an upper bound for the expression T2 can be directly obtained from the assumption

[TABLE]

The first expression can be written as:

[TABLE]

For expression T3:

[TABLE]

By substituting the bounds (52), (53), and (54) into (51) we obtain:

[TABLE]

We now take the middle expectation (see (49)) and apply it to inequality (55):

[TABLE]

Now let us find a bound on the quantity $\mathbb{E}\left[\|({\bf I}-\omega{\bf B}^{-1}{\bf Z}_{k})(x_{k}-x_{*})\|_{\mathbf{B}}\;|\;x_{k}\right]$. Note that from (48) and (47) we have that $\mathbb{E}\left[\|({\bf I}-\omega{\bf B}^{-1}{\bf Z}_{k})(x_{k}-x_{*})\|^{2}_{\mathbf{B}}\;|\;x_{k}\right]\leq\rho\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$. By using Remark 7 in the last inequality we obtain:

[TABLE]

By substituting (57) in (56):

[TABLE]

We take the final expectation (outermost expectation in the tower rule (49)) on the above expression to find:

[TABLE]

Using $r_{k}={\mathbb{E}\left[\|x_{k}-x_{*}\|_{{\bf B}}^{2}\right]}$ equation (59) takes the form:

[TABLE]

If we further substitute $p_{k}=\sqrt{r_{k}}$ and $\ell=\sqrt{\rho}$ the recurrence simplifies to:

[TABLE]

By unrolling the final inequality:

[TABLE]

Hence,

[TABLE]

The result is obtained by using V.I in the last expression. ∎

B.2 Proof of Corollary 1

By denoting $r_{k}=\mathbb{E}[\|x_{k}-x_{*}\|_{{\bf B}}]$ in (28) we obtain:

[TABLE]

Since $1-\rho^{k}\leq 1$ the result is obtained.

B.3 Proof of Theorem 2

In order to prove Theorem 2 we need to follow similar steps to the proof of Theorem 1. The main differences of the two proofs appear at the points that we need to upper bound the norm of the inexactness error ($\|\epsilon_{k}\|^{2}$). In particular instead of using the general sequence $\sigma_{k}^{2}\in\mathbb{R}$ we utilize the bound $q^{2}\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}$ from Assumption Assumption 1${b}$. Thus, it is sufficient to focus at the parts of the proof that these bound is used.

Similar to the proof of Theorem 1 we first decompose to obtain the equation (51). There, the expression T1 can be upper bounded from (53) but now using the Assumption Assumption 1${b}$ the expression T2 and T3 can be upper bounded as follows:

[TABLE]

[TABLE]

As a result by substituting the bounds (53), (60), and (61) into (51) we obtain:

[TABLE]

By following the same steps to the proof of Theorem 1 the equation (58) takes the form:

[TABLE]

We take the final expectation (outermost expectation in the tower rule (49)) on the above expression to find:

[TABLE]

The final result follows by unrolling the recurrence.

B.4 Proof of Theorem 3

Proof.

Similar to the previous two proofs by decomposing the update rule and using the innermost expectation of (49) we obtain equation (51). An upper bound of expression T1 is again given by inequality (53). For the expression T2 depending the assumption that we have on the norm of the inexactness error different upper bounds can be used. In particular,

1. (i)

If Assumption Assumption 1 holds then: $T2=\mathbb{E}[\|\epsilon_{k}\|^{2}_{\mathbf{B}}\;|\;x_{k},\mathbf{S}_{k}]\leq\sigma_{k}^{2}.$ 2. (ii)

If Assumption Assumption 1${b}$ holds then: $T2=\mathbb{E}[\|\epsilon_{k}\|^{2}_{\mathbf{B}}\;|\;x_{k},\mathbf{S}_{k}]\leq\sigma_{k}^{2}=q^{2}\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}.$ 3. (iii)

If Assumption Assumption 1${c}$ holds then: $T2=\mathbb{E}[\|\epsilon_{k}\|^{2}_{\mathbf{B}}\;|\;x_{k},\mathbf{S}_{k}]\leq\sigma_{k}^{2}=2q^{2}f_{S_{k}}(x_{k}).$

The main difference from the previous proofs, is that due to the Assumption Assumption 2 and tower property (49) the expression T3 will eventually be equal to zero. More specifically, we have that:

[TABLE]

Thus, in this case equation (55) takes the form:

[TABLE]

Using the above expression depending the assumption that we have we obtain the following results:

1. (i)

By taking the middle expectation (see (49)) and apply it to the above inequality:

[TABLE]

We take the final expectation (outermost expectation in the tower rule (49)) on the above expression to find:

[TABLE]

Using $r_{k}={\mathbb{E}\left[\|x_{k}-x_{*}\|_{{\bf B}}^{2}\right]}$ the last inequality takes the form $r_{k+1}\leq\rho r_{k}+{\bar{\sigma}}_{k}^{2}$. By unrolling the last expression: $r_{k}\leq\rho^{k}r_{0}+(\rho^{0}{\bar{\sigma}}_{k-1}^{2}+\rho{\bar{\sigma}}^{2}_{k-2}+\cdots+\rho^{k-1}{\bar{\sigma}}^{2}_{0})=\rho^{k}r_{0}+\sum_{i=0}^{k-1}\rho^{k-1-i}{\bar{\sigma}}^{2}_{i}.$ Hence,

[TABLE] 2. (ii)

For the case (ii) inequality (65) takes the form:

[TABLE]

and by taking the middle expectation (see (49)) we obtain:

[TABLE]

By taking the final expectation of the tower rule (49) and apply it to the above inequality:

[TABLE]

and the result is obtain by unrolling the last expression. 3. (iii)

For the case (iii) inequality (65) takes the form:

[TABLE]

and by taking the middle expectation (see (49)) we obtain:

[TABLE]

By taking the final expectation of the tower rule (49) to the above inequality:

[TABLE]

and the result is obtain by unrolling the last expression.

∎

Appendix C Useful Inequalities and Frequently Used Notation