Randomized Projection for Rank-Revealing Matrix Factorizations and Low-Rank Approximations Thanks:  Revised from Randomized QR with Column Pivoting for submission to SIGEST

Randomized projection allows us to reduce communication complexity at the cost of increasing uncertainty regarding latent 2-norms (and inner products) among columns in a large matrix. Randomized QR with Column Pivoting (RQRCP), shown in algorithm 4, harnesses this trade-off to obtain full blocks of pivot decisions that can be applied to $\boldsymbol{A}$ all at once. We do this by drawing a Gaussian Independent Identically Distributed (GIID) matrix, $\boldsymbol{\Omega}$, and compressing columns of $\boldsymbol{A}$ into a sample matrix, $\boldsymbol{B}=\boldsymbol{\Omega}\boldsymbol{A}$. The sample matrix retains enough information about the original matrix for us to obtain a full block of pivoting decisions, while requiring far less communication to do so. Because we construct the sample to contain far fewer rows than the original, it resides in a faster level of memory than the full matrix. Having a full block of pivots allows us to update the matrix that becomes $\boldsymbol{R}$ with blocked BLAS-3 operations in the same fashion as unpivoted QR, thus entirely eliminating the communication bottleneck of BLAS-2 operations encountered in the standard algorithm.

We show how rank-revealing decompositions are well suited to this approach because each pivot decision must merely avoid selecting a column containing a relatively small component orthogonal to the span of previous pivots. As many columns are often suitable for this purpose, the use of precise 2-norm computations in standard QRCP is unnecessarily. Our approach has been adopted in subsequent work in computing matrix factorization-based low-rank approximations on sequential and parallel platforms [11, 24, 25, 27, 37].

We also propose a sample update formula that reduces the number of BLAS-3 operations required to process a full matrix. Given a block of pivots, updating $\boldsymbol{R}$ with blocked Householder reflections requires two matrix-matrix multiplications. Without a formula to update the sample, we would need a new sample matrix after each block and one additional matrix-matrix multiply. Instead, we utilize computations that are needed to update $\boldsymbol{R}$ to also update $\boldsymbol{B}$ into a suitable sample for the next decision block.

RQRCP naturally extends to a truncated formulation, described in algorithm 5, that further reduces communication by avoiding trailing updates. For low-rank matrix approximations, this algorithm requires only one block matrix-matrix multiply per update. We accomplish this by storing reflector information in the compact WY notation [29], which allows us to construct each block of $\boldsymbol{R}$ without intermediate updates. Moreover, this algorithm serves as a key initial step in our approximation of the truncated SVD described in algorithm 6, which is a variant of Stewart’s QLP algorithm [34].

Section 2 will discuss the nature of the communication bottleneck and related approaches to address it. Section 3 will analyze sample-based pivoting as the maximization of expected utility and derive our main result, RQRCP. Section 4 will derive and explain our truncated algorithms for low-rank approximation. Section 5 will provide numerical experiments that explore the performance and decomposition quality of these approaches. Section 6 will offer concluding remarks.

Quintana-Ortí, Sun, and Bischof [30] were able to halve BLAS-2 operations with the insight that the trailing norm update does not require completing the full rank-1 update on each iteration. Instead, reflections can be gathered into blocks, as in QR. This method appears in algorithm 2.

At the end of iteration $j$, the algorithm has collected a block of reflectors $\boldsymbol{Y}^{(j)}$. Reflector $\boldsymbol{y}_{i}$, for $i\leq j$, appears in column $i$ from the diagonal down. This forms a block reflection $\boldsymbol{Q}^{(j)}=\boldsymbol{I}-\boldsymbol{Y}^{(j)}\boldsymbol{T}^{(j)}\boldsymbol{Y}^{(j)T}$, where $\boldsymbol{T}^{(j)}$ is an upper triangular $j\times j$ connection matrix that can be solved from $\boldsymbol{Y}^{(j)}$ so that $\boldsymbol{Q}^{(j)}$ is orthogonal. This algorithm must also collect each corresponding scaled inner product, $\boldsymbol{w}_{i}^{T}$, which appears as row $i$ in a matrix $\boldsymbol{W}^{(j)T}$. In the compact WY notation, this block of inner products is $\boldsymbol{W}^{(j)T}=\boldsymbol{T}^{(j)T}\boldsymbol{Y}^{(j)T}\boldsymbol{A}\boldsymbol{P}^{(j)}$, which provides enough information to update row $j$ alone and adjust trailing column norms to prepare for the next pivot selection

