Convergence Analysis of Inexact Randomized Iterative Methods

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.

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

Stochastic Quadratic Optimization 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.

We study the stochastic quadratic optimization problem

first proposed in for reformulating consistent linear systems

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

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 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 , problem (1) was solved via Stochastic Gradient Descent (SGD)The gradient is computed with respect to the inner product $\langle{\bf B}x,y\rangle$.:

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 :

In all methods $\omega>0$ is a fixed stepsize and ${\bf S}_{k}$ is sampled afresh in each iteration from distribution ${\cal D}$. See 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 , 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:

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):

first proposed in . 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

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 .

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

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:

It can be proved, , 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:

In the dual method was analyzed for the case of unit stepsize ($\omega=1$). Later in 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 .

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 . That is,

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).

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 , 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.

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.

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.

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

In order to keep the expression brief throughout the paper we defineIn 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}$.:

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

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

which has the same spectrum with the matrix $\nabla^{2}f(x)$ but at the same time is symmetric and positive semi-definiteNote 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 that $0\leq\lambda_{i}\leq 1$ for all $i\in[n]$.

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.

Let us devote this subsection on explaining how the inexactness can be introduced in the current exact update rules of SGDNote 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:

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 The Inexact Basic Method

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:

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.

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 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.

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

Inexact Randomized Block Coordinate Descent (iRBCD)In 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 ).: 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

For more papers related to Kaczmarz method (randomized, greedy, cyclic update rules) we refer the interested reader to . For the coordinate descent method (a.k.a Gauss-Seidel for linear systems) and its block variant, Randomized Block Coordinate Descent we suggest .

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 . 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 for Inexact Newton methods, for Inexact Proximal Point methods and 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 . 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 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 .

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 .

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 ).

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.

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.

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}}$.

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.

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.

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

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,

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:

Inspired from , 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 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.

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

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$ .

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.

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:

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

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

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.

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:

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:

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,

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.

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:

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). ∎

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:

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:

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.

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}.$ ∎

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.

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}$.

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 In 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:

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:

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,

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.

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.

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.

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.

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).

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

Convergence of InexactSP:

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.

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.

(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

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

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})$. ∎

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.

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.

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). ∎

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)Note 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 Recall 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 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.

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

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 . 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.

2 Inexactness and Block Size (iRBCD)

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).

3 Evaluation of iRBK

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 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 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.

Acknowledgements

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

References

Appendix A Technical Preliminaries

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:

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

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}^{+}$.

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

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

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

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

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

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

B.2 Proof of Corollary 1

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

B.3 Proof of Theorem 2

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:

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

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

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

The final result follows by unrolling the recurrence.

B.4 Proof of Theorem 3

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,

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:

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

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

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

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

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

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

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

and the result is obtain by unrolling the last expression.

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

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

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

and the result is obtain by unrolling the last expression.

Appendix C Useful Inequalities and Frequently Used Notation