Stochastic Reformulations of Linear Systems

,


Introduction
Linear systems form the backbone of most numerical codes used in academia and industry.With the advent of the age of big data, practitioners are looking for ways to solve linear systems of unprecedented sizes.The present work is motivated by the need to design such algorithms.As an algorithmic tool enabling faster computation, randomization is well developed, understood and appreciated in several fields, typically traditionally of a "discrete" nature, most notably theoretical computer science [28].However, probabilistic ideas are also increasingly and successfully penetrating "continuous" fields, such as numerical linear algebra [57,10,11,60,27,50,1], optimization [24,33,35,46,59,40], control theory [4,5,62], machine learning [52,47,21,14,44], and signal processing [7,22].
In this work we are concerned with the problem of solving a consistent linear system.In particular, consider the problem solve Ax = b, where 0 = A ∈ R m×n .We shall assume throughout the paper that the system is consistent, i.e., L def = {x : Ax = b} = ∅.Problem ( 1) is arguably one of the most important problems in linear algebra.As such, a tremendous amount of research has been done to design efficient iterative algorithms [51].However, surprisingly little is know about randomized iterative algorithms for solving linear systems.In this work we aim to contribute to closing this gap.

Stochastic reformulations of linear systems
We propose a fundamental and flexible way of reformulating each consistent linear system into a stochastic problem.To the best of our knowledge, this is the first systematic study of such reformulations.Stochasticity is introduced in a controlled way, into an otherwise deterministic problem, as a decomposition tool which can be leveraged to design efficient, granular and scalable randomized algorithms.
Parameters defining the reformulation.Stochasticity enters our reformulations through a userdefined distribution D describing an ensemble of random matrices S ∈ R m×q .We make use of one more parameter: a user-defined n×n symmetric positive definite matrix B. Our approach and underlying theory support virtually all thinkable distributions.The choice of the distribution should ideally depend on the problem itself, as it will affect the conditioning of the reformulation.However, for now we leave such considerations aside.
One stochastic reformulation in four disguises.Our reformulation of (1) as a stochastic problem has several seemingly different, yet equivalent interpretations, and hence we describe them here side by side.
1. Stochastic optimization problem.Consider the problem minimize f (x) where f S (x) = 1 2 (Ax − b) H(Ax − b), H = S(S AB −1 A S) † S , and † denotes the Moore-Penrose pseudoinverse.When solving the problem, we do not have (or do not wish to exercise, as it may be prohibitively expensive) explicit access to f , its gradient or Hessian.Rather, we can repeatedly sample S ∼ D and receive unbiased samples of these quantities at points of interest.That is, we may obtain local information about the stochastic function f S , such as the stochastic gradient ∇f S (x), and use this to drive an iterative process for solving (2).
2. Stochastic linear system.Consider now a preconditioned version of the linear system (1) given by solve where P = B −1 A E S∼D [H] is the preconditioner.The preconditioner is not assumed to be known explicitly.Instead, when solving the problem, we are able to repeatedly sample S ∼ D, obtaining an unbiased estimate of the preconditioner (not necessarily explicitly), B −1 A H, for which we coin the name stochastic preconditioner.This gives us access to an unbiased sample of the preconditioned system (3): B −1 A HAx = B −1 A Hb.As we shall see-in an analogy with stochastic optimization-the information contained in such systems can be utilized by an iterative algorithm to solve (3).
That is, we seek to find a fixed point of the mapping x → E S∼D Π B L S (x) .When solving the problem, we do not have an explicit access to the average projection map.Instead, we are able to repeatedly sample S ∼ D, and use the stochastic projection map x → Π B L S (x). 4. Probabilistic intersection problem.Note that L ⊆ L S for all S. We would wish to design D in such a way that a suitably chosen notion of an intersection of the sets L S is equal to L. The correct notion is what we call probabilistic intersection, denoted ∩ S∼D L S , and defined as the set of points x which belong to L S with probability one.This leads to the problem: As before, we typically do not have an explicit access to the probabilistic intersection when designing an algorithm.Instead, we can repeatedly sample S ∼ D, and utilize the knowledge of L S to drive the iterative process.If D is a discrete distribution, probabilistic intersection reduces to standard intersection.
All of the above formulations have a common feature: they all involve an expectation over S ∼ D, and we either do not assume this expectation is known explicitly, or even if it is, we prefer, due to efficiency or other considerations, to sample from unbiased estimates of the objects (e.g., stochastic gradient ∇f S , stochastic preconditioner B −1 A H, stochastic projection map x → Π B L S (x), random set L S ) appearing in the formulation.
Equivalence and exactness.We show that all these stochastic reformulations are equivalent (see Theorem 3.4).In particular, the following sets are identical: the set of minimizers of the stochastic optimization problem (2), the solution set of the stochastic linear system (3), the set of fixed points of the stochastic fixed point problem (4), and the probabilistic intersection (5).Further, we give necessary and sufficient conditions for this set to be equal to L. If this is the case, we say the the reformulation is exact (see Section 3.4).Distributions D satisfying these conditions always exist, independently of any assumptions on the system beyond consistency.The simplest, but also the least useful choice of a distribution is to pick S = I (the m × m identity matrix), with probability one.In this case, all of our reformulations become trivial.