Note, however, that this construction complicates reflector formation. As before, the next pivot index $p_{j+1}$ is selected and swapped into column $j+1$. Let this new column be $\boldsymbol{a}_{j+1}$. From rows $j+1$ down, elements of $\boldsymbol{a}_{j+1}$ have not been updated with the current block of reflectors. Before we can form $\boldsymbol{y}_{j+1}$, prior transformations must be applied to these rows from the formula $\boldsymbol{\hat{a}}_{j+1}=\boldsymbol{a}_{j+1}-\boldsymbol{Y}^{(j)}\boldsymbol{W}^{(j)T}(:,p_{j+1})$. An additional step is also required to compute reflector inner products because they must account for reflections that have not been applied to the trailing matrix. The adjusted formula for these inner products is

The reflector and inner product blocks are then updated:

Unfortunately, the remaining BLAS-2 operations $\boldsymbol{y}_{j+1}^{T}\boldsymbol{A}^{(j)}$ and $\boldsymbol{y}_{j+1}^{T}\boldsymbol{Y}^{(j)}$ in the inner product computation still dominate slow communication complexity for large matrices. The entire trailing matrix must pass from slow to fast memory once per iteration. Consequently, even heavily optimized implementations of blocked QRCP still run substantially slower than blocked QR on both sequential and parallel architectures.

Several mechanisms have been put forward to avoid repeating full passes over the trailing matrix on each iteration. Bischof [3] proposed pivoting restricted to local blocks and Demmel, et al. [7, 8] propose a procedure called Communication Avoiding Rank-Revealing QR (CARRQR). CARRQR proceeds by partitioning the trailing matrix into $\mathcal{P}$ subsets of columns that are processed independently and possibly simultaneously. From within each column subset, $b$ candidate pivots are selected using QRCP. Adjacent subsets of candidates are then combined to form $\frac{1}{2}\mathcal{P}$ subsets of $2b$ candidates. This procedure continues, using QRCP to filter $b$ candidates per subset followed by merging results into $\frac{1}{4}\mathcal{P}$ subsets and so on, until only one subset of $b$ candidates remains. The trailing matrix is then updated as before, with blocked reflections.

We now examine several practical constraints in implementing CARRQR. First, the reflectors $\boldsymbol{Y}$, inner products $\boldsymbol{W}^{T}$, and leading rows of $\boldsymbol{R}$ must be stored separately from the original matrix for each independently processed subset of columns. Furthermore, one must employ a version of QRCP that avoids the trailing update, because the final reflectors are unknown until the last selection stage. Any intermediate changes to the original columns would have to be undone before the final transformations can be correctly processed. In contrast, QRCP can be written to convert columns into reflectors, storing the results in the same array as the input on the strictly lower triangle portion of the matrix. Likewise, $\boldsymbol{R}$ can be stored in place of the input on the upper triangle.

Depending on the initial column partition, CARRQR performs up to 2 times as many inner products as QRCP per block iteration. Note that as the reflector index $j$ increases, the total number of inner products of the form $\boldsymbol{y}_{j+1}^{T}\boldsymbol{y}_{j+1}$, $\boldsymbol{y}_{j+1}^{T}\boldsymbol{Y}^{(j)}$ and $\boldsymbol{y}_{j+1}^{T}\boldsymbol{A}^{(j)}$ remains constant. Therefore, if the $i$th subset contains $n_{i}$ columns, $bn_{i}$ inner products will be required to produce $b$ candidates. Letting $n_{1}+n_{2}+\cdots+n_{\mathcal{P}}=n$ on the first stage of refinement, summing over all of the subsets gives $bn$ inner products to produce $\mathcal{P}$ subsets of $b$ candidates. QRCP requires the same computational complexity to produce $b$ final pivots. Assuming that the number of candidates is at least halved for each subsequent stage of refinement in CARRQR, it easily follows that no more than $2bn$ inner products will be computed in total.

