An Implicit Representation and Iterative Solution of Randomly Sketched Linear Systems
Vivak Patel, Mohammad Jahangoshahi, Daniel Adrian Maldonado
TL;DR
This paper introduces an implicit, iterative approach to randomized sketching for linear systems that addresses practical issues like sketch size selection and storage costs, while enhancing convergence analysis and rates.
Contribution
It proposes an implicit method for solving sketched linear systems that eliminates the need for pre-determined sketch size and reduces storage costs, connecting sketching with randomized iterative solvers.
Findings
Improved convergence theory for randomized iterative solvers under various sampling schemes.
Enhanced convergence rates with controlled computational and storage costs.
Validated approach on forty-nine different linear systems.
Abstract
Randomized linear system solvers have become popular as they have the potential to reduce floating point complexity while still achieving desirable convergence rates. One particularly promising class of methods, random sketching solvers, has achieved the best known computational complexity bounds in theory, but is blunted by two practical considerations: there is no clear way of choosing the size of the sketching matrix apriori; and there is a nontrivial storage cost of the sketched system. In this work, we make progress towards addressing these issues by implicitly generating the sketched system and solving it simultaneously through an iterative procedure. As a result, we replace the question of the size of the sketching matrix with determining appropriate stopping criteria; we also avoid the costs of explicitly representing the sketched linear system; and our implicit representation…
| Total Time and Effort Costs to Iteration | ||||||||||
| Platform | Computing | Update Costs | Network | |||||||
| Time | Total Effort | Iterate | Matrix | |||||||
| Sequential | No | |||||||||
| Shared Memory | No | |||||||||
| Distributed Memory | Yes | |||||||||
| Comparison of Estimated Theoretical Rates of Convergence | ||||
| Matrix Name | Estimated Rates by Result | |||
| Theorem 4.8 of Richtárik and Takác (2020) | Proposition 3 | |||
| deriv2 | ||||
| heat | ||||
| randsvd | ||||
| ursell | ||||
| wing | ||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
An Implicit Representation and Iterative Solution of Randomly Sketched Linear Systems
Vivak Patel, Mohammad Jahangoshahi & Daniel Adrian Maldonado
Abstract
Randomized linear system solvers have become popular as they have the potential to reduce floating point complexity while still achieving desirable convergence rates. One particularly promising class of methods, random sketching solvers, has achieved the best known computational complexity bounds in theory, but is blunted by two practical considerations: there is no clear way of choosing the size of the sketching matrix apriori; and there is a nontrivial storage cost of the sketched system. In this work, we make progress towards addressing these issues by implicitly generating the sketched system and solving it simultaneously through an iterative procedure. As a result, we replace the question of the size of the sketching matrix with determining appropriate stopping criteria; we also avoid the costs of explicitly representing the sketched linear system; and our implicit representation also solves the system at the same time, which controls the per-iteration computational costs.
Additionally, our approach allows us to generate a connection between random sketching methods and randomized iterative solvers (e.g., randomized Kaczmarz method, randomized Gauss-Seidel). As a consequence, we exploit this connection to (1) produce a stronger, more precise convergence theory for such randomized iterative solvers under arbitrary sampling schemes (i.i.d., adaptive, permutation, dependent, etc.), and (2) improve the rates of convergence of randomized iterative solvers at the expense of a user-determined increases in per-iteration computational and storage costs. We demonstrate these concepts on numerical examples on forty-nine distinct linear systems.
1 Introduction
Over the past few decades, randomized linear system solvers have become popular as they have the potential to reduce floating point complexity or maintain limited memory footprints, while still achieving desirable convergence rates (e.g., Strohmer and Vershynin, 2009; Woodruff, 2014). In particular, the noniterative class of randomized linear system solvers, based on random matrix sketching (see Woodruff, 2014), have exceptionally low computational complexities, at least in theory. Unfortunately, the theoretical promise of these random matrix sketching solvers is blunted by their practical limitations: there is no clear way of choosing the size of the sketching matrix and there is a nontrivial storage cost of the projected system (Mahoney, 2016). In fact, the practical challenges of random matrix sketching solvers have prevented them from being fully embraced by the numerical optimization community (e.g., Nocedal, 2018).
In this work, we begin to address these two primary practical issues of random matrix sketching, which we recall are: the challenge of choosing the size of the sketching matrix, and the challenge of storing the projected system. Our main insight is to recast the separate sketch-then-solve core of random sketching methods into an equivalent, iterative sketch-and-solve, in which the sketching matrix is generated incrementally without being explicitly stored and the system is incrementally solved from the implicitly derived sketched matrix.111It is worth mentioning that the random sketch solvers have been used iteratively in a different sense (e.g., see Gower and Richtárik, 2015): the noniterative scheme is simply repeated in order to get better convergence properties. We are not doing this, but rather turning the noniterative scheme into an iterative one. As a result of our approach, (1) we can implicitly grow the size of the sketching matrix until a user-determined stopping criteria is reached without having to determine the size of the sketching matrix apriori; (2) we implicitly represent the sketched system without having to explicitly store the projected system, which allows us to avoid the cost of storing the projected system; and (3) we can naturally implement random sketching solvers within distributed and parallel computing paradigms. Thus, our approach of converting the usual sketch-then-solve procedure to a sketch-and-solve procedure begins to address the aforementioned practical challenges of random matrix sketching.
Moreover, our approach provides a bridge between the newer concerns around sketching-based solvers and more classical areas of applied mathematics research such as stopping criteria. One such bridge is the placement of random sketching methods and (what we will call) base randomized iterative methods222We will be more precise about what we refer to as base methods. For now, such methods are exemplified by randomized Kaczmarz (Strohmer and Vershynin, 2009) and randomized Gauss-Seidel (Leventhal and Lewis, 2010). on a single spectrum of procedures, which has several immediate consequences.
First, the number of rows of the sketching matrix that results in the solution (this number is a random quantity) connects to an alternative rate-of-convergence result for general base randomized iterative methods that guarantees a rate-of-convergence less than one for arbitrary sampling schemes—even for underdetermined systems (Theorem 5). Consequently, our results complement and improve on previous results in several ways. In particular, we allow for arbitrary sampling schemes, not just sampling schemes that are independent and identically distributed as in Gower and Richtárik (2015) (Lemma 4.2), Richtárik and Takác (2020) (Theorem 4.8), Zouzias and Freris (2013) (Theorem 3.4), and Ma et al. (2015) (Equation 3.10). Moreover, our results do away with the exactness assumption (see Richtárik and Takác, 2020, Assumption 2), and precisely characterize the inexactness that can occur for arbitrary sampling schemes (Theorems 5 and 6). Additionally, our results define convergence on a maximal set—effectively, a set occurring with probability one for sampling schemes of interest—, which builds on the work of Chen and Powell (2012). As example applications of our results, we supply rates of convergence with probability one for random permutation sampling methods (Proposition 2) and independent, identically distributed sampling schemes (asymptotically, see Proposition 3). As a more interesting application of our results, we specify generic conditions for the convergence of a broad class of adaptive schemes (see Subsection 4.5), which can account for the maximum residual scheme, the maximum distance scheme, schemes that randomize over a greedy subset, and schemes that are greedy over randomized subsets (Motzkin and Schoenberg, 1954; Gubin et al., 1967; Lent, 1976; Censor, 1981; Nutini et al., 2016; Bai and Wu, 2018; Haddock and Ma, 2019). We note that the rates that we provide as examples are rather loose in comparison to results that are specialized to each case, yet our results often supply information that is not available in these other results as discussed above.
Second, we can generate a series of “intermediate” procedures between sketching methods and base methods that trade-off between computational resources (e.g., floating-point operations, storage) and rates of convergence. Thus, we can take a sketching method and reduce its computational footprint in exchange for a slower rate of convergence, or increase the computational footprint of base methods to improve their rate of convergence (Algorithm 2). Moreover, these “intermediate” procedures can be readily parallelized as we discuss in Section 2.
Finally, by shifting our perspective from improving the sketch-then-solve procedure to improving the performance of base methods, we find that our approach is a randomized orthogonalization procedure in the row space of the coefficient matrix of the linear system. Thus, by presenting our approach from this latter perspective, we will simplify the introduction and the related theory of our approach. Now, before pursuing this further, we reiterate our main contributions.
First, we turn the typical sketch-then-solve noniterative random sketching solver into an iterative, sketch-and-solve method, which lays a foundation for addressing the previously enumerated practical challenges of random sketching solvers: there is no clear way of choosing the size of the sketching matrix apriori; and there is a nontrivial storage cost of the sketched system. 2. 2.
Second, through our approach, we place random sketching methods and base randomized iterative methods (e.g., randomized Kaczmarz, randomized Gauss-Seidel, and Sketch-and-Project (Gower and Richtárik, 2015)) on a single spectrum of methods. 3. 3.
Third, owing to this connection, we are able to generate “intermediate” methods between random sketching and base methods, which can trade-off between computational resources and rates of convergence. 4. 4.
Fourth, owing to this connection, we use the geometric implications of random sketching methods to develop an alternative rate-of-convergence result for general base methods for arbitrarily determined systems and arbitrary sampling schemes, which advances the with-probability-one results of Chen and Powell (2012), generalizes the deterministic cyclic results in Bai and Liu (2013); Galántai (2005); Wallace and Sekmen (2014), complements the mean-squared error results of Richtárik and Takác (2020), and accounts for a litany of adaptive methods considered in Motzkin and Schoenberg (1954); Gubin et al. (1967); Lent (1976); Censor (1981); Nutini et al. (2016); Bai and Wu (2018); Haddock and Ma (2019). 5. 5.
Finally, we provide a generic set of conditions for characterizing a broad class of adaptive methods, and, from these conditions, prove convergence and rate-of-convergence results for a number of classical and emerging adaptive methods in the literature under a unified framework (see Subsection 4.5).
The remainder of this paper is organized as follows. In Section 2, we introduce our procedure; we state the connection between our procedure and random sketching methods, which allows us to convert the less practical sketch-then-solve approach to our sketch-and-solve approach; and, finally, we introduce our general algorithm and variants for low-memory environments, shared memory environments, distributed memory environments, and large, sparse, structured linear systems. In Sections 3 and 4, we develop the convergence theory for the two methodological extremes—sketching and base methods—leaving the intermediate, more complex cases to future work, and discuss particular examples. In Section 5, we test our algorithms on forty-nine distinct linear systems. In Section 6, we conclude this work and preview future efforts.
2 Our Procedure
While our motivating application is to address the practicality of random sketching methods, our approach is best introduced from the perspective of base randomized iterative methods. Here, we review the basic formulation of randomized iterative methods (Subsection 2.1), which we then use to heuristically introduce our general procedure (Subsection 2.2). We then refine our procedure for the case of rank-one methods, such as Randomized Kaczmarz and Randomized Gauss-Seidel, which allows us to restate random sketching from a sketch-then-solve procedure to a sketch-and-solve procedure (Subsection 2.3). We conclude this section with comments on algorithmic refinements for parallel platforms (Subsection 2.4.1), limited memory platforms (Subsection 2.4.2), and, for structured linear systems, limited communication platforms (Subsection 2.4.3).
2.1 A Brief Overview
Let and be the coefficient matrix and constant vector, respectively. Assuming consistency, our goal is to determine an , not necessarily unique, such that
[TABLE]
In a base randomized iterative approach, a sequence of iterates is generated that has the form
[TABLE]
where are independent random variables, which we call residual projection matrices (RPM). The RPM defines the base technique which is being used. To make this formulation concrete, we give several examples of randomized iterative methods that have this formulation.
Randomized Kaczmarz.
Let denote the row of and let denote the standard basis vector of dimension . Define the random variable such that
[TABLE]
Now, given an independent copy of at each , define the RPM, Then, using Eq. 2,
[TABLE]
which is the Randomized Kaczmarz method of Strohmer and Vershynin (2009).
Randomized Gauss-Seidel.
Let denote the column of and let denote the standard basis vector of dimension . Define a random variable such that
[TABLE]
Now, given an independent copy of at each , define the RPM, Then, using Eq. 2,
[TABLE]
which is the Randomized Gauss-Seidel method of Leventhal and Lewis (2010).
Randomized Block Coordinate Descent.
Let be a subset of . Let whose columns are the -dimensional standard basis vectors whose non-zero components correspond to the indices in . Let be a partition , and define a random variable that randomly selects a partition in . Given an independent copy of at each , define the RPM, . Then, using Eq. 2,
[TABLE]
which is a version of the randomized block coordinate descent method specified by (Gower and Richtárik, 2015, Equation 3.14).
Sketch-and-Project.
Let be a sequence of sketching matrices with columns. Define the RPM to be . Then, using Eq. 2,
[TABLE]
which is the general sketch-and-project method (Gower and Richtárik, 2015, Equation 2.2).
2.2 A Heuristic Derivation
Here, given a strategy for defining , we consider how to augment the randomized iterative method with prior information in order to improve convergence. For this purpose, we propose defining a sequence of matrices (discussed below) and modify Eq. 2 to be
[TABLE]
Of course, can simply be absorbed by ; however, our goal is to augment a randomized iterative method. For this reason, we will keep these two quantities separate.
The main question now is how to choose . Our guiding principle is that should minimize some measure of error between and . However, implementing this guiding principle requires (1) choosing an appropriate error measure and (2) handling the fact that is unknown. In order to convey the intuition behind our procedure, we now state the heuristics that we use to make these choices.
Choosing an Error Measure.
Temporarily, suppose is known, and suppose we choose the error as our measure. Then, we must minimize the difference between the next iterate and . While this error metric might have merit, solving it is a convex optimization problem that is as difficult to solve as the original linear system. Therefore, we will need an error measure which gives an explicit representation for . Hence, one sensible choice is to use the Mahalanobis norm,
[TABLE]
where is a positive definite, symmetric matrix.
Compensating for the Unknown Solution.
Now, we consider the task of compensating for the unknown . For a fixed and for all , let . Then, is related to by
[TABLE]
where we have made use of Eq. 3. Using Eq. 5, we can rewrite Eq. 4 as
[TABLE]
To find an optimal , we differentiate the right hand side and set the quantity equal to zero, which, explicitly is
[TABLE]
Clearly, is positive semi-definite, so the solution to such a system will be the minimizer of the original objective function. However, Eq. 6 may have many possible solutions or may fail to be consistent. In the case of nonunique solutions, we arbitrarily choose the solution with the smallest Frobenius norm. In the case of an inconsistent system, we arbitrarily choose the solution that minimizes the Frobenius norm of the residual and has the minimal Frobenius norm. In both cases, a straightforward calculation gives
[TABLE]
where represents the Moore-Penrose Pseudo-inverse. Using Eq. 7 with Eq. 5, we have the following recursion
[TABLE]
From Eq. 7 and Eq. 8, it is clear that if were known, then the remaining unknown quantities could be determined.
Our Procedure.
Since is unknown, we use the following heuristic procedure instead. First, we let , where is the -dimensional identity matrix. Then, we recursively define and according to Eqs. 7 and 8. To summarize, given , we let , let , and define
[TABLE]
where
[TABLE]
and
[TABLE]
To interpret the terms in the above procedure, we begin by ignoring (i.e., set it to the identity). In this case, and its role in updating to is familiar: serves to map the residual onto the row space of , thereby ensuring that satisfies . If we now consider the role of , we see that it is an orthogonal projector that “weights” the behavior of to ensure that satisfy for . We will see these interpretations clearly and formally when we focus on the case of rank-one next.
We pause here momentarily to discuss the relationship between our procedure, as specified by Eqs. 9, 10 and 11, and the sketch-and-project method in Gower and Richtárik (2015) and Richtárik and Takác (2020). At first glance, it may seem that our procedure is a special case of sketch-and-project with adaptive choices of the inner product at each iteration of the sketch-and-project update. Unfortunately, an effort to recast our approach as a special case of sketch-and-project breaks down at two fundamental points. First, the adaptive choices of the sketch-and-project inner product would have to be the inverse of , which are orthogonal projection matrices. As a result, the inverse is ill-defined and the inner product is ill-defined. Of course, this can be rectified by allowing for a pseudo-metric, but this then results in the second major point of difficulty: the theory presented in Gower and Richtárik (2015) and Richtárik and Takác (2020) relies on the determinism and invertibility of the matrix defining the metric space to prove convergence. Thus, sketch-and-project, without a substantial investment, cannot readily include our approach. On the other hand, we can state sketch-and-project as a base randomized iterative approach, as shown in Subsection 2.2, and then improve on it with our procedure via Eqs. 9, 10 and 11.
2.3 Rank-One Refinements and Random Sketching
By choosing and , Eqs. 9, 10 and 11 describe an orthogonal projection procedure for typical randomized iterative procedures. However, because our goal is to improve the practicality of random sketching methods, we will need to focus on a particular refinement of the general procedure that occurs when are rank-one matrices, that is, when there exist pairs of vectors such that for each . In this case, Eqs. 10 and 11 become
[TABLE]
and
[TABLE]
Moreover, if we substitute Eq. 12 into Eq. 9, we recover
[TABLE]
It follows from Eqs. 14 and 13 that in the case of a rank-one RPM, the left singular vector of the RPM is not important. To give some explicit examples, recall that rank-one RPM methods include the important special cases of randomized Kaczmarz and Gauss-Seidel.
Randomized Kaczmarz with Orthogonalization.
Let denote the row of and let denote the standard basis vector of dimension . Define the random variable arbitrarily taking values in . Now, given an independent copy of at each , the randomized Kaczmarz method has rank-one RPM, Then, using Eqs. 14 and 13, the randomized Kaczmarz method with orthogonalization is
[TABLE]
when , or is and otherwise.
Randomized Gauss-Seidel with Orthogonalization.
Let denote the column of and let denote the standard basis vector of dimension . Define a random variable arbitrarily taking values in . Now, given an independent copy of at each , the randomized Gauss-Seidel method has rank-one RPM, Then, using Eqs. 14 and 13, the randomized Gauss-Seidel method with orthogonalization is
[TABLE]
when , or is and otherwise.
Again, we see from the two preceding examples that the left singular vector of the rank-one RPM does not play a role in the updates for our procedure. As we now explain, this observation is critical for converting the impractical, noniterative randomized sketch-then-solve methods into iterative randomized sketch-and-solve methods.
Recall that the fundamental sketch-then-solve procedure is to construct a specialized matrix , then generate and solve the smaller, sketched problem (see Woodruff, 2014, Ch. 1).333We note that the typical formulation considers linear regression rather than a linear system. The special matrix , called the sketching matrix, can be generated in a variety of ways such as making each entry an independent, identically distributed Gaussian random variable (Indyk and Motwani, 1998), or by setting the columns of as uniformly sampled columns (with replacement) of the appropriately-dimensioned identity matrix (Cormode and Muthukrishnan, 2005).
In order to convert the usual sketch-then-solve procedure into our sketch-and-solve procedure, we simply set to the transposed rows of , which we will rigorously demonstrate in Section 3. Of course, this requires that we have a streaming procedure for generating arbitrarily many rows of . For concreteness, we show how to do this for the two sketching strategies just mentioned.
Random Gaussian Sketch.
In the random Gaussian sketch, the entries of the sketching matrix, , are independent, standard normal random variables. Accordingly, we let be independent, -dimensional standard normal vectors. We see that if has rows, then and
[TABLE]
have the same distribution.
Count Sketch.
Fix , and let be drawn from the standard basis vectors with replacement. Define a sequence of Rademacher random variables which are independent and independent of . The count sketch sketching matrix, , is specified by
[TABLE]
which is a matrix whose entries are either , [math] or . Generally, the choice of is the topic of substantial theory and consideration (Cormode and Muthukrishnan, 2005; Clarkson and Woodruff, 2017). Owing to the fact that we have a streaming procedure, we do not need to worry too much about . Therefore, we generate as follows:
Generate a count sketch matrix with small. In our experiments below, we used . 2. 2.
To generate a , pop a row of the matrix and set it to . 3. 3.
Once the count sketch matrix is exhausted, regenerate a new count sketch matrix with the same . Repeat.
From this strategy, (a) if we let denote a sequence of independent count sketch matrices, (b) denote the remainder of an integer divided by and incremented by one, and (c) we let denote the standard basis vectors of , then for all .
Thus, if we let RPMStrategy() define a generic user-defined procedure for choosing , then this observation gives us Algorithm 1 for (1) converting the sketch-then-solve procedure into a sketch-and-solve procedure, and (2) adding orthogonalization to such base methods as randomized Kaczmarz and randomized Gauss-Seidel.
2.4 Algorithmic Refinements Considering the Computing Platform
Algorithm 1 implicitly assumes the traditional sequential programming paradigm. However, the performance of the algorithm can be improved by taking advantage of parallel computing architectures. Here, we will consider a handful of important computing architecture abstractions and how our procedure can adapt to different configurations. In Subsection 2.4.1, we will consider the case of a parallel computing architecture for which the communication overhead, which is proportional to the dimension , is not a limiting factor. For this subsection, the problems that we have in mind come from data and imaging sciences, where and is reasonably sized. In Subsection 2.4.2, we consider a similar class of problems where the communication of -sized vectors is acceptable and , but that is so large that storing and manipulating a matrix in is burdensome. Finally, in Subsection 2.4.3 we will consider problems in which computational overhead becomes a bottleneck for scalability, but that we have structured systems that will allow us to circumvent this issue. For this ultimate subsection, the problems that we have in mind here come from the solution of systems of differential equations (e.g., Dongarra and Sørensen, 1986).
2.4.1 Asynchronous Parallelization on Shared and Distributed Memory Platforms
First, when we are using a matrix sketch for RPMStrategy(), one of the expensive components of the computation is determining . Fortunately, in our sketch-and-solve procedure, this expensive computation can be trivially asynchronously parallelized on a shared memory platform when
the data within the rows are stored together, and 2. 2.
the RPMStrategy() generates that are either independent (e.g., the Gaussian Strategy) or can be grouped into independent subsets (e.g., the Count-Sketch strategy).
When these two requirements are met, each processor can generate its own independently of the other processors, and evaluate . It can then simply write the resulting row to an address reserved for performing the iterate and matrix updates by the master processor. Importantly, this procedure does not require locking any of the rows of , and the reserved addresses can use fine grained locks to prevent any wasted calculations.
Similarly, in our sketch-and-solve procedure, computing can be trivially asynchronously parallelized on a distributed memory platform using a Fork-join model, when
the rows of are distributed across the different storages, and 2. 2.
the RPMStrategy() generates such that have independent groups of components (e.g., the Gaussian Strategy and the Count-Sketch strategy).
When these two requirements are met, each processor can generate its own and operate on the local rows of . It can then simply pass the resulting row to the master processor which performs the iterate and matrix updates. For each iteration, a scattering and gathering of the data is performed but no other data exchange is required.
Table 1 summarizes the time and total computational costs of computing and from and in the following context: (1) the sequential platform refers to the case where there is a single processor with a sufficiently large memory to store the system, and perform the necessary operations in Algorithm 1; (2) the shared memory platform assumes that there are processors that share a sufficiently large memory. One of the processors is dedicated to performing the iterate and matrix updates, while the remaining processors compute ; (3) the distributed memory architecture assumes that there are processors each with a sufficient memory capacity. The rows of are split evenly or nearly evenly amongst of the processors, and each process only manipulates its local information about and . Finally, master processor is dedicated to performing the iterate and matrix updates.
2.4.2 Memory-Reduced Procedure
Another notable aspect of Algorithm 1 (and its aforementioned parallel variants described above) is that it must store and manipulate the matrix at each iteration, which is clearly expensive when is large or is excessive when is comparable to or greater than . This difficulty motivates a partial orthogonalization approach, as described in Algorithm 2. In this approach, a user-defined parameter specifies the number of -dimensional vectors needed to implicitly store an approximate representation of (based on Theorem 1). With this implicit representation, the cost of computing reduces to ,444If replace in the calculation of , then the cost of computing is (see Golub and Van Loan, 2012, Ch. 5.2). which, consequently, reduces the overall cost of updating to to . Moreover, because is implicitly represented by a dimesnional vectors in , there is no notable additional computational cost incurred for updating to . Thus, an entire iteration incurs a computational cost plus the cost of computing , which can be mollified under the strategies above in shared memory or distributed memory platforms.
Remark 1**.**
Algorithm 2* is an efficient implementation of the partial orthogonalization procedure and, as a result, at , seems to only recover row-action base randomized iterative methods as specified by Eq. 51. A less efficient algorithm based on directly applying Eqs. 12 and 13 with the appropriate low memory modification would recover all rank-one base randomized iterative methods when .*
2.4.3 Optimizing Communication Overhead. Structured Systems
In the above approaches, we take for granted that is not so large such that communicating vectors is acceptable during the procedure. However, for many problems coming from the solution of differential equations (e.g., see Dongarra and Sørensen, 1986), and are of the same order and are so large that communicating vectors at arbitrary points during the procedure is impossible. Fortunately, linear system problems in this class are highly sparse and structured (Saad, 2003, Ch. 2). A simple example is the case where is a square, banded system with nonzero bandwidth for some ; that is, if and the remaining can take arbitrary values.
For such sparse and structured problems, our methodology can be efficiently implemented across a distributed memory platform with processors under some additional qualifications. However, to understand these qualifications, let us first introduce some notation and concepts that define the communication pattern across the nodes.
Suppose somehow that we distribute the equations of our linear system of interest across nodes. Figure 1 shows how the coefficient matrix of a banded system with bandwidth can be distributed across five nodes. Note, in this example, the entries of the constant vector would be stored on the same processor as the corresponding rows of the coefficient matrix. Moreover, we need a way of tracking which components of are manipulated by each node: let be the set of indices of the components of with nonzero coefficients at node in the distributed system for . In our example, , , , , and . Finally, for any vector and any set over the indices of , let be the vector whose elements are the elements of indexed by .
From this example and from our discussion in Subsection 2.4.1 of distributing the RPMStrategy(), we can use the local rows of at Node 1 and a Gaussian sketch to generate a such that are arbitrarily valued and . Thus, our vector is highly sparse and can be generated locally on the node. However, following Algorithm 1, the next step of computing requires computing the product between and , which, in a naive implementation, would require storing a dense matrix and computing a global matrix-vector product. Such a required computation raises several concerns, which we detail and address in the following enumeration.
Given that is relatively large to the computing environment, is storing a matrix even feasible? Generally, the answer will be that storing such a matrix is infeasible. However, by exploiting the properties of (see Theorem 1), we will approximately and implicitly store as , which is a collection of orthonormal vectors. 2. 2.
Even if we use in place of , will the resulting implicit matrix-vector product and update of incur prohibitive communication costs? To answer these questions completely, we will need to specify how the implicit matrix-vector product will be computed and how will be stored. Here, we will compute the implicit matrix-vector product by using twice-iterated classical Gram-Schmidt (Algorithm 4), which was shown to be numerically stable in the seminal work of Giraud et al. (2005). Owing to this calculation pattern, we can store in a distributed fashion across the processors, which we detail below along with the communication cost of the synchronization of .
To understand the costs associated with computing from the orthonormal vectors in and the vector , we will characterize the support of (i.e., index set of its nonzero entries).
Lemma 1**.**
Let and let . Let be a set of orthonormal vectors (hence, ), and let for . If denotes the result of Algorithm 4 applied to over the set then
[TABLE]
where .
Proof.
Letting denote the matrix whose columns are elements of the orthonormal set, we recall that classical Gram-Schmidt generates . Thus, twice iterated Gram-Schmidt can be written as
[TABLE]
which is expected in exact arithmetic. Thus, we can consider classical Gram-Schmidt and ignore the iteration to compute the support of . For any ,
[TABLE]
where we use the fact that if then . For a contradiction, suppose such that
[TABLE]
Then, and for . Using the above formula for , , which is a contradiction. ∎
At iteration , Lemma 1 states that the support of will depend on the support of , which, in turn, has elements whose support depend on (a subset of) . Moreover, if has elements whose combined support cover , which will be necessary to solve the system,555Note, if the combined supports of the elements of do not cover all of , then some components of our iterates, will not be updated. it is possible that the support of will be all of (ignoring any trivial independence in the system). Thus, it appears that we will eventually have to store vectors in whose support is all of . Naively, we may think that we need a faithful copy of at each node in the system, which incurs prohibitive communication costs as the support of tends to . While this is true, a careful inspection of Gram-Schmidt and the nonzero patterns of suggest a less naive approach, which we now detail.
We begin by supposing that on a processor , only are stored on the node for every . Immediately, we have eliminated the need for synchronizing all of on each processor. Instead, we need only to synchronize those components of in for all . Thus, we have that our synchronization costs will depend on the maximum overlap, , between two processors, which, formally, is
[TABLE]
Now, we can understand the precise nature of this synchronization by inspecting Algorithm 4. If for some , , then
[TABLE]
From Eq. 20, we see that we must communicate the values of to all nodes such that , and we must communicate the inner products to all nodes. The resulting number of floating point values that must be communicated (counting each replicate to a node individually) during the first iteration of Algorithm 4 is
[TABLE]
where for (see the notation in Lemma 1). For the second iteration of Algorithm 4, we must broadcast inner products that are partially computed (using some ordering that respects the non-associative property of floating point complexity) on each node to the remaining nodes. Thus, the number of floating point values that must be communicated (counting each replicate to a processor individually) to ensure synchronization is
[TABLE]
which we can bound by
[TABLE]
Noting that represents the maximum shared indices between two nodes and that represents the maximum number of nodes that overlap, the first term in the bound can be controlled by the ordering choice of the differential equations that generate the system, but a discussion of this topic is beyond the scope of this work. Algorithm 5 summarizes a simple version of the procedure described here. We can also modify this algorithm to the low memory context of Algorithm 1 by limiting the number of vectors that can be stored in .
3 Convergence Theory for Orthogonalization
Here, we prove that the complete orthogonalization approach (i.e., Algorithm 1) converges to the solution under a variety of sampling RPM strategies. In Subsection 3.1, we establish a collection of core results that are useful in characterizing the behavior of our procedure. A key feature of these core results is that they will rely on a stopping time , which will depend on the random variables . Therefore, in Subsection 3.2, we characterize under common probabilistic relationships between the elements of . All statements hold with probability one unless stated otherwise.
3.1 Core Results
We establish two key results. First, we establish that our procedure is an orthogonalization procedure: that is, the matrices project the current search direction onto a subspace that is orthogonal to previous search directions. Second, we characterize the limit point of our iterates, , in terms of a true solution of the linear system and the subspace generated by the rank-one RPMs, .
Theorem 1**.**
Let be an arbitrary sequence in , and let and for . Now, let and be defined recursively as in Eq. 13. Then, for , is an orthogonal projection matrix onto .
Proof.
We will prove the result by induction. For the base case, , . It follows that is an orthogonal projection onto since and . Now suppose that the result holds for . If then there is nothing to show. Therefore, for the remainder of this proof, suppose .
First, we show that is a projection matrix by verifying that by direct calculation. Making use of the recursive definition of and the induction hypothesis that ,
[TABLE]
Second, we use the fact that a projection is orthogonal if and only if it is self-adjoint to show that is an orthogonal projection. By induction, because is an orthogonal projection, , and so
[TABLE]
Finally, let be in the range of and we can decompose into the components and such that , and . We will show that , which characterizes the range of as being all vectors orthogonal to . To show this note that because is a projection matrix, we have that
[TABLE]
By construction and so . Using the induction hypothesis, we then have that . Moreover, because by construction, . Then, using the recursive definition of , we have that
[TABLE]
Therefore, and, by Eq. 26, . We now decompose into and where and . By the induction hypothesis, . Therefore, and such that . Finally, using the recursive formulation of and ,
[TABLE]
Thus, we have shown that the range of is orthogonal to . ∎
From Theorem 1, we see that our procedure is an orthogonalization procedure just like quasi-Newton methods (Nocedal and Wright, 2006, Ch. 8) and conjugated direction methods (Hestenes, 2012). As a consequence, we have the following common and insightful characterization of the iterates of such an orthogonalization procedure.
Corollary 1**.**
In addition to the setting of Theorem 1, let be arbitrary and let be defined according to Eq. 14. For any , .
Proof.
We again proceed by induction. Because , the case of follows by recursion formula, Eq. 14. Now suppose that the result holds up to some . Note, by the recursion formula
[TABLE]
Therefore, . Now, using the induction hypothesis,
[TABLE]
Second, when , then . Consequently,
[TABLE]
Now suppose . By Theorem 1, is an orthogonal projection onto . Hence, , which is contained in ∎
Corollary 1 demonstrates that, as is common with orthogonalization procedures, the iterates are in a subspace generated by the initial iterate and the search directions . For deterministic procedures, such a characterization is usually sufficient and the next step would be to demonstrate that the iterates are the closest points to the true solutions within the given subspace. However, for a procedure in which the subspace is randomly generated, there is substantially more nuance. In order to be conscientious of space, we will not go through the litany of issues, but rather skip to the appropriate definitions and characterizations.
First, we begin by defining the maximal possible subspace that can be generated by a random quantity . Let be a random variable defined on a space , and let
[TABLE]
Moreover, we define the subspace such that and (hence, ). Correspondingly, let denote the orthogonal projection matrix onto a subspace . The following result characterizes .
Lemma 2**.**
For as defined in Eq. 32, is the smallest subspace of such that .
Proof.
First, we verify that . Suppose that . Then,
[TABLE]
However, we know that for any such that , and with probability one, which is a contradiction. Hence, .
Now suppose there is a proper subspace of , , such that . Let denote the subspace orthogonal to relative to . Then, for any , which implies that . However, since , . Thus, is the smallest subspace such that . ∎
Second, we must define when the maximal possible subspace of can be achieved by a sequence of random variables , which may or may not be related to . Note, by not requiring a relationship between and our next result is particularly general and applies to a variety of situations, from the case in which are independent copies of to the case where have complex dependencies. Now, let be random variables defined on , and let be a stopping time defined by
[TABLE]
Using this notation, we have the following fundamental characterization result of the limit points of .
Theorem 2**.**
Let be a random variable, and let , and be as defined above (see Eq. 32). Moreover, let be random variables such that for all , and let be as defined in Eq. 34. Let be arbitrary and , and let and be defined as in Eqs. 14 and 13. On the event ,
For any , and . 2. 2.
If admits a solution (not necessarily unique), then
[TABLE]
Proof.
Recall that Therefore, by the definition of , on the event that . Therefore, by Theorem 1, is an orthogonal projection onto and its null space is .
We now proceed by induction. Because and with probability one (by hypothesis), . Therefore, by the recursion equations, Eqs. 14 and 13, and . Suppose now that and for . Again, by hypothesis, . Therefore, . By the recursion equations, Eqs. 14 and 13, and .
To establish the second part of the result, we must first establish that for any ,
[TABLE]
We will prove this by induction. For ,
[TABLE]
by the recursion equations, Eq. 14. Noting that and by using Eq. 13, we conclude that . Now suppose that this relationship holds for some . Again, using Eq. 14,
[TABLE]
Using the induction hypothesis, and Eq. 13,
[TABLE]
With this result established and noting that is a projection onto (i.e., ), on the event ,
[TABLE]
∎
With Theorem 2 in hand, the natural subsequent question is when the limit point of the iterates is actually a solution to the original system. This question is addressed in the following corollary.
Corollary 2**.**
Under the setting of Theorem 2, on the event , if and only if .
Proof.
Recall that . Because , . Moreover, by the definition of , . Therefore, . Now, using the characterization in Theorem 2,
[TABLE]
Similarly, because ,
[TABLE]
Setting these two quantities equal to each other, we conclude that if and only if . Clearly, if then . So, what we have left to show is that implies .
Let denote the Moore-Penrose pseudo-inverse of , and recall that is a projection onto . Moreover, . Therefore, since if then , if then
[TABLE]
∎
Corollary 2 provides criteria on the initial condition and on to determine when our procedure will solve the linear system. However, we would rarely have a way of choosing the initial condition apriori such that the requirement of Corollary 2 holds. Thus, the alternative is to design and so that , which would guarantee that on the event . It is worth reiterating that we have made very limited assumptions about the relationships between and and amongst . This is important because it allows us to apply the preceding results to a variety of common relationship patterns between and . In the next subsection, we explore some specific relationships and whether these relationships will result in .
3.2 Common Sampling Patterns
Theorem 2 supplies a general result about the behavior of any sampling methodology on the solution of the system using Eqs. 13 and 14, yet it does not suggest a precise sampling methodology. Generally, the sampling methodology choice will depend on both the hardware environment and the nature of the problem. For example, a random permutation sampling methodology will limit the parallelism achievable in Algorithm 5. On the other hand, a random permutation sampling methodology might be well-advised in a sequential setting where very little known is about the coefficient matrix . Thus, the precise sampling scheme should depend on the hardware environment and should exploit the structure of the problem.
Despite this, in practice, there are two general sampling schemes that form a basis for more problem and hardware specific sampling schemes: random permutation sampling and independent and identically distributed sampling. The former sampling pattern is exemplified by randomly permuting the equations of the linear system. More concretely, let be the standard basis; let be a random variable with nonzero probability on each element of the basis; let be random variables sampled from without replacement (until the set is exhausted, then we repopulate the set with its original elements and repeat the sampling without replacement). The following statement provides a simple characterization of this sampling scheme.
Lemma 3**.**
Let . Let be a random variable such that
[TABLE]
Moreover, let be random variables sampled from without replacement (and once the set is exhausted, we repopulate the set with its original elements and repeat sampling without replacement). Then . Moreover, for every initialization if , which holds if .
Proof.
First, note that . Therefore,
[TABLE]
In turn, because , is at most .
By Corollary 2, if and only if where satisfies . Now, given that and , if then . Therefore, for any initialization. The final claim is straightforward. ∎
The second sampling scheme, independent and identically distributed sampling, is exemplified by randomly sampling equations from the system with uniform discrete probability. However, we do not need to limit ourselves to sampling from a finite population of elements. As the next result shows, we can do much more.
Proposition 1**.**
Suppose that are independent, identically distributed random variables. There exists a such that
[TABLE]
Moreover, and where and .
Proof.
First, we show that there exists such that for any nontrivial, proper subspace , , which implies Eq. 46 when we take to be the relative orthogonal compliment to the span of a unit vector . Suppose there is no such . Then, for every , there is a nontrivial subspace such that . Let be the smallest integer between [math] and such that
[TABLE]
For , let be an -dimension subspace with . Note, by Lemma 2, . Therefore, let be an -dimensional subspace with . Given that and are distinct and the inclusion-exclusion principle,
[TABLE]
However, this is contradicts the minimality of since is arbitrary and . Thus, we conclude that such a exists.
It follow from Eq. 46 that for any k,
[TABLE]
Therefore, we can bound by a negative binomial distribution. In particular,
[TABLE]
∎
In light of the two preceding results, we may be convinced that there is a gap between the convergence properties between random permutation sampling and the independent and identically distributed sampling. However, by modifying the structure of the rank-one RPM, we can find more intermediate cases. The next result demonstrates this behavior with a somewhat contrived example, and we will leave more complex cases to future work.
Theorem 3**.**
Suppose are i.i.d. random variables such that the entries of are independent, identically distributed subgaussian random variables with mean zero and unit variance. Then, there exists a depending only on the distribution of the entries of such that for .
Proof.
Let denote a () random matrix whose entries are independent and identically distributed subgaussian random variables with zero mean and unit variance. As a consequence of (Rudelson and Vershynin, 2009, Theorem 1.1), there exists a that depends on the distribution of the entries such that for all , . At iteration , let denote the matrix whose rows are given by . Then, by hypothesis, has entries that are independent, identically distributed subgaussian random with zero mean and unit variance. Therefore, there exists a depending only on the distribution of the entries in such that for . ∎
4 Convergence Theory for Base Methods
In the previous section, we proved convergence for the complete orthogonalization method (i.e., Algorithm 1) and explored some specific sampling patterns. Here, we will consider the extreme opposite of the complete orthogonalization method: the “base” randomized iterative approach (e.g., Randomized Kaczmarz). That is, we consider when is a rank one matrix of one of two general classes.
In the first class, we consider Algorithm 2 in the case . In this case, Eq. 14 supplies the simplified iteration scheme,
[TABLE]
which encompasses randomized Kaczmarz, when is a random draw from the standard basis vectors in , as shown in Subsection 2.1.
Unfortunately, Eq. 51 would not include randomized Gauss-Seidel. This motivates the second class, which has the closely related iteration
[TABLE]
In this class, we recover randomized Gauss-Seidel if we choose randomly from the standard basis vectors in , as shown in Subsection 2.1.
While these two classes are distinct, we will see that their analysis is nearly identical and is intimately related to the analysis of the complete orthogonalization method. Our analysis offers two highlights: (1) we can prove convergence with probability one for arbitrary sampling schemes—only the i.i.d. case is considered in Zouzias and Freris (2013); Gower and Richtárik (2015); Richtárik and Takác (2020); and (2) we can provide rates of convergence with probability one which complements the mean-squared-error results of Zouzias and Freris (2013); Gower and Richtárik (2015); Richtárik and Takác (2020). Our main approach is an extension of Meany’s inequality (see Subsection 4.1) combined with stopping time arguments, as derived in Subsections 4.2 and 4.3. We then explore some common, non-adaptive sampling patterns in Subsection 4.4. To conclude, we develop a general framework for the analysis of adaptive sampling schemes, and provide concrete examples from the literature (see Subsection 4.5).
4.1 An Extension of Meany’s Inequality
Here, we will derive an extension of Meany’s Inequality Meany (1969), which, under a different extension, has recently been used to study the convergence rate of row-action solvers including the a block-variant of the Kaczmarz method Bai and Liu (2013). We begin by stating a geometric lemma derived by Meany (1969), and follow it with the extension, which closely follows Meany’s original proof with several modifications.
Lemma 4** (Meany (1969)).**
Let with . Write where belongs to the space spanned by and is perpendicular to . Let be the matrix whose columns are , and let be the matrix whose columns are . Then,
[TABLE]
Theorem 4**.**
Let be unit vectors in for some . Let . Let denote all matrices where the columns of are the vectors that are a maximal linearly independent subset. Then
[TABLE]
where
[TABLE]
Proof.
The proof proceeds by induction. For the case , both sides of the inequality are zero and so the result holds. Now suppose that the result holds for . To prove the case , we need the following additional notation.
Let ; let denote a maximal linearly independent subset of the unit vectors that achieve the minimum determinant; let be the matrix whose columns are ; and let
[TABLE]
For a unit vector , let denote the component of in , and let denote the component of orthogonal to . Moreover, let . Then, by the induction hypothesis,
[TABLE]
Similarly, write where and is perpendicular to .
Case A: Suppose that . Then . Moreover, since ,
[TABLE]
Thus, the result holds when .
Case B: Suppose that . Then,
[TABLE]
where we have made use of and are colinear, implying that their inner product is equal to the product of their norms. Finally, since ,
[TABLE]
Case B(1): Suppose that . Then,
[TABLE]
where, in the last line, we use Lemma 4 and, since , , which, in turn, implies .
Case B(2): Suppose that . Since , then . Using these inequalities and Eq. 57,
[TABLE]
Therefore,
[TABLE]
Applying this relationship to Eq. 62,
[TABLE]
where, in the last line, we use Lemma 4 and, since , , which, in turn, implies .
Therefore, from Cases A, B(1) and B(2), we conclude that the result holds. ∎
4.2 Main Convergence Result for Row-Action Methods
Recall that is a random variable and is a sequence of random variables taking value in chosen such that .777Again, we can avoid this requirement and consider set inclusions below. However, this generalization will require additional, cumbersome notation and there is no practical reason for considering this case. We will now define a sequence of stopping times where ,
[TABLE]
and, if , we define
[TABLE]
else . As an aside, it is worthwhile to note the commonalities between the definition of and the stopping time from Eq. 34.
Moreover, whenever the stopping times are finite, we will define the collection, , for , that contains all matrices whose columns are maximal linearly independent subsets of
[TABLE]
Moreover, define
[TABLE]
Note, it follows by Hadamard’s inequality that .
Theorem 5**.**
Suppose admits a solution (not necessarily unique). Let be a random variable valued in , and let and be defined as above (see Eq. 32). Moreover, let be random variables such that for all . Let be arbitrary and let be defined as in Eq. 51. Then, for any , on the event ,
[TABLE]
where are defined in Eq. 69 and . Therefore, for any ,
[TABLE]
where ; and where we are on the event .
Proof.
From the basic iteration stated in Eq. 51, we have
[TABLE]
Iterating on this relationship, we conclude
[TABLE]
Moreover, by assumption, with probability one, which implies that . Therefore,
[TABLE]
and .
Note, when is finite, then the span of is . Therefore, on the event , Theorem 4 implies that
[TABLE]
We now proceed by induction. Suppose Eq. 70 holds for some . Using Eq. 74, for ,
[TABLE]
Now, when , the conditions of Theorem 4 are satisfied. Therefore,
[TABLE]
By applying the induction hypothesis, we conclude that Eq. 70 holds on the event .
Now, for an orthogonal projection matrix, , . The bound on follows by applying this fact and the definition of . ∎
As an analogue of Corollary 2, we have the following characterization of whether solves the system .
Corollary 3**.**
Under the setting of Theorem 5, on the events and , if and only if .
Proof.
By Theorem 5, and on the events and ,
[TABLE]
Therefore, , which implies if and only if . Clearly, if , then . Now, since , if , then follows from Eq. 43. ∎
4.3 Main Convergence Result for Column-Action Methods
For the family of methods specified by Eq. 52, we will follow an almost identical proof except on the residual rather than the error. Specifically, if we let , then Eq. 52 implies
[TABLE]
Thus, we will see two changes in the proof. First, we will see that see that for column-action methods will take the place of for row-action methods. Second, we already see that in Eq. 79 has taken the place of in Eq. 72. Owing to this latter issue, we will need to specify analogues of , and .
Let be a random variable, and let
[TABLE]
Just as generalized the null space of under the action of an -dimensional random variable from the left, we see that is a generalization of the left null space of under the action of a -dimensional random variable from the right. Analogously, just as restricted the row space of under the action of an -dimensional random variable from the left, we see that is a restriction of the column space of under the action of a -dimensional random variable from the right. Finally, we let denote the subspace that is orthogonal to such that is the column space of .
With these new definitions, we may proceed just as we do in Subsection 4.2. For a random variable , let be a sequence of random variables in such that . We will now define a sequence of stopping times where ,
[TABLE]
and, if , we define
[TABLE]
else .
Moreover, whenever the stopping times are finite, we will define a collection, for , that contains all matrices whose columns are maximal linearly independent subsets of
[TABLE]
We can then define just as we do in Eq. 69. For completeness, we will define it again here so that we reference the appropriate definitions. Define
[TABLE]
Theorem 6**.**
Suppose admits a solution (not necessarily unique). Let be a random variable valued in , and let and be defined as above (see Eq. 80). Moreover, let be random variables such that for all . Let be arbitrary, let be defined as in Eq. 52, and define for . Then, for any , on the event ,
[TABLE]
where are defined in Eq. 84 and . Therefore, for any ,
[TABLE]
where ; and where we are on the event .
Proof.
Iterating on Eq. 79, we conclude
[TABLE]
Moreover, by assumption, with probability one, which implies . Therefore,
[TABLE]
and
Note, when is finite, then the span of is . Therefore, on the event , Theorem 4 implies that
[TABLE]
We now proceed by induction. Suppose Eq. 85 holds for some . Using Eq. 88, for ,
[TABLE]
Now, when , the conditions of Theorem 4 are satisfied. Therefore,
[TABLE]
By applying the induction hypothesis, we conclude that Eq. 85 holds on the event . The second part of the result follows readily. ∎
We have the following characterization of whether solves the system .
Corollary 4**.**
Under the setting of Theorem 6, on the events , and , if and only if .
Proof.
On the events and , Theorem 6 implies
[TABLE]
It straightforwardly follows that if and only if .
Moreover, by construction of , we have that . Thus,
[TABLE]
Since the left null space of is orthogonal to the column space of , and is in the column space of because is consistent, we have that . ∎
4.4 Common, Non-Adaptive Sampling Patterns
Just as for Theorem 2, Theorems 5 and 6 are general results that characterizes convergence for any sampling scheme. Following the discussion in Subsection 3.2, the sampling scheme should depend on the hardware environment and the problem setting. Despite this, the two sampling patterns studied in Subsection 3.2 form a foundation for most sampling schemes in practice and warrant a precise analysis. After this analysis, certain adaptive schemes have become popular and are also analyzed in a generic manner. We will focus on the case of row-action methods (corresponding to Theorem 5) as the column-action results (corresponding to Theorem 6) are nearly identical.
The first result provides a proof of convergence when we sample without replacement from a finite population. We note that the result is quite general and does not depend on the nature of the sampling without replacement or the dependency of the samples whenever the finite population is exhausted. As a result, the bounds are loose, which may be unsatisfying. Should particular sampling patterns become sufficiently important to warrant a more detailed analysis, we will do so in future work.
Proposition 2**.**
Let and be defined as in Lemma 3. Then, under the setting of Theorem 5,
* for all , and* 2. 2.
.
Moreover, are uniformly bounded by that depends on . Therefore, with probability one,
[TABLE]
Proof.
By the definition of in Lemma 3, . Moreover, by the definitions of , we are sampling from without replacement. Then, we are guaranteed that spans if . Now, suppose that at iteration , are exhausted. Then, to ensure that is contained in , we need to exhaust and then the entire set . Since , we need at most more iterations from to achieve . Therefore, . Now, let denote all matrices whose columns are maximal linearly independent subsets of
[TABLE]
Then, . Therefore,
[TABLE]
It is clear, by Hadamard’s inequality, that . Hence, . The result follows by Theorem 5. ∎
It is worth pausing here to compare our approach in Proposition 2 to previous results for cyclic row-action methods (e.g., Kaczmarz (1993),888This is a translated copy of Kaczmarz’s original article, which is published in German (Karczmarz, 1937). algebraic reconstruction technique (Gordon et al., 1970), cyclic block Kaczmarz). Our use of Meany’s inequality to analyze such methods is not novel: Meany’s inequality has been used previously to analyze deterministic row-action methods (Galántai, 2005; Bai and Liu, 2013; Wallace and Sekmen, 2014) with even more sophisticated refinements of Meany’s inequality than what we have here, and a detailed comparison of Meany’s inequality and other approaches to analyzing these deterministic variants can be found in Dai and Schön (2015). However, our use of Meany’s inequality generalizes these deterministic approaches as it (1) allows for an arbitrary transformation (via ) of the original system, which has borne out to be a fruitful approach vis-à-vis matrix sketching Woodruff (2014); and (2) allows for the benefits of random cyclic sampling, which many have observed to be the most productive route in practice and there is mounting evidence in adjacent fields that random cyclic sampling does indeed have practical benefits (Lee and Wright, 2019; Wright and Lee, 2020). While our generalizations are valuable, further improvements are to be found by marrying our randomization framework with the more nuanced refinements of Meany’s inequality found in Galántai (2005) and Bai and Liu (2013), which we leave to future efforts.
The next result revisits the case of independent and identically distributed sampling. The result makes intuitive sense as, for such a situation, we should expect the difference in the stopping times to be independent and identically distributed, which, results in the natural conclusion that are also independent and identically distributed. Moreover, we show that eventually, the rate of convergence is almost controlled by with probability one. We again stress here that the generality of the results naturally makes them quite loose, and we discuss this further after the result.
Proposition 3**.**
Let and be defined as in Proposition 1. Then, under the setting of Theorem 5, almost surely for all , and are independent and identically distributed such that . Hence, for all and ,
[TABLE]
where . Moreover, .
Remark 2**.**
In the proof below, we also compute the probability for each for which the conclusion of the preceding result holds. Thus, we can also make the usual “high-probability” statements without any additional effort.
Proof.
Again, our main workhorse will be (Durrett, 2010, Theorem 4.1.3). By this result, conditioned on , are independent and identically distributed. By this property, conditioned on , is independent of and have the same distribution for all . We conclude then that since is a function of , then are independent and identically distributed. We now conclude that Eq. 70 holds with probability one by applying Theorem 5. For any , by Markov’s inequality and independence,
[TABLE]
Since , the Borel-Cantelli lemma implies that the probability that the product of is eventually less than is one. ∎
Here, we again take a moment to compare this result to the results of Richtárik and Takác (2020). Namely, we are interested in how the rate of convergence of Proposition 3 compares with the rate of convergence result in Richtárik and Takác (2020). To make this comparison, we numerically estimate the theoretical rates of convergence proposed by our result and the result of Richtárik and Takác (2020) on five matrices from the MatrixDepot (as described in Section 5). We show these comparisons in Table 2. We show these comparisons in Table 2. As expected, the results of Richtárik and Takác (2020), which are specialized to the i.i.d. case and apply on average, are much tighter than our general results that apply to more than just i.i.d. case and hold with probability one.
4.5 Adaptive Sampling Schemes
To bookend this section, we discuss how our results can be applied to a broad set of adaptive methods that make use of the residual information at a given iterate whether deterministically (e.g., Motzkin and Schoenberg (1954); Gubin et al. (1967); Lent (1976); Censor (1981)) or randomly (e.g., Nutini et al. (2016); Bai and Wu (2018); Haddock and Ma (2019)). In Subsection 4.5.1, we will begin with some formalism to establish a general class of adaptive methods, and we then prove convergence and a rate of convergence for such methods. In Subsection 4.5.2, we provide concrete examples at the end.
4.5.1 A General Class and Analysis of Adaptive Methods
To be rigorous, let and let be an adaptive procedure for generating according to the following procedure: for ,
[TABLE]
Remark 3**.**
While we will focus on the base methods of type Eq. 51, methods of the type Eq. 52 can be handled analogously.
While Eq. 99 is quite general, the vast majority of adaptive schemes make further restrictions that we abstract in the following definitions.
Definition 1** (Markovian).**
For a fixed integer , an adaptive procedure, , is -Markovian if the conditional distribution of given is equal to the conditional distribution of given . If a procedure is -Markovian, we will frequently call it Markovian.
A consequence of the -Markovian property is that we can write as . In the case of a -Markovian adaptive procedure, we will simply write . The -Markovian property is readily satisfied for a number of common procedures analyzed in the literature (e.g., maximum residual, maximum distance, etc.), which may suggest that the -Markovian notion is irrelevant for general . We contend though, that procedures that are memory-sensitive may be more apt to make use of the -Markovian property for . For example, to demonstrate its potential value, consider a procedure that selects the equations with the top residuals, pulls them into memory, and simply cycles through them deterministically or randomly. Then this simple procedure would be -Markovian. However, owing to the lack of such procedures in the literature, we will focus on the -Markovian case for which we can write , and note that the results and definitions are readily extendable.
The next definition establishes another key property of these adaptive schemes that rely on residuals.
Definition 2** (Magnitude Invariance).**
Let represent the set of solutions to , and let represent the projection of a vector onto ,999Since is a flat, is not guaranteed to be a linear operator. then an adaptive procedure, , is magnitude invariant if, for any and any , the distribution of is equal to the distribution of
[TABLE]
The magnitude invariance of a number of adaptive methods often follows from the following simple calculation that we state as a lemma for future reference.
Lemma 5**.**
Let and let . Then, for any , if then
[TABLE]
If the hypothesis holds with a strict inequality, then so does the conclusion.
Proof.
Note, . Therefore, . From the hypothesis and , . Also owing to , we can replace the inequalities with strict inequalities. ∎
Furthermore, the magnitude invariance property has hidden within it an additional feature: the projection of onto the null space is irrelevant (as we might expect for a procedure depending on the residual). As a result, we can, without losing generality, focus our discussion to that are in the row space of , which has a unique intersection with at a point that we denote . Furthermore, the magnitude invariance property allows us to focus specifically on the Euclidean unit sphere around , which we denote by . This will be essential to the next definition.
The final definition ensures that if Eq. 99 makes too much progress along one particular subspace, then it must have a nonzero probability of exploring an orthogonal subspace relative to, roughly, the row space of . Before stating this definition, we need to be slightly careful here with using the row space of : if the rows of can be partitioned into two sets that are mutually orthogonal and is initialized in the span of one of these subsets, then we will never need to visit the other set and, consequently, we will never observe the entire row space of . To account for this, we can focus on the restricted row space,
[TABLE]
This definition may seem unnecessary as we can account for this (more generally) via by an appropriate choice of . However, in our previous statements, we defined before specifying . Here, we would need to know in order to define and, thus, appropriately. Fortunately, an examination of the preceding results shows that this ordering is not important and the results hold even if is defined given or even future iterates. With this explanation in hand, we can now state the final definition.
Definition 3** (Exploratory).**
Let and define accordingly. An adaptive procedure, , is exploratory if for any proper subspace , there exists such that
[TABLE]
Remark 4**.**
If magnitude invariance does not hold, then we could specify the exploratory property to hold for any point in that is distinct from . For this modified definition of the exploratory property, the results below would still hold. Then, why should we keep the magnitude invariance property? It is out of practicality. The magnitude invariance property allows us to restrict the verification of the exploratory property to the unit ball, and then we can apply it to any iterate regardless of its distance to the solution.
For a Markovian, magnitude invariant and exploratory adaptive scheme, , we will need one assumption before stating the result.
Assumption 1**.**
Let denote the set of matrices whose columns are normalized, maximal linearly independent subsets of
[TABLE]
where are arbitrary vectors. Suppose, for this choice of ,
[TABLE]
Remark 5**.**
As we will see, Assumption 1 is sufficient for us to uniformly treat the many examples in the literature that are selecting equations or, more generally, are of the form in Lemma 3, rather than generating linear combinations of them. In the case of linear combinations, we could refine this assumption to account for the nature of the linear combinations as we do in Proposition 3.
Theorem 7**.**
Suppose admits a solution (not necessarily unique); let denote the set of all solution, and be the projection onto this flat. Let and let be defined as above (see Eq. 102). Moreover, let be a -Markovian, magnitude invariant and exploratory adaptive procedure satisfying Assumption 1 that generates and according to Eq. 99 and so that for all . Then, there exist an increasing sequence of stopping times such that , where:
* is the event of iterates that terminate finitely to a solution of ; that is,*
[TABLE] 2. 2.
* is the event of iterates that infinitely converge to a solution of ; that is,*
[TABLE]
Moreover, on , has finite expectation for such that . Similarly, on , has finite expectation for all .
Proof.
Without loss of generality, we will assume . We will consider the nontrivial case where . Note, by the construction of , it must hold then . To prove the result, we will make three claims of the following rough nature and purpose, which we will make precise below.
Finite termination can only occur at a point if and only if is parallel to . We will use this claim to specify the set . 2. 2.
For the first time the span of the iterate errors, , fails to (non-trivially) increase in dimension, the corresponding up to this iterate span the subspace. As a result, with an appropriate definition of , we will apply Theorem 5 to prove a multiplicative decrease in the iterate errors by a factor of . 3. 3.
Finally, we show that the first time that the span of the iterate errors fails to (non-trivially) increase in dimension must be finite with probability one and have bounded expectation. By combining the first claim with this claim, we have the property specified by the event . By combining this claim with the second claim, we have the property specified by the event . By this claim alone, we have that .
To establish our claims, we need some additional notation. Let be an arbitrary finite stopping time and define
[TABLE]
and . Furthermore, define
[TABLE]
Note, corresponds to the first time that the span of the iterate errors, starting at , fails to non-trivially increase in dimension. It will often be more succinct to specify the non-trivial cases by an indicator variable given by
[TABLE]
By Eq. 99, we can readily replace in the definition of with . We now state and prove our claims precisely.
Claim 1: Suppose . We claim that if and only if .
Note, this claim readily follows from
[TABLE]
which, in turn, follows from Eq. 99.
Claim 2: Suppose is finite and define . We claim that
[TABLE]
We first note that for any by Eq. 99. Therefore, we see that the span of is contained in . To show that is included in the span of , note that, by the definition of and by Eq. 99,
[TABLE]
Moreover, the nonzero terms on the generating set on the right hand side of Eq. 113 must be linearly independent, as anything else would contradict the minimality of . We are left to show that is in the span of . To do this, we perform Gram-Schmidt on the generating set in Eq. 113 starting with . Denote the remaining vectors in this set where . Then, by the definition of , . Therefore, there exist constants such that
[TABLE]
If , we see that the claim follows. For a contradiction, suppose that . Then can be written as a linear combination of vectors that are orthogonal to . This would imply then that , which contradicts the definition of . Hence, we see that the claim holds.
Claim 3: For any finite stopping time , is finite with probability one and has bounded expectation.
To show this, we define a sequence of stopping times. Define
[TABLE]
and
[TABLE]
By the definition of , can only take values in . Moreover, at each , we must either observe or . Hence, at most, we see that can only take values in where . Thus, if we show that each is finite and has bounded expectation, then must be finite and have bounded expectation. By the magnitude invariance, Markovian and exploratory properties, we conclude that
[TABLE]
Therefore, we see that is finite and has bounded expectation.
Conclusion: From these three claims we can now prove the result by induction.
Base Case
Define . On this event, if we take and define to be the corresponding . On , is finite and has finite expectation by Claim 3. Then, we can define, as a subset of ,
[TABLE]
and to be its relative complement on .
Note,
By Claim 1, is equivalent to the event up to a measure zero set. 2. 2.
By Claim 2, Theorem 5 with , and Assumption 1, is contained in the event on which
[TABLE]
up to a measure zero set.
Induction Hypothesis
Let . On the event , we let and, for the correspondingly defined , we can define . Furthermore, on , is finite and has finite expectation. We can define, as a subset of ,
[TABLE]
and to be its relative complement on .
Further,
is equivalent to the event up to a measure zero set. 2. 2.
is contained in the event on which
[TABLE]
up to a measure zero set.
Generalization
On the event , we let and, for the correspondingly defined , we can define . On , is finite and has finite expectation by Claim 3. Then, we can define, as a subset of ,
[TABLE]
and to be its relative complement on .
By Claim 1, is equivalent to the event up to a measure zero set. 2. 2.
By Claim 2, Theorem 5 with , and Assumption 1, is contained in the event on which
[TABLE]
up to a measure zero set.
Therefore, by the induction claims,
[TABLE]
and
[TABLE]
and . ∎
4.5.2 Applying our General Theory to Specific Adaptive Schemes
To demonstrate the utility of Theorem 7, we show that a number of classical and recent methods satisfy Definitions 2, 1, 3 and 1. In fact, we will show that a stronger version of Definition 3 holds for these methods, which allows us to explicitly upper bound the elements of (when they are defined).
Proposition 4**.**
Suppose admits a solution . Let and let be defined as above (see Eq. 102). Suppose that we define and according to Eq. 99 for the following adaptive methods
the maximum residual method (see Agmon, 1954, Section 4)); 2. 2.
the maximum distance method (see Agmon, 1954, Section 3); 3. 3.
the Greedy Randomized Kaczmarz method (see Bai and Wu, 2018, Method 2); 4. 4.
the Sampling Kaczmarz-Motzkin method (see Haddock and Ma, 2019, Page 4).
Then, for each of the above methods, there exists a such that the conclusions of Theorem 7 hold. Moreover, there exists a constant such that for any finite (as specified in Theorem 7), .
Remark 6**.**
Greedy Randomized Kaczmarz is an example of methods that deterministically determine a threshold over residuals; select the equations whose residuals surpass this threshold; and then randomly select from this set. For this more general class, so long as the threshold satisfies the magnitude invariance property and the random selection does not give any equation less than zero probability, then the result applies to this more general class. Similarly, Sampling Kaczmarz-Motzkin is an example of methods that randomly determine a set of equations; and then deterministically select from this subset of equations based on the residual values. So long as the random subset of equations does not give any equation less than zero probability (that is not already satisfied), then the result will apply to this more general class as well.
Remark 7**.**
Our partial orthogonalization methods (see Algorithm 2) do not satisfy the -Markovian property, as the partial orthogonalizations have a dependence on every preceding iterate.
For each method, we show that it satisfies Definitions 2, 1, 3 and 1. In fact, for each method, we will show that a stronger version of Definition 3 holds. We will start by establishing several general facts that will be useful in the discussion of each method.
Lemma 6**.**
Let and define as in Eq. 102. Then,
[TABLE]
where is the Euclidean unit sphere around the zero vector.
Proof.
For each , we see that
[TABLE]
else and and we would have a contradiction since . By continuity, we see that we can construct an open ball around each , , such that
[TABLE]
for all . Now, is an open cover of , which is a compact space. Hence, there is a finite subcover given by . It follows that since each belongs to one of the elements in the subcover, then
[TABLE]
Therefore . ∎
Lemma 7**.**
Let and define as in Eq. 102. Let . Let be the matrices whose columns are normalized, maximal linearly independent vectors from . Then
[TABLE]
Proof.
There are only a finite number of matrices in up to column permutations. Therefore, we can choose the that minimizes . By Hadarmard’s inequality, , which implies that . ∎
Maximum Residual Method.
In the maximum residual method, is the standard basis vector in , , that solves
[TABLE]
-Markovian: It follows from the definition of the maximum residual method that it only relies on the current iterate to evaluate . Therefore, it is -Markovian.
Magnitude Invariance: By Lemma 5, it follows that is magnitude invariant.
Exploratory: Consider any . Then, . Therefore, we have that the only equations whose residuals are non-zero are the ones such that , and there is at least one such equation by Lemma 6. Therefore,
[TABLE]
That is, we satisfy the exploratory property in a stronger manner:
[TABLE]
With these three properties verified and by Lemma 7, the conditions of Theorem 7 are satisfied and the result holds. The only thing left to show is that are bounded by some . By the proof of Theorem 7, it is enough to bound the conditional expectations of in Eq. 117. Given that for all ,
[TABLE]
Hence, . Thus, .
Maximum Distance Method.
In the maximum distance method, is the standard basis vector in that solves
[TABLE]
-Markovian: It follows from the definition of the maximum distance method that it only relies on the current iterate to evaluate . Therefore, it is -Markovian.
Magnitude Invariance: Note, Lemma 5 still holds if we were to divide by the norm squared of . It follows that the maximum distance method is magnitude invariant.
Exploratory: Just as in the maximum residual method, if that is orthogonal to a subspace , then for any . Moreover, by Lemma 6, there is at least one equation such that for all . Hence, the maximum distance method satisfies a stronger version of the exploratory condition, namely,
[TABLE]
By the same argument as above, Theorem 7 follows. Similarly, .
Greedy Randomized Kaczmarz.
In Bai and Wu (2018) (Method 2), a residual threshold is selected given by
[TABLE]
Then, from the set of equations whose residual surpasses this threshold (which is shown to at least contain the equation selected by the maximum distance method), an equation is selected by a probability proportional to the equation’s residual squared.
-Markovian: Given that the threshold relies only on the current iterate value and that the random selection criteria only relies on the current residual, it follows that the Greedy Randomized Kaczmarz method is -Markovian.
Magnitude Invariance: Suppose . For , let . Then, by Lemma 5,
[TABLE]
which implies that the threshold is magnitude invariant. Similarly, we can show that the selection probabilities are magnitude invariant (we look at the preceding calculation, but only for a nonempty subset of the equations).
Exploratory: Let be a nontrivial subspace. Then for any , we saw that any equations for which have a zero residual. Therefore, the only equations with nonzero residuals are those that not orthogonal to . Since the threshold is bounded away from zero, only equations that are not orthogonal to can be in the subset. Therefore,
[TABLE]
By the same argument as above, Theorem 7 follows. Similarly, .
Sampling Kaczmarz-Motzkin.
In Haddock and Ma (2019) (Page 4), a subset of equations are randomly selected, and then the equation with the maximum residual is selected from this subset.
-Markovian: The Sampling Kaczmarz-Motzkin method only relies on the current residual to sample. As a result, it is -Markovian.
Magnitude Invariance: The distribution of the initial subsetting is independent and identical at each iteration. Therefore, conditioned on a given subset, we choose the maximum residual. By Lemma 5, this last step is magnitude invariant. Moreover, since the random subsetting is independent and identical at each iteration, it too is magnitude invariant. Therefore, the entire procedure is magnitude invariant.
Exploratory: Let be a nontrivial subspace. Then, for any , we have shown that there exists a such that . Therefore, so long as the probability of selecting this equation is nonzero, then we are guaranteed that there is some choice of such that
[TABLE]
Let be the smallest inclusion probability for any equation in the random subset. Then, it follows that
[TABLE]
For the Sampling Kaczmarz-Motzkin method, the minimum inclusion probability is at least , which corresponds to random sampling without replacement of subsets of size .
With these three properties verified and by Lemma 7, the conditions of Theorem 7 are satisfied and the result holds. The only thing left to show is that are bounded by some . By the proof of Theorem 7, it is enough to bound the conditional expectations of in Eq. 117. Supposing that for all ,
[TABLE]
Hence, .
5 Numerical Experiments
Here, we present a variety of numerical experiments to study the practicality of our approach in a sequential computing environment. Specifically, we test forty-nine systems with five hundred equations and five hundred unknowns. The coefficients are generated from forty-nine built-in matrices found in the MatrixDepot package for the Julia programming language (Zhang and Higham, 2016). The solution to the equation is then generated using a standard, multi-variate normal vector. The constant vector is generated by the product of the two. Then, using the generated coefficient matrix and the generated constant vector, we solve the systems by varying the sample-generation method (i.e., the generation of and ) and the solver. The sample generation method is either produced by the count-sketch approach, the Gaussian approach, by uniformly sampling the equations of the matrix with replacement, or by uniformly sampling the equations of the matrix without replacement. The solver is either a base method, the complete method, an intermediate method with , or an intermediate method with . Finally, we initialize .
We recorded the wall clock time to improve the initial residual norm by a factor of ten with an upper bound of three seconds (note, a single iteration of a base method requires approximately seconds, which allows the base method on the order of iterations on a system). If the temporal upper bound is reached before a ten fold improvement in the initial residual norm is observed, the wall clock times is reported as . Inherently, this metric results in a disadvantage for complete orthogonalization methods as such methods pay more for marginal improvements, but generate precise solutions with fewer iterations. However, with an eye towards solving much larger systems that require using a parallel or distributed environment, this metric of time-to-ten-fold improvement is the appropriate choice as the complete method would not be appropriate in such environments owing to the high communication costs that would be incurred. For the count-sketch sampling method, the wall clock times are reported in Section 5. For the remaining sampling approaches, the wall clock times are reported in the appendix.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Agmon [1954] Shmuel Agmon. The relaxation method for linear inequalities. Canadian Journal of Mathematics , 6:382–392, 1954.
- 2Bai and Liu [2013] Zhong-Zhi Bai and Xin-Guo Liu. On the meany inequality with applications to convergence analysis of several row-action iteration methods. Numerische Mathematik , 124(2):215–236, 2013.
- 3Bai and Wu [2018] Zhong-Zhi Bai and Wen-Ting Wu. On greedy randomized kaczmarz method for solving large sparse linear systems. SIAM Journal on Scientific Computing , 40(1):A 592–A 606, 2018.
- 4Censor [1981] Yair Censor. Row-action methods for huge and sparse systems and their applications. SIAM review , 23(4):444–466, 1981.
- 5Chen and Powell [2012] Xuemei Chen and Alexander M Powell. Almost sure convergence of the kaczmarz algorithm with random measurements. Journal of Fourier Analysis and Applications , 18(6):1195–1214, 2012.
- 6Clarkson and Woodruff [2017] Kenneth L Clarkson and David P Woodruff. Low-rank approximation and regression in input sparsity time. Journal of the ACM (JACM) , 63(6):1–45, 2017.
- 7Cormode and Muthukrishnan [2005] Graham Cormode and Shan Muthukrishnan. An improved data stream summary: the count-min sketch and its applications. Journal of Algorithms , 55(1):58–75, 2005.
- 8Dai and Schön [2015] Liang Dai and Thomas B Schön. On the exponential convergence of the kaczmarz algorithm. IEEE Signal Processing Letters , 22(10):1571–1574, 2015.