Stochastic algorithms
Besides proposing a family of stochastic reformulations of the linear system (1), we also propose three stochastic algorithms for solving them: Algorithm 1 (basic method), Algorithm 2 (parallel/minibatch method), and Algorithm 3 (accelerated method).Each method can be interpreted naturally from the viewpoint of each of the reformulations.
Basic method.Below we list some of the interpretations of Algorithm 1 (basic method), which performs updates of the form where S k ∼ D is sampled afresh in each iteration.The method is formally presented and analyzed in Section 4.
1. Stochastic gradient descent.Algorithm 1 can be seen as stochastic gradient descent [48], with fixed stepsize, applied to (2).At iteration k of the method, we sample S k ∼ D, and compute ∇f S k (x k ), which is an unbiased stochastic approximation of ∇f (x k ).We then perform the step where ω > 0 is a stepsize.
Let us note that in order to achieve linear convergence it is not necessary to use any explicit variance reduction strategy [52,21,23,9], nor do we need to use decreasing stepsizes.This is because the stochastic gradients vanish at the optimum, which is a consequence of the consistency assumption.Surprisingly, we get linear convergence in spite of the fact that we deal with a nonfinite-sum problem (2), and without the need to assume boundedness of the stochastic gradients, and without f being strongly convex.To the best of our knowledge, this is the first linearly convergent accelerated method for stochastic optimization without requiring strong convexity.This beats the minimax bounds given by Srebro [58].This is because (2) is not a black-box stochastic optimization objective; indeed, we have constructed it in a particular way from the underlying linear system (1).
2. Stochastic Newton method.However, Algorithm 1 can also be seen as a stochastic Newton method.At iteration k we sample S k ∼ D, and instead of applying the inverted Hessian of f S k to the stochastic gradient (this is not possible as the Hessian is not necessarily invertible), we apply the B-pseudoinverse.That is, we perform the step where ω > 0 is a stepsize, and the B-pseudoinverse of a matrix M is defined as One may wonder why methods (7) and (8) are equivalent; after all, the (stochastic) gradient descent and (stochastic) Newton methods are not equivalent in general.However, in our setting it turns out that the stochastic gradient ∇f S k (x) is always an eigenvector of (∇ 2 f S k (x)) † B , with eigenvalue 1 (see Lemma 3.1).
Stochastic Newton-type methods were recently developed and analyzed in the optimization and machine learning literature [39,38,41,29].However, they are design to solve different problems, and operate in a different manner.

Stochastic proximal point method.
If we restrict our attention to stepsizes satisfying 0 < ω ≤ 1, then Algorithm 1 can be equivalently (see Theorem A.3 in the Appendix) written down as That is, ( 9) is a stochastic variant of the proximal point method for solving (2), with a fixed regularization parameter [49].The proximal point method is obtained from (9) by replacing f S k with f .If we define the prox operator of a function ψ : R n → R with respect to the B-norm as prox B ψ (y) , then iteration ( 9) can be written compactly as 4. Stochastic fixed point method.From the perspective of the stochastic fixed point problem (4), Algorithm 1 can be interpreted as a stochastic fixed point method, with relaxation.We first reformulate the problem into an equivalent form using relaxation, which is done to improve the contraction properties of the map.We pick a parameter ω > 0, and instead consider the equivalent fixed point problem x = E S∼D ωΠ B L S (x) + (1 − ω)x .Now, at iteration k, we sample S k ∼ D, which enables us to obtain an unbiased estimate of the new fixed point mapping, and then simply perform one step of a fixed point method on this mapping: 5. Stochastic projection method.Algorithm 1 can also be seen as a stochastic projection method applied to the probabilistic intersection problem (5).By sampling S k ∼ D, we are one of the sets defining the intersection, namely L S k .We then project the last iterate onto this set, in the B-norm, followed by a relaxation step with relaxation parameter ω > 0. That is, we perform the update This is a randomized variant of an alternating projection method.Note that the representation of L as a probabilistic intersection of sets is not given to us.Rather, we construct it with the hope to obtain faster convergence.
An optimization algorithm utilizing stochastic projection steps was developed in [30].For a comprehensive survey of projection methods for convex feasibility problems, see [2].
Parallel method.A natural parallelization strategy is to perform one step of the basic method independently τ times, starting from the same point x k , and average the results: where S 1 k , . . ., S τ k are independent samples from D (recall that φ ω is defined in (6)).This method is formalized as Algorithm 2, and studied in Section 5.1.Betrayed by our choice of the name, this method is useful in scenarios where τ parallel workers are available, allowing for the τ basic steps to be computed in parallel, followed by an averaging operation.
From the stochastic optimization viewpoint, this is a minibatch method.Considering the SGD interpretation (7), we can equivalently write (12) in the form . This is minibatch SGD.Iteration complexity of minibatch SGD was first understood in the context of training support vector machines with the hinge loss [61].Complexity under a lock-free paradigm, in a different setting from ours, was first studied in [36].Notice that in the limit τ → ∞, we obtain gradient descent.It is therefore interesting to study the complexity of the parallel method as a function τ .Of course, this method can also be interpreted as a minibatch stochastic Newton method, minibatch proximal point method and so on.
From the probabilistic intersection point of view, method (12) can be interpreted as a stochastic variant of the parallel projection method.In particular, we obtain the iterative process That is, we move from the current iterate, x k , towards the average of the τ projection points, with undershooting (if ω < 1), precisely landing on (if ω = 1), or overshooting (if ω > 1) the average.Projection methods have a long history and are well studied [12,3].However, much less is known about stochastic projection methods.
Accelerated method.In order to obtain acceleration without parallelization-that is, acceleration in the sense of Nesterov [34]-we suggest to perform an update step in which x k+1 depends on both x k and x k−1 .In particular, we take two dependent steps of Algorithm 1, one from x k and one from x k−1 , and then take an affine combination of the results.That is, the process is started with x 0 , x 1 ∈ R n , and for k ≥ 1 involves an iteration of the form where the matrices {S k } are independent samples from D, and γ ∈ R is an acceleration parameter.
Note that by choosing γ = 1 (no acceleration), we recover Algorithm 1.This method is formalized as Algorithm 3 and analyzed in Section 5.2.Our theory suggests that γ should be always between 1 and 2.
In particular, for well conditioned problems (small ζ), one should choose γ ≈ 1, and for ill conditioned problems (large ζ), one should choose γ ≈ 2.