Despite increased computational complexity, CARRQR is intended to benefit from better memory utilization and better parallel scalability. If each column subset is thin enough to fit in fast memory, then slow communication is eliminated between iterations of $j$. The only remaining slow communication transmits pivot candidates between stages of refinement.

Unfortunately, writing and tuning CARRQR is nontrivial. We implemented this algorithm and found that it ran slightly slower than the LAPACK implementation of algorithm 2, called DGEQP3, on a shared-memory parallel machine. We believe this was mainly due to inefficient parallelization in the final stages of refinement. We assigned each column subset to a different processor, which then worked independently to produce candidates. This approach was attractive because it did not require communication between processors during each filtration stage. However, despite communication efficiency, this technique can only engage as many processors as there are column subsets. Most processors are left idle during the final stages of refinement. A second problem with this approach occurres when the matrix is too tall. In such cases, it is not possible to select column subsets that are thin enough to fit into fast memory. An efficient implementation would need alternative or additional workload-splitting tactics to use all processors at every stage of refinement.

As we will discuss in the next section, the method we propose also gathers pivots into blocks, which are then applied to the trailing matrix. Our method, however, improves performance by reducing both the communication and computational complexity needed to form a block of pivots.

Bayesian inference is not often used in numerical linear algebra. In this case, it allows us to obtain an exact expression for the uncertainty in column norms, which rigorously frames the amount of padding $p$ that is needed to satisfy a relative error bound with a desired probability of success.

We begin by expressing a single column $\boldsymbol{a}$ as its 2-norm multiplied by a unit vector $\boldsymbol{q}$, so that $\boldsymbol{a}=\left\lVert\boldsymbol{a}\right\rVert_{2}\boldsymbol{q}$. Because a GIID matrix is invariant in distribution under independent orthogonal transformations, we can construct an orthogonal matrix $\boldsymbol{Q}=[\boldsymbol{q}\;\boldsymbol{Q}_{\bot}]$ and write $\boldsymbol{\Omega}$ as

where both the leading column $\boldsymbol{\hat{\omega}}$ and remaining columns $\boldsymbol{\hat{\Omega}}$ are GIID. It easily follows that each element of the sample column, $\boldsymbol{b}=\left\lVert\boldsymbol{a}\right\rVert_{2}\boldsymbol{\hat{\omega}}$, is normally distributed with mean $0$ and latent variance $\left\lVert\boldsymbol{a}\right\rVert_{2}^{2}$.

Inferring variance from a normal distribution with a known mean is a standard problem in Bayesian statistics. The likelihood probability distribution is written $p(\boldsymbol{b}\mid\left\lVert\boldsymbol{a}\right\rVert_{2}^{2})\equiv\mathcal{N}(\boldsymbol{b}\mid 0,\left\lVert\boldsymbol{a}\right\rVert_{2}^{2}\boldsymbol{I})$ where $\boldsymbol{I}$ is the $\ell\times\ell$ identity. Jeffreys proposed the maximally uninformative prior for an unknown variance, $p(\left\lVert\boldsymbol{a}\right\rVert_{2}^{2})\equiv\left\lVert\boldsymbol{a}\right\rVert_{2}^{-2}$, which is invariant under scaling and power transformations [20]. We apply Bayes’ theorem, $p(\left\lVert\boldsymbol{a}\right\rVert_{2}^{2}\mid\boldsymbol{b})\propto p(\boldsymbol{b}\mid\left\lVert\boldsymbol{a}\right\rVert_{2}^{2})p(\left\lVert\boldsymbol{a}\right\rVert_{2}^{2})$, to obtain the posterior distribution. Normalization results in the inverse gamma distribution

Consequently, we can cast each pivoting decision as the maximizer of expected utility, where utility is taken to be the latent 2-norm squared

Since expected utility is monotonic in the sample column norms, we simply choose the maximum as in QRCP.

More importantly, the posterior distribution clearly relates sample rank $\ell$ to uncertainty. In order to capture the probability distribution of the latent relative error, we can change variables to express the latent column norm as a fraction $\phi$ of the expectation, $\left\lVert\boldsymbol{a}\right\rVert_{2}^{2}=\frac{\phi}{\ell-2}\left\lVert\boldsymbol{b}\right\rVert_{2}^{2}$. This gives analytic expressions for both the probability distribution function and the cumulative distribution function of the relative scaling

