Iterative Row Sampling
Mu Li, Gary L. Miller, Richard Peng
Introduction
Overview
Finding such a is equivalent to reducing the size of a regression problem involving since:
This means finding a shorter approximation to the matrix , and solving a regression problem on this approximation gives a solution within of the minimum.
The main framework of our algorithm is iterative in nature and relies on the two-way connection between row sampling and estimation of sampling probabilities. A crude approximation to , allows us to compute equally crude approximations of sampling probabilities, while such probabilities in turn lead to higher quality approximations. The computation of these sampling probabilities can in turn be sped up using a high quality approximation of . Our algorithm is based on observing that as long as has smaller size, we have made enough progress for an iterative algorithm. A single step in this algorithm consists of computing a small but crude approximation , finding a higher quality approximation to , and using this approximation to find estimates of sampling probabilities of the rows of . This leads to a tail-recursive process that can also be viewed as an iterative one where the calls generate a sequence of gradually shrinking matrices, and sampling probabilities are propagated back up the sequence. An example of such a sequence is given in Figure 1.
We will term the creation of the coarse approximation as reduction, and the computation of the more accurate approximation based on it recovery. As in the figure, we will label the matrices that we generate, as well as their approximations using the indices .
reduction: creates a smaller version of , with fewer rows either by a projection or a coarser row sampling process. Equilvaent to moving rightwards in the diagram.
recovery: finds a small, high quality approximation of , using information obtained from , , and . This is done by estimating leverage scores and is equivalent to moving leftwards in the diagram.
Preliminaries
This sampling process can be formalized in several ways, leading to similar results both theoretically and experimentally [IW12]. We will treat it as a blackbox that takes a set of probabilities over the rows of and samples them accordingly. It keeps row or with probability , and rescales it appropriately so the expected value of this row is preserved. The two key properties of that we will use repeatedly are:
It returns with at most rows.
Its running time can be bounded by .
The convergence of sampling relies on matrix Chernoff bounds, which can be viewed as generalizations of single variate concentration bounds. Necessary conditions on the probabilities can be formalized in several ways, with the most common being statistical leverage scores. Although these values have been studied in statistics, their use in algorithms is more recent. To our knowledge, their first use in a limited row sampling setting was in spectral sparsification of graphs [SS08]. The most general definition of -norm leverage scores is based on the row norms of a basis of the column space of . However, significant simpliciations are possible when , and this alternate view is crucial for in our algorithm. As a result, we will state the relevant convergence results for Sample separately in Sections 4 and 5.
They show that statistical leverage scores are closely associated the probabilities needed for row sampling, and give algorithms that efficiently approximate these values. We will also formalize an observation implicit in previous results that both the sampling and estimation algorithms are very robust. The high error-tolerance of these algorithms makes them ideal as core routines to build iterative algorithms upon.
where is the -th row of .
To our knowledge, the first near tight bounds for row sampling using statistical leverage scores were given in [AW02], and various extensions and simplifications were made since [RV07, Ver09, Har11, AT11]. They can be stated as follows:
The importance of statistical leverage scores can be reflected in the following fact, which implies that we can obtain with rows.
(see e.g. [SS08]) Given matrix , and let be the leverage score w.r.t. . Assume has rank , then
Although it is tempting to directly obtain high quality approximations of the leverage scores, their computation also requires a high quality approximation of , leading us back to the original problem of row sampling. Our way around this issue relies on the robustness of concentration bounds such as Lemma 4.1. Sampling using even crude estimates on leverage scores can lead to high quality approximations [DDH+09, DM10, DMMS11, DMIMW12, AT11]. Therefore, we will not approximate directly, and instead obtain a sequence of gradually better approximations. The need to compute sampling probabilities using crude approximations leads us to define a generalization of statistical leverage scores.
If and satisfies:
Then for any vector we have:
This representation leads to faster algorithms for estimating stretch using the Johnson-Lindenstrauss transform. This tool is used in a variety of settings from estimating effective resistances [SS08] to more generally leverage scores [DMIMW12]. We will use the following randomized projection theorem:
2 Reductions and Recovery
Our reduction and recovery processes are based on projecting to one with fewer rows, and moving the estimates on the projection back to the original matrix. Our key operation is to combine every rows into rows, where and are set to and respectively. By padding with additional rows of zeros, we may assume that the number of rows is divisible by . We will use to denote the number of blocks, and use the notation to index into the th block. Our key step is then a -reduction of the rows:
A -reduction of describes the following procedure:
For each block , pick to be a random Gaussian matrix with entries picked independently from and compute \mathbf{A}\textrm{\tiny\downarrow}_{(b)}=\mathbf{U}_{(b)}\mathbf{A}_{(b)}.
Concatenate the blocks \mathbf{A}\textrm{\tiny\downarrow}_{(b)} together vertically to form \mathbf{A}\textrm{\tiny\downarrow}.
We first show that projections preserve the stretch of blocks w.r.t. . This can be done by bounding the effect of on the norm of each column of . It follows directly from properties of the Johnson-Lindenstrauss projections described in Lemma 4.5, and we’ll give its proof in Appendix B.
Assume for some constant and let \mathbf{A}\textrm{\tiny\downarrow} be a -projection of . For any constant there exists a constant such that
holds for all block with probability at least .
However, generalized stretches w.r.t. and \mathbf{A}\textrm{\tiny\downarrow} are evaluated under the norms given by the inverses of these matrices, and (\mathbf{A}\textrm{\tiny\downarrow}^{T}\mathbf{A}\textrm{\tiny\downarrow})^{+}. As a result, we need to bound the operator bound between these two pseudoinverses, which we obtain using the following lemma.
Let and be symmetric positive semi-definite matrices and let be the orthogonal projection operator onto the range space of . Then:
This is straightforward when both and are full rank, or share the same null space. However, as pseudo-inverses do not act on the null space, it is crucial that we’re only considering vectors of the form . This Lemma is proven in Appendix B. Combining it with bounds in the other direction allows us to bound the distortion caused by switching reference from to \mathbf{A}\textrm{\tiny\downarrow}.
For any constant , there exists a constant , such that with probability at most , we have for each row of \mathbf{A}\textrm{\tiny\downarrow}, denoted by \mathbf{a}\textrm{\tiny\downarrow}_{i}, satisfies
Proof Denote by the -th block of and \mathbf{A}\textrm{\tiny\downarrow}_{(b)} the corresponding block in \mathbf{A}\textrm{\tiny\downarrow}, by Lemma 4.9,
Since each consists of independent random variables chosen from , is distributed as . This gives:
Further note that \mathbf{a}\textrm{\tiny\downarrow}_{i} is completely contained within the range space of \mathbf{A}\textrm{\tiny\downarrow}^{T}\mathbf{A}\textrm{\tiny\downarrow}. Therefore for all , \mathbf{P}\mathbf{a}\textrm{\tiny\downarrow}_{i}=\mathbf{a}\textrm{\tiny\downarrow}_{i} and:
For any constant , there exists a setting of constants such that for any , we have with probability at least
3 Iterative Algorithm
Assume and satisfy the following condition
Then for any constant , there is a setting of the constants such that
holds with probability at least .
Proof The given condition implies that satisfies the condition needed for Lemma 4.6 with . Let the constants , then with probability at least we have:
Since is a projection of , we can index corresponding blocks in them. Apply Corollary 4.12, then with probability at least , we have
holds for all blocks with probability at least , which is equal to by the definition of .
and has rows, each being a scaled copy of some row of ,
Proof We first show correctness via. induction backwards on . Define , we show that has rows and satisfies
with probability at least for each .
As is a constant, one -projection decreases the number of rows by a factor of . After projections, we get that has rows. Therefore the base case where follows from .
For the inductive step, we assume that the inductive hypothesis holds for and try to show it for . As was set to , we have:
This allows us to invoke Lemma 4.13, which combined with Lemma 4.1 gives that with probability , the inductive hypothesis also holds for .
We will use this Fact with one of or being , in which case it gives:
If ,
If ,
Let be an matrix of rank , and be its dual norm such that . Then an matrix is an -well-conditioned basis for the column space of if the columns of span the column space of and:
.
where is the -th row of and is a constant depending only on . Then with probability at least , returns satisfying
The estimation of and sampling can then be done in a way similar to Section 4. When , we can compute approximations using -stable distributions [Ind06] in a way analogous to Section 4.2.1. of Clarkson et al. [CDMI+12]. When , we will use the -norm as a surrogate at the cost of more rows. As all of our calls to Sample will be using probabilities estimated via. the same matrix, we will the estimation of -norm leverage scores and sampling as a single blackbox.
For any constant , there exist an algorithm that given a , such that is a -well-conditioned basis for , returns a matrix with probability at least such that:
Proof We start by show that is close to the identity matrix as an operator. Also, since is a full rank matrix and is a projection operator onto the column space of , we have Taking pseudoinverses of the given condition on gives:
Substituting it into then gives:
This allows us to infer that , and for any vector , .
Next we find values of and that meet the requirements of a well-conditioned basis given in Definition 5.2. Let be the dual norm for which satisfies .
First consider the case where . We can view all entries of the matrix as a vector of length vector. Apply Fact 5.1 gives:
Which gives and therefore . For the second part, given any vector , we have
Now we consider the case where similarly.
Combining the bounds from these two cases on gives that is a well-conditioned basis, where
It can be checked that in both cases the stated bound on holds.
And the number of rows in can be bounded by:
Proof It suffices to verify both conditions of Definition 5.2 holds for . For , we can treat its th power as a summation over the columns of and get:
Proof of Lemma 5.6: By the guarantees of RowCombineL2 given in Theorem 4.14, we can set its constants so that with probability at least we have:
And the number of rows in can be bounded by
To bound the runtime, we first show that if has rows, then in the next iteration the number of rows in can be bounded by where is a constant based on . When , Lemma 5.6 gives that there exist some constant such that the new row count can be bounded by:
So and suffices. Similarly for the case where , it can be checked that
Gives that the new row count can be bounded by where . As we can set to any arbitrary constant, the constant in front of its exponent can also be removed.
2 Fewer Rows by Iterating Again
A closer look at the proof of Lemma 5.8 shows that a significant increase in the number of rows comes from dividing by the term in the exponent of . As a result, the row count can be further reduced if the leverage scores are computed via. a -norm approximation where is between and . For simplicity we only show this improvement for the case where .
We will start by proving a generalization of Lemma 5.5.
This leads to a result analogous to Lemma 5.6. Solving for the fixed point of this process allows us to prove Theorem 2.2.
Proof of Theorem 2.2: Similar to the proof of Lemma 5.8, in each iteration we reduce the number of rows from to . Ignoring terms in and , we have that the number of rows converges doubly exponentially towards:
This is minimized when , giving . As we can set , the number of rows in can be bounded by for any constant .
References
Appendix A Properties and Estimation of Stretch and Leverage Scores
We now give proofs for estimating leverage scores and row sampling that we stated in Sections 4 and 5.
Proof of Fact 4.4: Note that both stretch and the Frobenius norm acts on the rows independently. Therefore it suffices to prove this when has a single row, aka. . In this case the cyclic property of trace gives:
Proof of Lemma 4.3: The condition given implies that the null spaces of and are identical, giving:
Applying this to the vector gives:
A.2 Estimation of Generalized Stretch
Based on this fact, we can estimate these scores using randomized projections in a way that’s by now standard [SS08, DMIMW12]. Pseudocode of our estimation algorithm is shown in Algorithm 4, while the error analysis is nearly identical to the ones given in Section 4 of [SS08] and Section 3.2. of [DMIMW12].
We remark that can be replaced by any matrix whose product with its transpose equals to . nee candidate for this is , and using it would avoid computing the power of a matrix. However, from a theoretical point of view both of these operations take time, and we omit this extra step for simplicity.
where the last inequality is due to the assumption that . By a suitable choice of constants in this can be made the above probability less than , taking a union bound over the rows gives Part 1.
A.3 Estimation of pp-Norm Leverage Scores
We will estimate the values of using similar dimensionality reduction theorems. Specifically, we utilize a result on -stable distributions first shown by Indyk [Ind06].
Note that the result from [Ind06] was only stated in terms of obtaining approximations, which leads to a factor of on the leading term. However, these bounds can be obtained analogously by the fact that -stable distributions have bounded derivative.
The guarantees on the output then follows from computing Lemma 5.3, while the total running time follows from the cost of evaluating .
Appendix B Deferred Proofs from Section 4
Where the last inequality is by with the assumption that . If we to and substitute , we get:
Proof of Lemma 4.10: Consider an orthonormal basis for the range space of , . Since , this basis can be extended to an orthonormal basis to the range space of by adding . It suffices to prove the claim under this basis system. Here and can be rewritten as by proper rotation:
where and are strictly positive definite. Furthermore, since is positive semi-definite we have that is also positive semi-definite. For any vector , gives a vector that’s non-zero only in the first entries. Let this part be . Then evaluating becomes solving the following system:
The second equation gives . Substituting it into the first one gives:
Note that this is the same as taking the partial Cholesky factorization onto the range space of . Combining things gives:
Since both and are positive definite, we have and therefore:
holds for every .