Complexity
The complexity of our methods is completely described by the spectrum of the (symmetric positive semidefinite) matrix Let W = UΛU be the eigenvalue decomposition of W, where U = [u 1 , . . ., u n ] are the eigenvectors, With all of the above reformulations we associate the same condition number Alg.
Table 1: Summary of the main complexity results.In all cases, x * = Π B L (x 0 ) (the projection of the starting point onto the solution space of the linear system)."Complexity" refers to the number of iterations needed to drive "Quantity" below some error tolerance > 0 (we suppress a log(1/ ) factor in all expressions in the "Complexity" column).In the table we use the following expressions: where • is the spectral norm, λ max is the largest eigenvalue of W and λ + min is the smallest nonzero eigenvalue of W. Note that, for example, ζ is the condition number of the Hessian of f , and also the condition number of the stochastic linear system (3).Natural interpretations from the viewpoint of the stochastic fixed point and probabilistic intersection problems are also possible.As one varies the parameters defining the reformulation (D and B), the condition number changes.For instance, choosing S = I with probability one gives ζ = 1.
Exact formula for the evolution of expected iterates.We first show (Theorem 4.3) that after the canonical linear transformation x → U B 1/2 x, the expected iterates of the basic method satisfy the identity where x * is an arbitrary solution of the linear system (i.e., x * ∈ L).This identity seems to suggest that zero eigenvalues cause an issue, preventing convergence of the corresponding elements of the error to zero.Indeed, if λ i = 0, then (15) implies that B for three different stepsize choices is summarized in the first three rows of Table 1.In particular, the default stepsize ω = 1 leads to the complexity 1/λ + min , the long stepsize ω = 1/λ max gives the improved complexity λ max /λ + min , and the optimal stepsize ω = 2/(λ max + λ + min ) gives the best complexity 0.5 + 0.5λ max /λ + min .However, if we are interested in the convergence of the larger quantity E x k − x * 2 B (L2 convergence), it turns out that ω = 1 is the optimal choice, leading to the complexity 1/λ + min .
Parallel and accelerated methods.The parallel method improves upon the basic method in that it is capable of faster L2 convergence.We give a complexity formula as a function of τ , recovering the complexity the 1/λ + min rate of the basic method in the τ = 1 case, and achieving the improved asymptotic complexity λ max /λ + min as τ → ∞ (recall that λ max ≤ 1, whence the improvement).Because of this, λ max is the quantity driving parallelizability.If λ max is close to 1, then there is little or no reason to parallelize.If λ max is very small, parallelization helps.The smaller λ max is, the more gain is achieved by utilizing more processors.
With the correct stepsize (ω) and acceleration (γ) parameters, the accelerated method improves the complexity λ max /λ + min achieved by the basic method to λ max /λ + min .However, this is established for the quantity We conjecture that the L2 convergence rate of the accelerated method (for a suitable choice of the parameters ω and γ) is 1/λ + min .

Stochastic preconditioning
We coin the phrase stochastic preconditioning to refer to the general problem of designing matrix B and distribution D such that the appropriate condition number of W is well behaved.For instance, one might be interested in minimizing (or reducing) the condition number 1/λ + min if the basic method with unit stepsize is used, and the quantity we wish to converge to zero is either (see Lines 1, 4 and 5 of Table 1).On the other hand, if we can estimate λ max , then one may use the basic method with the larger stepsize 1/λ max , in which case we may wish to minimize (or reduce) the condition number λ max /λ + min (see Line 2 of Table 1).One possible approach to stochastic preconditioning is to choose some B and then focus on a reasonably simple parametric family of distributions D, trying to find the parameters which minimize (or reduce) the condition number of interest.The distributions in this family should be "comparable" in the sense that the cost of one iteration of the method of interest should be comparable for all distributions; as otherwise comparing bounds on the number of iterations does not make sense.
To illustrate this through a simple example, consider the family of discrete uniform distributions over m vectors in R m (that is, S 1 , . . ., S m ∈ R m×1 ), where the vectors themselves are the parameters defining the family.The cost of one iteration of the basic method will be proportional to the cost of performing a matrix-vector product of the form S A, which is comparable across all distributions in the family (assuming the vectors are dense, say).To illustrate this approach, consider this family, and further assume that A is n × n symmetric and positive definite.Choose B = A. It can be shown that 1/λ + min is maximized precisely when {S i } correspond to the eigenvectors of A. In this case, 1/λ + min = n, and hence our stochastic preconditioning strategy results in a condition number which is independent of the condition number of A. If we now apply the basic method to the stochastic optimization reformulation, we can interpret it as a spectral variant of stochastic gradient descent (spectral SGD).Ignoring logarithmic terms, spectral SGD only needs to perform n matrix vector multiplications to solve the problem.While this is not a practical preconditioning strategy-computing the eigenvectors is hard, and if we actually had access to them, we could construct the solution directly, without the need to resort to an iterative scheme-it sheds light on the opportunities and challenges associated with stochastic preconditioning.
All standard sketching matrices S can be employed within our framework, including the count sketch [6] and the count-min sketch [8].In the context to this paper (since we sketch with the transpose of S), S is a count-sketch matrix (resp.count-min sketch) if it is assembled from random columns of [I, −I] (resp I), chosen uniformly with replacement, where I is the m × m identity matrix.
The notion of importance sampling developed in the last 5 years in the randomized optimization and machine learning literature [45,64,42,40] can be seen a type of stochastic preconditioning, somewhat reverse to what we have outlined above.In these methods, the atoms forming the distribution D are fixed, and one is seeking to associate them with appropriate probabilities.Thus, the probability simplex is the parameter space defining the class of distributions one is considering.
Stochastic preconditioning is fundamentally different from the idea of randomized preconditioning [50,1], which is based on a two-stage procedure.In the first step, the input matrix is randomly projected and an good preconditioning matrix is extracted.In the second step, an iterative least squares solver is applied to solve the preconditioned system.
Much like standard preconditioning, different stochastic preconditioning strategies will need to be developed for different classes of problems, with structure of A informing the choice of B and D. Due to its inherent difficulty, stochastic preconditioning is beyond the scope of this paper.

Notation
For convenience, a table of the most frequently used notation is included in Appendix D. All matrices are written in bold capital letters.By Range (M) and Null (M) we mean the range space and null space of matrix M, respectively.Given a symmetric positive definite matrix B ∈ R n×n , we equip R n with the Euclidean inner product defined by x, h B def = x Bh.We also define the induced norm: x, x B .The short-hand notation • means • I , where I is the identity matrix.We shall often write x M for matrix M ∈ R n×n being merely positive definite; this constitutes a pseudonorm.

Further Connections to Existing Work
In this section we outline several connections of our work with existing developments.We do not aim to be comprehensive.