respectively. The numerator of the cumulative distribution is the upper incomplete gamma function and normalization results in the regularized gamma function.

Rank-revealing decompositions must avoid selecting columns that are already well approximated by components in the span of previous pivots, i.e columns with small trailing norms. As such, we only care about the probability that a sample column radically overestimates the true column norm. The CDF plotted in fig. 1 shows the probability that the relative scaling $\phi$ falls below a specified upper bound $\tau$ for several choices of $\ell$. Our experiments provide good results using padding $p=8$ with a block size $b=32$ so that the effective sample rank satisfies $8<\ell\leq 40$ for each pivot decision. Note that this analysis also holds for linear combinations of columns in $\boldsymbol{A}$ and $\boldsymbol{B}$, provided that those linear combinations are independent of $\boldsymbol{B}$. An in-depth reliability and probability analysis has appeared in Xiao, Gu, and Langou [37]. In particular, they show that in the case of decaying singular values in the matrix $\boldsymbol{A}$, a nearly optimal low-rank approximation can be computed with the QR factorization with a slight modification in the strategy used to choose $\boldsymbol{P}$.

We now show how the sample matrix $\boldsymbol{B}=\boldsymbol{\Omega}\boldsymbol{A}$ can be used to select a full block of $b$ pivots. If QRCP is performed on the sample matrix, then at iteration $j$ we can examine $\boldsymbol{B}$ as a partial factorization. We let $\boldsymbol{P}^{(j)}$ be the aggregate permutation so far and represent the accumulated orthogonal transformations applied to $\boldsymbol{B}$ as $\boldsymbol{U}^{(j)}$, with corresponding intermediate triangular factor $\boldsymbol{S}^{(j)}$, as shown below. We also consider a partial factorization of $\boldsymbol{A}$ using the same pivots that were applied to $\boldsymbol{B}$. The corresponding factors of $\boldsymbol{A}$ are $\boldsymbol{Q}^{(j)}$ and $\boldsymbol{R}$.

Both $\boldsymbol{S}_{11}^{(j)}$ and $\boldsymbol{R}_{11}^{(j)}$ are upper triangular. $\boldsymbol{\Omega}$ can then be expressed as elements $\boldsymbol{\hat{\Omega}}$ in the bases given by $\boldsymbol{U}^{(j)}$ and $\boldsymbol{Q}^{(j)}$:

Noting that $\boldsymbol{B}\boldsymbol{P}^{(j)}=\boldsymbol{\Omega}\boldsymbol{A}\boldsymbol{P}^{(j)}$, we have

If $\boldsymbol{S}_{11}^{(j)}$ is nonsingular, then both $\boldsymbol{\hat{\Omega}}_{11}^{(j)}$ and $\boldsymbol{R}_{11}^{(j)}$ are also nonsingular. It follows that $\boldsymbol{\hat{\Omega}}_{11}^{(j)}=\boldsymbol{S}_{11}^{(j)}\boldsymbol{R}_{11}^{(j)-1}$ is upper triangular and $\boldsymbol{\hat{\Omega}}_{21}^{(j)}=0$. In other words, we have implicitly formed a QR factorization of $\boldsymbol{\Omega}\boldsymbol{Q}^{(j)}$ using the same orthogonal matrix $\boldsymbol{U}^{(j)}$. Finally, the trailing matrix in the sample simplifies to $\boldsymbol{S}_{22}^{(j)}=\boldsymbol{\hat{\Omega}}_{22}^{(j)}\boldsymbol{R}_{22}^{(j)}$, which is a sample of the trailing matrix $\boldsymbol{R}_{22}^{(j)}$ using the compression matrix $\boldsymbol{\hat{\Omega}}_{22}^{(j)}$. If the permutation $\boldsymbol{P}^{(j)}$ were independent of the sample, then $\boldsymbol{Q}^{(j)}$ would be formed from $\boldsymbol{A}$, independent of $\boldsymbol{\Omega}$. Likewise, $\boldsymbol{U}^{(j)}$ only depends on the leading $j$ columns of the sample, which are independent of the trailing columns. Thus $\boldsymbol{\hat{\Omega}}_{22}^{(j)}$ would be GIID. As such, we may use column norms of $\boldsymbol{S}_{22}^{(j)}$ to approximate column norms of $\boldsymbol{R}_{22}^{(j)}$ when we select the $(j+1)$st pivot, so that a full block of pivots can be selected without interleaving any references to $\boldsymbol{A}$ or $\boldsymbol{R}$ memory. We note, however, that the permutation may exhibit a subtle dependence on the pivots, which we briefly discuss in section 3.3.

Once $b$ pivots have been selected from the sample matrix $\boldsymbol{B}$, the corresponding columns of $\boldsymbol{A}$ are permuted and processed all at once, as is done in BLAS-3 QR, thus reducing both the communication and computational complexity associated with selecting a block of $b$ pivots by a factor of $\ell/m$. More significantly, if the sample matrix $\boldsymbol{B}$ fits in fast memory, then slow communication between consecutive pivot decisions is eliminated within each block iteration. As in BLAS-3 QR, the remaining communication costs are due to the matrix-matrix multiplications needed to perform block reflections. As such, RQRCP satisfies the BLAS-3 performance standard.

Algorithm 3 outlines the full procedure for a sample-based block permutation. It is acceptable for very low-rank approximations, wherein the required sample is small enough to maintain communication efficiency. That is, when the desired approximation rank $k$ is small enough to be a single block, $b=k$. For larger approximations we will resort to a more comprehensive algorithm that includes a sample update formulation that subsumes this version. Since the single-sample algorithm illuminates the performance advantage gained from this approach, we examine it first.

The bias of an estimator is the difference between the expected value of the estimator and the true value of the quantity being estimated. In this case, sample column norms are used to estimate true column norms. As illustrated in the following thought experiment, a subtle form of bias occurs when we use the maximum to select a pivot.

Suppose Jack and Jill roll one six-sided die each. The expected outcome for each of them is $3.5$. We then learn that Jill rolled $5$, which was greater than Jack’s roll. This new information reduces Jack’s expectation to $2.5$, because it is no longer possible for him to have rolled $5$ or $6$. Similarly, the act of selecting the largest sample norm creates a dependency with remaining samples by truncating their plausible outcomes.

This effect becomes less pronounced, however, if the decision is more likely to be determined by the true value being estimated rather than a chance sample outcome. For example, now suppose that Jill rolls three dice and Jack rolls one. If we only know each person’s sum, we can still infer the number of dice each person rolled and select the expected maximum as before. In $90.7\%$ of trials, Jill’s sum will be $7$ or greater, driven by the fact that three dice were rolled. In most cases, this leaves Jack’s potential outcomes unconstrained and his expectation unbiased after we observe the maximum. Analogously, when a selected column exhibits a norm that is an order of magnitude greater than others, the bias effect is negligible and we may proceed as though the pivot decision was independent of the sample.

The original version of this paper included analysis of this distribution truncation effect [9]. For our purposes, we proceed as though the progression of sample updates is independent of the pivot decisions, which is also the approach Xiao, et al. [37] take. This assumption is potentially problematic when multiple trailing columns have similar norms, but any such column provides a suitable pivot in that scenario.

The sample matrix in algorithm 3 is formulated to have rank $\ell=k+p$, where $k$ was both the desired approximation rank and the permutation block size, $b$. As $k$ increases, however, it becomes inefficient to simply increase the sample rank $\ell$. In the extreme case, a full decomposition would require a sample just as big as the original matrix. If we require a decomposition with a larger rank than that which can be efficiently sampled and blocked, that is if the sample rank cannot exceed $\ell=b+p$, but we require $k>b$, then we need to update the sample matrix after each block. Martinsson [28] developed an approach in which one simply processes each subsequent block by drawing a new random matrix $\boldsymbol{\Omega}$ and applying it to each new trailing matrix. We propose a sample update formulation that does not require multiplying the trailing matrix by a new compression matrix and reduces BLAS-3 communication in the overall factorization by at least one third.