Randomized Kaczmarz method, with relaxation and acceleration
Let B = I, and choose D as follows: S = e i with probability p i = A i: the condition number is Basic method.In this setup, Algorithm 1 simplifies to For ω = 1, this reduces to the celebrated randomized Kaczmarz method (RK) of Strohmer and Vershyinin [60].For ω > 1, this is RK with overrelaxation -a new method not considered before.Based on Theorem 4.6, for ω ∈ [1/λ max , ω * ] the iteration complexity of Algorithm 1 is . This is an improvement on standard RK method (with ω = 1), whose complexity depends on Trace A A instead of λ max .Thus, the improvement can be as large as by a factor n.
Accelerated method.In the same setup, Algorithm 3 simplifies to This is accelerated RK method with overrelaxation -a new method not considered before.Based on Theorem 5.3, for the parameter choice ω = 1/λ max and γ = 2/(1 + ζ −2 ), the iteration complexity of this method is .
If we instead choose ω = 1 and γ = 2/(1 + ζ −2 ), the iteration complexity gets slightly worse: To the best of our knowledge, this is the best known complexity for a variant of RK.Let us remark that an asynchronous accelerated RK method was developed in [25].
The randomized Kaczmarz method, its variants have received considerable attention recently [31,65,32,43], and several connections to existing methods were made.Kaczmarz-type methods in a Hilbert setting were developed in [37].

Basic method with unit stepsize
The method x k+1 ← Π B L S k (x k ) was first proposed and analyzed (under a full rank assumption on A) in [17].Note that in view of (10), this is the basic method with unit stepsize.However, it was not interpreted as a method for solving any of the reformulations presented here, and as a result, all the interpretations we are giving here also remained undiscovered.Instead, it was developed and presented as a method for finding the unique solution of (1).

Duality
As we have seen, all three methods developed in this paper converge to a specific solution of the linear system (1), namely, to the projection of the starting point onto the solution space: x * = Π B L (x 0 ).Therefore, our methods solve the best approximation problem The "dual" of the basic method with unit stepsize was studied in this context in [18].The Fenchel dual of the best approximation problem is an unconstrained concave quadratic maximization problem of the form max y∈R m D(y), where the dual objective D depends on A, b, B and x 0 .In [18] it was shown that the basic method with unit stepsize closely related to a dual method (stochastic dual subspace ascent) performing iterations of the form where S k ∼ D, and λ k is chosen in such a way that the dual objective is as large as possible.Notice that the dual method in each iteration performs exact "line search" in a random subspace of R m spanned by the columns of the random matrix S k , and passing through y k .In particular, the iterates of the basic method with unit stepsize arise as affine images of the iterates of the dual method: In a similar fashion, it is possible to interpret the methods developed in this paper as images of appropriately designed dual methods.In [18], the authors focus on establishing convergence of various quantities, such as dual objective, primal-dual gap, primal objective and so on.They obtain the complexity 1/λ + min , which is identical to the rate we obtain here for the basic method with unit stepsize.However, their results require a stronger assumption on D (their assumption implies exactness, but not vice versa).We perform a much deeper analysis from the novel viewpoint of stochastic reformulations of linear systems, include a stepsize, and propose and analyze parallel and accelerated variants.
In the special case when S is chosen to be a random unit coordinate vector, (18) specializes to the randomized coordinate descent method, first analyzed by Leventhal and Lewis [24].In the special case when S is chosen as a random column submatrix of the m × m identity matrix, (18) specializes to the randomized Newton method of Qu et al. [41].Randomized coordinate descent methods are the state of the art methods for certain classes of convex optimization problems with a very large number of variables.The first complexity analysis beyond quadratics was performed in [54,35,46], a parallel method was developed in [47], duality was explored in [55] and acceleration in [14].

Randomized gossip algorithms
It was shown in [17,26] that for a suitable matrix A encoding the structure of a graph, and for b = 0, the application of the randomized Kaczmarz and randomized block Kaczmarz methods to (17) lead to classical and new randomized gossip algorithms developed and studied in the signal processing literature, with new insights and novel proofs.Likewise, when applied in the same context, our new methods lead to new parallel and accelerated gossip algorithms.

Empirical risk minimization
Regularized empirical risk minimization (ERM) problems are optimization problems of the form where f i is a loss function and g a regularizer.Problems of this form are of key importance in machine learning and statistics [53].Let In this setting, the ERM problem (19) becomes equivalent to (17).While quadratic regularizers similar to g are common in machine learning, zero/infinity loss functions are not used.For this reason, this specific instance of ERM was not studied in the machine learning literature.Since all our methods solve (17), they can be seen as stochastic algorithms for solving the ERM problem (19).
Since there is no reason to expect that any of our methods will satisfy Ax k = b for any finite k, the ERM objective value can remain to be equal to +∞ throughout the entire iterative process.From this perspective, the value of the ERM objective is unsuitable as a measure of progress.

Matrix inversion and quasi-Newton methods
Given an invertible matrix A, its inverse is the unique solution of the matrix equation AX = I.In [20] the authors have extended the "sketch and project" method [17] to this equation.In each iteration of the method, one projects the current iterate matrix X k , with respect to a weighted Frobenius norm, onto the sketched equation S k AX = S k I.This is a similar iterative process to the basic method with unit stepsize.The authors of [20] prove that the iterates of method converge to the inverse matrix at a linear rate, and detail connections of their method to quasi-Newton updates and approximate inverse preconditioning.A limited memory variant of the stochastic block BFGS method has been used to develop new efficient stochastoc quasi-Newton methods for emprical risk minimization problems appearing in machine learning [16].
It is possible to approach the problem AX = I in the same way we approach the system (1) in our paper, writing down stochastic reformulations, and then developing new variants of the sketch and project method [20]: a basic method with a stepsize, and parallel and accelerated methods.This would lead to the development of new variants of stochastic quasi-Newton rules, notably parallel and accelerated block BFGS.We conjecture that these rules will have superior performance to classical BFGS in practice.
Similar extensions and improvements can be done in relation to the problem of computing the pseudoinverse of very large rectangular matrices [19].

Stochastic Reformulations of Linear Systems
In this section we formally derive the four stochastic formulations outlined in the introduction: stochastic optimization, stochastic linear system, stochastic fixed point problem and probabilistic intersection.Along the way we collect a number of results and observations which will be useful in the complexity analysis of our methods.