The update formula we derive is an extension of the implicit update mechanism described in the previous section. Both algorithm 4 and algorithm 5, the truncated variation, will proceed in blocks of pivots. Bracket superscripts denote the results of a computation that occurred on the indicated block-iteration. At entry to the first block-iteration, the sample is $\boldsymbol{B}^{[0]}=\boldsymbol{\Omega}^{[0]}\boldsymbol{A}^{[0]}$, where $\boldsymbol{A}^{[0]}$ is the original matrix. At the end of block-iteration $J$, the sample will be in the transformed state

just as in eq. 1. $\boldsymbol{S}_{11}^{[J]}$ is the leading upper triangle from the partial factorization of the sample and $\boldsymbol{S}_{22}^{[J]}$ gives the trailing sample columns. Likewise $\boldsymbol{R}_{11}^{[J]}$ and $\boldsymbol{A}^{[J]}$, respectively, give the leading upper triangle and trailing columns that would be obtained by factorizing the original matrix with the same pivots. By absorbing the transformations $\boldsymbol{U}^{[J]T}$ and $\boldsymbol{Q}^{[J]}$ into $\boldsymbol{\Omega}^{[J]}$, we obtained an effective compression matrix $\boldsymbol{\hat{\Omega}}_{22}^{[J]}$, which had already been implicitly applied to the trailing columns: $\boldsymbol{S}_{22}^{[J]}=\boldsymbol{\hat{\Omega}}_{22}^{[J]}\boldsymbol{A}^{[J]}$. The difficulty is $\boldsymbol{S}_{22}^{[J]}$ only has rank $p$. In order to construct a rank $\ell=b+p$ sample of the trailing matrix $\boldsymbol{A}^{[J]}$, we need to include $\boldsymbol{\hat{\Omega}}_{12}^{[J]}$ in the updated compression matrix

In other words, the new compression matrix $\boldsymbol{\Omega}^{[J]}$ is simply $\boldsymbol{U}^{[J]T}\boldsymbol{\Omega}^{[J-1]}\boldsymbol{Q}^{[J]}$, with the leading $b$ columns removed. This new compression matrix does not need to be explicitly formed or applied to the trailing columns $\boldsymbol{A}^{[J]}$. Instead, we form the result implicitly by removing $\boldsymbol{\hat{\Omega}}_{11}^{[J]}\boldsymbol{R}_{12}^{[J]}$ from $\boldsymbol{S}_{12}^{[J]}$. Both $\boldsymbol{R}_{11}^{[J]}$ and $\boldsymbol{R}_{12}^{[J]}$ will be computed in blocked matrix multiply operations using the previous $b$ pivots of $\boldsymbol{A}$. Since $\boldsymbol{\hat{\Omega}}_{11}^{[J]}$ can then be recovered from $\boldsymbol{S}_{11}^{[J]}$, we can avoid any direct computations on $\boldsymbol{\Omega}$. We only need to update the first $b$ rows of $\boldsymbol{B}$, which gives us the sample update formula

Full RQRCP, described in algorithm 4, can be structured as a modification to blocked BLAS-3 QR. The algorithm must simply interleave processing blocks of reflectors with permutations obtained from each sample matrix, then update the sample as above. To obtain a modified version that employs repeated sampling simply replace the sample update with a new sample of the trailing matrix, $\boldsymbol{B}^{[J]}=\boldsymbol{\Omega}^{[J]}\boldsymbol{A}^{[J]}$.

When QRCP is applied to the sample matrix $\boldsymbol{B}$, only a partial decomposition is necessary. The second argument $b$ in the subroutine call $\texttt{QRCP}(\boldsymbol{B}^{[J]},b)$ indicates that only $b$ column permutations are required. It is relatively simple to modify the QR algorithm to avoid any unnecessary computation. After sample pivots have been applied to the array containing both $\boldsymbol{A}$ and $\boldsymbol{R}$, we perform QR factorization on the new leading $b$ columns of the trailing matrix. Although this is stated as returning $\boldsymbol{Q}^{[J]}$ for convenience, it can be implemented efficiently with the blocked Householder reflections described in section 2.1. We can then apply reflectors to the trailing matrix and form the sample update $\boldsymbol{B}^{[J]}$ to prepare for the next iteration.

Our technique is analogous to the method Quintana-Ortí, Sun, and Bischof used to halve BLAS-2 operations in QRCP. In their version of QRCP, all reflector inner products are computed, but rows and columns are only updated as needed. In order to compute correct reflector inner products, without having updated the trailing matrix, we need to formulate blocked reflector compositions

The corresponding reflector inner products $\boldsymbol{W}^{T}=\boldsymbol{T}^{T}\boldsymbol{Y}^{T}\boldsymbol{A}$ become

If we store these reflector inner products, then we can construct any submatrix of the accumulated transformation $\boldsymbol{\hat{A}}^{[J]}=\boldsymbol{A}-\boldsymbol{Y}^{[J]}\boldsymbol{W}^{[J]T}$ as needed. Columns that are selected by sample pivots are constructed just before becoming the next reflectors and corresponding rows of $\boldsymbol{R}$ are constructed just before being used to update the sample, as outlined in algorithm 5.

TRQRCP naturally extends to an approximation of the truncated SVD by following the QLP method proposed by Stewart. The QLP decomposition proceeds by first applying QRCP to obtain $\boldsymbol{A}\boldsymbol{P}_{0}=\boldsymbol{Q}_{0}\boldsymbol{R}$. Then the right triangular matrix $\boldsymbol{R}$ is factored again using an LQ factorization $\boldsymbol{P}_{1}\boldsymbol{R}=\boldsymbol{L}\boldsymbol{Q}_{1}$, where row-pivoting is an optional safeguard (otherwise $\boldsymbol{P}_{1}=\boldsymbol{I}$), giving the factorization $\boldsymbol{A}=(\boldsymbol{Q}_{0}\boldsymbol{P}_{1}^{T})\boldsymbol{L}(\boldsymbol{Q}_{1}\boldsymbol{P}_{0}^{T})$. The diagonal elements of $\boldsymbol{L}$ approximate the singular values of $\boldsymbol{A}$. Huckaby and Chan [18] provide convergence analysis.

The approximate truncated SVD we propose (TUXV) simply adapts low-rank versions of the steps in QLP. The rank-$k$ approximation that results is exactly the same as the truncated approximation that would be obtained if QLP had been processed to completion using RQRCP, without secondary row-pivoting, and then truncated to a rank-$k$ approximation.

We begin by using TRQRCP to produce $k$ left reflectors, which defines the initial left orthogonal matrix $\boldsymbol{U}^{(0)}$. Superscript $(0)$ refers to the initial state of an array upon entry to the first iteration of the main loop. We can compare our results to what would have been obtained from full RQRCP-based QLP:

The asterisk denotes additional pivoting produced from the full factorization. The first $k$ pivots in $\boldsymbol{P}^{(0)}$ and corresponding reflectors in $\boldsymbol{U}^{(0)}$ are the same, as are corresponding rows in $\boldsymbol{R}^{(0)}$ modulo additional column permutations. We reverse these permutations to construct the $k\times n$ matrix $\boldsymbol{Z}^{(0)}=\boldsymbol{R}^{(0)}\boldsymbol{P}^{(0)T}$ so that

Taking the LQ factorization from $\boldsymbol{Z}^{(0)}$, so that $\boldsymbol{L}^{(1)}\boldsymbol{V}^{(1)T}=\boldsymbol{Z}^{(0)}$, instead of from $\boldsymbol{R}^{(0)}$ simply absorbs the permutation $\boldsymbol{P}^{(0)T}$ into the definition of $\boldsymbol{V}^{(1)T}$ so that

For consistency with the following algorithm, we label the $k\times k$ connecting matrix $\boldsymbol{L}^{(1)}_{11}$ as $\boldsymbol{X}^{(0)}$. We will discuss the connecting matrix further after explaining the rest of the algorithm. Provided that no secondary row-pivoting is considered, the leading $k$ reflectors in $\boldsymbol{V}^{(1)}$ and $\boldsymbol{V}^{(1*)}$ are identical because they are only computed from the leading $k$ rows of $\boldsymbol{Z}^{(0)}$. At this point, the rank-$k$ approximation of RQRCP-based QLP would require $\boldsymbol{L}^{(1*)}_{21}$, which is unknown. Fortunately, the leading $k$ columns of $\boldsymbol{U}^{(0*)}\boldsymbol{L}^{(1*)}$ can be reconstructed with one matrix multiply. We label this $m\times k$ matrix $\boldsymbol{Z}^{(1)}$, which is