Projections
The projection onto L = {x : Ax = b} is given by Note that for B = I, we get A † I = A (AA ) † = A † , and hence the I-pseudoinverse reduces to the standard Moore-Penrose pseudoinverse.The B-pseudoinverse satisfies

Stochastic functions
Let D be an arbitrary distribution over m × q matrices, which we shall denote as S. We shall write S ∼ D to say that S is drawn from D. We shall often refer to matrix expressions involving S, A and B. In order to keep the expressions brief throughout the paper, it will be useful to define and = A S(S AB −1 A S) † S A.
Notice that B −1 Z is the projection, in the B-norm, onto Range B −1 A S .In particular, Given S ∼ D, we define the stochastic (random) function By combining ( 25) and ( 23), this can be also written in the form For each x, h ∈ R n we have the expansion are the gradient and Hessian of f S with respect to the B-inner product, respectively. 1In view of ( 23) and ( 27), the gradient can also be written as Identities (29) in the following lemma explain why algorithm (6) can be equivalently written as stochastic gradient descent (7), stochastic Newton method (8), stochastic fixed point method (10), and stochastic projection method (11).For instance, the identity (∇ 2 f S )∇f S (x) = ∇f S (x) means that the stochastic gradients of f S are eigenvectors of the stochastic Hessian ∇ 2 f S , corresponding to eigenvalue one.Lemma 3.1.For all x ∈ R n , we have Moreover, If L S is the set of minimizers of f S , then L ⊆ L S , and Finally, for all x ∈ R n we have the identity Proof.Pick any x * ∈ L. First, we have = x − ∇f S (x).To establish (29), it now only remains to consider the two expressions involving the Hessian.We have = ∇f S (x), and = ∇f S (x).
If x ∈ L, then by picking x * = x in (28), we see that x ∈ L S .It remains to show that the sets defined in (i)-(iv) are identical.Equivalence between (i) and (ii) follows from (28).Now consider (ii) and (iii).Any x * ∈ L belongs to the set defined in (iii), which follows immediately by substituting b = Ax * .The rest follows after observing the nullspaces are identical.In order to show that (iii) and (iv) are equivalent, it suffices to compute Π B L S (x) and observe that Π B L S (x) = x if and only if x belongs to the set defined in (iii).
It remains to establish (31).In view of (i), it suffices to show that x − ∇f S (x) ∈ L S .However, from (29) we know that x − ∇f S (x) = Π B L S (x) ∈ L S .

Stochastic reformulation
In order to proceed, we shall enforce a basic assumption on D. This is an assumption on D since a suitable distribution satisfying it exists for all A ∈ R m×n and B 0. Note that if the assumption holds, then E [H] is symmetric and positive semidefinite.We shall enforce this assumption throughout the paper and hence will not henceforth refer to it.
Example 1.Let D be the uniform distribution over unit basis vectors in R m .That is, S = e i (the ith unit basis vector in R m ) with probability 1/m.Then where α i = 1/ A i:

2
B −1 for i = 1, 2, . . ., m.If A has nonzero rows, then E [H] 0. In this paper we reformulate the linear system (1) as the stochastic optimization problem Under Assumption 3.2, the expectation in ( 32) is finite for all x, and hence f is well defined.The following is a direct consequence of Lemma 3.1.We shall use these formulas throughout the paper.Lemma 3.3 (Representations of f ).Function f defined in (32) can be represented in multiple ways: and for any x * ∈ L we can write Since E [H] 0, f is a convex quadratic function.Moreover, f is nonnegative.The gradient and Hessian of f (with respect to the B-inner product) are given by and respectively, where x * is any point in L.
The set of minimizers of f , denoted X , can be represented in several ways, as captured by our next result.It immediately follows that the four stochastic formulations mentioned in the introduction are equivalent.
Proof.As f is convex, nonnegative and achieving the value of zero (since L = ∅), the sets in (i) are all identical.We shall now show that the sets defined in (ii)-(iv) are equal to that defined in (i).Using the formula for the gradient from (35), we see that {x : ∇f (x) = 0} = {x : ), which shows that (i) and (ii) are the same.Equivalence of (i) and (iii) follows by taking expectations in (29) to obtain = E x − Π B L S (x) .It remains to establish equivalence between (i) and (iv).Let and let X be the set from (iv).For easier reference, let ξ S (x) . The following three probabilistic events are identical: Therefore, if x ∈ X , then the random variable ξ S (x) is equal to zero with probability 1, and hence x ∈ X .Let us now establish the reverse inclusion.First, let 1 [ξ S (x)≥t] be the indicator function of the event [ξ S (x) ≥ t].Note that since ξ S (x) is a nonnegative random variable, for all t ∈ R we have the inequality Now take x ∈ X and consider t > 0. By taking expectations in (38), we obtain which implies that Prob(ξ S (x) ≥ t) = 0. Now choose t i = 1/i for i = 1, 2, . . .and note that the event [ξ S (x) > 0] can be written as Therefore, by the union bound, Prob(ξ S (x) > 0) = 0, which immediately implies that Prob(ξ S (x) = 0) = 1.From (37) we conclude that x ∈ X .That X does not depend on B follows from representation (iv).

Exactness of the Reformulations
In this section ask the following question: when are the stochastic formulations (2), ( 3), ( 4), ( 5) equivalent to the linear system (1)?This leads to the concept of exactness, captured by the following assumption.
That is, X = L.
We do not need this assumption for all our results, and hence we will specifically invoke it when needed.For future reference in the paper, it will be useful to be able to draw upon several equivalent characterizations of exactness.Theorem 3.6 (Exactness).The following statements are equivalent: Proof.Choose any x * ∈ L. We know that L = x * + Null (A).On the other hand, Theorem 3.4 says that X = x * + Null (E [Z]).This establishes equivalence of (i) and (ii).If (ii) holds, then ), proving (ii).We now show that (ii) and (iv) are equivalent.First, note that ).It remains to combine these observations.We now list two sufficient conditions for exactness.Lemma 3.7 (Sufficient conditions).Any of these conditions implies that Assumption 3.5 is satisfied: , we have exactness by applying Theorem 3.6.Finally, (ii) implies statement (iv) in Theorem 3.6, and hence exactness follows.

Basic Method
We propose solving (33) by Algorithm 1.In the rest of this section we offer several equivalent interpretations of the method.

Algorithm 1 Basic Method
Remark 1.Since S is random, the matrices H = S(S AB −1 A S) † S and Z = A HA are also random.At iteration k of our algorithm, we sample matrix S k ∼ D and perform an update step.It will be useful to to use the notation In the rest of this section we analyze Algorithm 1.

Condition number of the stochastic reformulation
Recall that the Hessian of f is given by Since B −1 E [Z] is not symmetric (although it is self-adjoint with respect to the B-inner product), it will be more convenient to instead study the spectral properties of the related matrix Note that this matrix is symmetric, and has the same spectrum as be the eigenvalue decomposition of W, where U = [u 1 , . . ., u n ] ∈ R n×n is an orthonormal matrix composed of eigenvectors, and Λ = Diag (λ 1 , λ 2 , . . ., λ n ) is a diagonal matrix of eigenvalues.Assume without loss of generality that the eigenvalues are ordered from largest to smallest: We shall often write λ max = λ 1 to denote the largest eigenvalue, and λ min = λ n for the smallest eigenvalue.
Proof.Since B −1/2 ZB −1/2 is symmetric positive semidefinite, so is its expectation W, implying that λ i ≥ 0 for all i.Further, note that B −1/2 ZB −1/2 is a projection matrix.Indeed, it is the projection (in the standard I-norm) onto Range B −1/2 A S .Therefore, its eigenvalues are all zeros or ones.Since the map X → λ max (X) is convex, by Jensen's inequality we get It follows from Assumption 3.5 that λ max > 0. Indeed, if we assume that λ i = 0 for all i, then from Theorem 3.6 and the fact that Null (W) = Range (u i : λ i = 0) we conclude that Null AB −1/2 = R n , which in turn implies that Null (A) = R n .This can only happen if A = 0, which is a trivial case we excluded from consideration in this paper by assumption.Now, let j be the largest index for which λ j > 0. We shall often write λ + min = λ j .If all eigenvalues {λ i } are positive, then j = n.
We now define the condition number of problem (32) to be the quantity Lemma 4.2 (Quadratic bounds).For all x ∈ R n and x * ∈ L we have and Moreover, if Assumption 3.5 holds, then for all x ∈ R n and x * = Π B L (x) we have Proof.In view of ( 34) and ( 40), we obtain a spectral characterization of f : where x * is any point in L. On the other hand, in view of ( 35) and ( 40), we have Inequality (42) follows by comparing ( 45) and ( 46), using the bounds λ + min λ i ≤ λ 2 i ≤ λ max λ i , which hold for i for which λ i > 0.
We now move to the bounds involving norms.First, note that for any x * ∈ L we have The upper bound follows by applying the inequality , then in view of (21), we have B 1/2 (x − x * ) ∈ Range B −1/2 A .Applying Lemma B.1 to (47), we get the lower bound.
Remark 2. Bounds such as those in Lemma 4.2 are often seen in convex optimization.In particular, if φ : R n → R is a µ-strongly convex and L-smooth function, then µ(φ(x) − φ * ) ≤ 1 2 ∇φ(x) 2 ≤ L(φ(x) − φ * ) for all x ∈ R n , where φ * = min x φ(x).In our case, the optimal objective value is zero.The presence of B-norm is due to us defining gradients using the B-inner product.Moreover, it is the case that f is λ max -smooth, which explains the upper bound.However, f is not necessarily µ-strongly convex for any µ > 0, since E [Z] is not necessarily positive definite.However, we still obtain a nontrivial lower bound.

Convergence of expected iterates
We now present a fundamental theorem precisely describing the evolution of the expected iterates of the basic method.1.Let x * ∈ L be chosen arbitrarily.Then Moreover, by transforming the error via the linear mapping h → U B 1/2 h, this can be written in the form which is separable in the coordinates of the transformed error: Finally, 2. Assumption 3.5 hold and let x * = Π B L (x 0 ).Then for all i = 1, 2, . . ., n, Moreover, where the rate is given by ρ(ω) Note that all eigenvalues of W play a role, governing the convergence speeds of individual elements of the transformed error vector.Under exactness, and relative to the particular solution x * = Π B L (x 0 ), the expected error E [x k − x * ] converges to zero at a linear rate.The proof of the theorem is provided in Section 4.3.
Remark 3. Having established (48), perhaps the the most obvious way of analyzing the method is by taking B-norms on both sides of identity (48).This way we obtain the estimate and σ max (•) denotes the largest singular value.This gives the inequality which can be directly compared with (51).We now highlight two differences between these two bounds.
The first approach gives a more detailed, information, as the identity in ( 51) is an exact formula for the norm of the expected error.Moreover, while in view of (54), we have ρ(ω) = max i: The two bounds are identical if λ min > 0, but they differ otherwise.In particular, as long as λ min = 0, we have ρ(ω) ≥ 1 for all ω, which means that the bound (55) does not guarantee convergence.
The following result, characterizing convergence of the expected errors to zero, is a straightforward corollary of Theorem 4.3.
Corollary 4.4 (Necessary and sufficient conditions for convergence).Let Assumption 3.5 hold.Choose any x 0 ∈ R n and let x * = Π B L (x 0 ).If {x k } are the random iterates produced by Algorithm 1, then the following statements are equivalent: We first start with a lemma.
Lemma 4.5.Let Assumption 3.5 hold.Consider arbitrary x ∈ R n and let We now proceed with the proof of Theorem 4.3.The iteration of Algorithm 1 can be written in the form e k+1 = ( where e k = x k −x * .Multiplying both sides of this equation by B 1/2 from the left, and taking expectation conditional on e k , we obtain Taking expectations on both sides and using the tower property, we get We now replace B −1/2 E [Z] B −1/2 by its eigenvalue decomposition UΛU (see (40)), multiply both sides of the last inequality by U from the left, and use linearity of expectation to obtain Unrolling the recurrence, we get (49).When this is written coordinate-by-coordinate, (50) follows.Identity (51) follows immediately by equating standard Euclidean norms of both sides of (49).If x * = Π B L (x 0 ), then from Lemma 4.5 we see that Using this in (50) gives (52).Finally, inequality (53) follows from ≤ ρ k (ω) The last identity follows from the fact that i u i u i = UU = I.

Choice of the stepsize / relaxation parameter
We now consider the problem of choosing the stepsize (relaxation) parameter ω.In view of ( 53) and ( 54), the optimal relaxation parameter is the one solving the following optimization problem: In the next result we solve the above problem.
where j is such that λ j = λ + min .Note that ρ j (ω) = ρ n (ω) for ω ∈ {0, ω * }.From this we deduce that ρ j ≥ ρ n on (−∞, 0], ρ j ≤ ρ n on [0, ω * ], and ρ j ≥ ρ n on [ω * , +∞), obtaining (58).We see that ρ is decreasing on (−∞, ω * ], and increasing on [ω * , +∞).The remaining results follow directly by plugging specific values of ω into (58).Theorem 4.6 can be intuitively understood in the following way.By design, we know that λ max ≤ 1.If we do not have a better bound on the largest eigenvalue, we can simply choose ω = 1 to ensure convergence.If we have a stronger bound available, say λ max ≤ U < 1, we can pick ω = 1/U , and the convergence rate will improve.The better the bound, the better the rate.However, using a stepsize of the form ω = 1/U where U is not an upper bound on λ max is risky: if we underestimate the eigenvalue by a factor of 2 or more, we can not guarantee convergence.Indeed, if U ≤ λ max /2, then 1/U ≥ 2/λ max and hence ρ(ω) ≥ 1.Beyond this point, information about λ + min is useful.However, the best possible improvement beyond this only leads to a further factor of 2 speedup in terms of the number of iterations.Therefore, one needs to be careful about underestimating λ max .
Example 2 (Random vectors).An important class of methods is obtained by restricting S to random vectors.In this case, and thus ω * = 2/(λ + min + λ max ) ≥ 2. This means that in this case we can always safely choose the relaxation parameter to be ω = 2.This results in faster rate than the choice ω = 1.

L2 convergence
In this section we establish a bound on E x k − x * 2 B , i.e., we prove L2 convergence.This is a stronger type of convergence than what we get by bounding B .Indeed, for any random vector x k we have the inequality (see Lemma 4.1 in [17]) Hence, L2 convergence also implies that the quantity We shall first establish an insightful lemma.The lemma connects two important measures of success: Lemma 4.7.Choose x 0 ∈ R n and let {x k } ∞ k=0 be the random iterates produced by Algorithm 1, with an arbitrary relaxation parameter ω ∈ R. Let x * ∈ L. Then we have the identities Moreover, and Proof.Recall that Algorithm 1 performs the update . From this we get In a similar vein, = ( establishing (62).Taking expectation in (64), we get Taking expectation in (62), we get It remains to take expectation again.
In our next result we utilize Lemma 4.7 to establish L2 convergence of the basic method.Theorem 4.8 (L2 convergence).Let Assumption 3.5 hold and set x * = Π B L (x 0 ).Let {x k } be the random iterates produced by Algorithm 1, where the relaxation parameter satisfies 0 < ω < 2.
(i) For k ≥ 0 we have (ii) The average iterate The best rate is achieved when ω = 1.
(ii) By summing up the identities from (63), we get 2ω(2 − ω) k−1 t=0 φ t = r 0 − r k .Therefore, Not that in part (i) we give both an upper and a lower bound on

Convergence of expected function values
In this section we establish a linear convergence rate for the decay of E [f (x k )] to zero.We prove two results, with different quantitative (speed) and qualitative (assumptions and insights gained) qualities.
The complexity in the first result (Theorem 4.9 ) is disappointing: it is (slightly) worse than quadratic in the condition number ζ.However, we do not need to invoke Assumption 3.5 (exactness).In addition, this result implies that the expected function values decay monotonically to zero.
The optimal rate is achieved for ω = 1/ζ, in which case we get the bound Proof.Let S ∼ D be independent from S 0 , S 1 , . . ., S k and fix any x * ∈ L. Then we have Taking expectations, conditioned on x k (that is, the expectation is with respect to S k ), we can further write where We shall now bound α k from below and β k from above in terms of f (x k ).Using the inequality E [Z] λ max B (this follows from B −1/2 E [Z] B −1/2 λ max I), we get = λ max f (x k ).
On the other hand, ≥ 2λ + min f (x k ).Substituting the bounds for α k and β k into (68), we get . It remains to unroll the recurrence.
We now present an alternative convergence result (Theorem 4.10), one in which we do not bound the decrease in terms of the initial function value, f (x 0 ), but in terms of a somewhat larger quantity.This allows us to provide a better convergence rate.For this result to hold, however, we need to invoke Assumption 3.5.Note also that this result does not imply that the expected function values decay monotonically.
Theorem 4.10 (Convergence of expected function values).Choose x 0 ∈ R n , and let {x k } ∞ k=0 be the random iterates produced by Algorithm 1, where the relaxation parameter satisfies 0 < ω < 2.
(ii) Now let Assumption 3.5 hold.For x * = Π B L (x 0 ) and k ≥ 0 we have The best rate is achieved when ω = 1. Proof.
Remark 4. Theorems 4.9 and 4.10 are complementary.In particular, the complexity result given in Theorem 4.9 (for the last iterate) holds under weaker assumptions.Moreover, Theorem 4.9 implies monotonicity of expected function values.On the other hand, the rate is substantially better in Theorem 4.10.Also, Theorem 4.10 applies to a wider range of stepsizes.
It is also possible to obtain other convergence results as a corollary.For instance, one can get a linear rate for the decay of the norms of the gradients as a corollary of Theorems 4.9 and 4.10 using the upper bound in Lemma 4.2.

Parallel and Accelerated Methods
In this section we propose and analyze parallel and accelerated variants of Algorithm 1.

Algorithm 2 Parallel Method
Average the results For brevity, we only prove L2 convergence results.However, various other results can be obtained as well, as was the case for the basic method, such as convergence of expected iterates, expected function values and average iterates.
Theorem 5.1.Let Assumption 3.5 hold and set x * = Π B L (x 0 ).Let {x k } ∞ k=0 be the random iterates produced by Algorithm 2, where the relaxation parameter satisfies 0 < ω < 2/ξ(τ ), and ξ(τ where ρ(ω, τ For any fixed τ ≥ 1, the optimal stepsize choice is ω(τ ) def = 1/ξ(τ ) and the associated optimal rate is Proof.Recall that Algorithm 2 performs the update Next, we can write Since Z ki B −1 Z ki = Z ki , and because Z ki and Z kj are independent for i = j, we can further write where we have used the estimate E , which follows from the bound W 2 ≤ λ max W. Plugging (74) into (73), and noting that x k − x * 2 E[Z] = 2f (x k ), we obtain: ≤ ρ(ω, τ ) The inequality involving f is shown in the same way as in Theorem 4.10.

Accelerated method
In this section we develop an accelerated variant of Algorithm 1. Recall that a single iteration of Algorithm 1 takes the form x k+1 = φ ω (x k , S k ), where We have seen that the convergence rate progressively improves as we increase ω from 1 to ω * , which is the optimal choice.In particular, with ω = 1 we have the complexity Õ(1/λ + min ), while choosing ω = 1/λ max = 1/λ max or ω = ω * leads to the improved complexity Õ(λ max /λ + min ) = Õ(ζ).In order to obtain further acceleration, we suggest to perform an update step in which x k+1 depends on both x k and x k−1 .In particular, we take two dependent steps of Algorithm 1, one from x k and one from x k−1 , and take an affine combination of the results.This, the process is started with x 0 , x 1 ∈ R n , and for k ≥ 1 involves an iteration of the form where the matrices {S k } are independent samples from D, and γ ∈ R is an acceleration parameter.
Note that by choosing γ = 1 (no acceleration), we recover Algorithm 1.This method is formalized as Algorithm 3.

Algorithm 3 Accelerated Method
1: Parameters: distribution D from which to sample matrices; positive definite matrix B ∈ R n×n ; stepsize/relaxation parameter ω > 0; acceleration parameter γ > 0 As we shall see, by a proper combination of overrelaxation (choice of ω) with acceleration (choice of γ), Algorithm 3 enjoys the accelerated complexity of Õ( √ ζ).We start with a lemma describing the evolution of the expected iterates.
Proof.Taking expectation on both sides of (78), we get After subtracting x * from both sides, using (77), and replacing b by Ax * , we get where Z k = A S k (S k AB −1 A S k ) † S k A. We now use the tower property and linearity of expectation, we finally obtain: We can now state our main complexity result.Note that the optimal choice of parameters, covered in case (i), leads to a rate which depends on the square root of the condition number.
(3) 0 < ωλ i < 1.In this case we have where the last inequality follows from the assumption µ < ωλ + min .Therefore, we can apply Lemma C.1, using which we can deduce the bound As |C 0 | + |C 1 | depends on i, we shall write Putting everything together, for all k ≥ 2 we have finishing the proof.
Note that we do not have a result on L2 convergence.We have tried to obtain an accelerated rate in the L2 sense, but were not successful.We conjecture that such a result can be obtained.

Conclusion
We have developed a generic scheme for reformulating any linear system as a stochastic problem, which has several seemingly different but nevertheless equivalent interpretations: stochastic optimization problem, stochastic linear system, stochastic fixed point problem, and probabilistic intersection problem.
While stochastic optimization is a broadly studied field with rich history, the concepts of stochastic linear system, stochastic fixed point problem and probablistic intersection appear to be new.
We give sufficient, and necessary and sufficient conditions for the reformulation to be exact, i.e., for the solution set of the reformulation to exactly match the solution set of the linear system.To the bets of our knowledge, this is is first systematic study of stochastic reformulations of linear systems.Further, we have developed three algorithms-basic, parallel and accelerated methods-to solve the stochastic reformulations.We have studied the convergence of expected iterates, L2 convergence, converge of a Cesaro average of all iterates, and convergence of f .Our methods recover an array of existing randomized algorithms for solving linear systems in special cases, including several variants of the randomized Kaczmarz method [60], randomized coordinate descent [24], and all the methods developed in [17,18].
Our work can be extended in several ways.One of the most promising of these is stochastic preconditioning, which refers to the generic problem of fine-tuning the formulations (by designing the distributions D and matrix B) to the structure of A. We conjecture that specific highly efficient methods can be designed in this way.Accelerated convergence in the L2 sense remains an important open problem.
Last but not least, we hope that this work provides bridge across several communities: numerical linear algebra, stochastic optimization, machine learning, computational geometry, fixed point theory, applied mathematics and probability theory.We hope that our work may inspire further progress at the boundaries of these fields.

3 .
Stochastic fixed point problem.Let Π B L S (x) denote the projection of x onto L S def = {x : S Ax = S b}, in the norm x B def = √ x Bx.Consider the stochastic fixed point problem

Assumption 3 . 2 (
Finite mean).The random matrix H has a mean.That is, the matrix E S∼D [H] = E S∼D S(S AB −1 A S) † S has finite entries.

Theorem 4 . 3 (
Convergence of expected iterates).Choose any x 0 ∈ R n and let {x k } be the random iterates produced by Algorithm 1.

Lemma 5 . 2 (
Expected iterates).Let x * be any solution of Ax = b and let r k def = E [x k − x * ].Then for all k we have the recursion
which does not change with k.However, it turns out that under exactness we have u i B 1/2 (x 0 − x * ) = 0 whenever λ i = 0 if we let x * to be the projection, in the B-norm, of x 0 onto L (Theorem 4.3).This is then used to argue (Corollary 4.4) that E [x k − x * ] B converges to zero if and only if 0 < ω < 2/λ max .The choice of stepsize issue is discussed in detail in Section 4.4.The main complexity results obtained in this paper are summarized in Table1.The full statements including the dependence of the rate on these parameters, as well as other alternative results (such as lower bounds, ergodic convergence) can be found in the theorems referenced in the table.L2 convergence.The rate of decay of the quantity E [x k − x * ] 2