Iterative Refinement and Oversampling for Low Rank Approximation

TL;DR

This paper extends iterative refinement to low rank matrix approximation, linking it to oversampling techniques, and demonstrates that using sparse sketches enables very fast, low-cost approximations with near-optimal accuracy for certain matrices.

Contribution

It introduces a novel iterative refinement approach for low rank approximation, connecting it to oversampling, and shows how sparse sketch matrices enable efficient, very low rank approximations with improved accuracy.

Findings

Sparse sketch matrices enable sublinear cost low rank approximation.

Oversampling improves accuracy to near-optimal levels.

The method is effective for matrices with fast decaying singular values.

Abstract

Iterative refinement is particularly popular for numerical solution of linear systems of equations. We extend it to Low Rank Approximation of a matrix (LRA) and observe close link of the resulting algorithm to oversampling techniques, commonly used in randomized LRA algorithms. We elaborate upon this link and revisit oversampling and some efficient randomized LRA algorithms. Applied with sparse sketch matrices they run significantly faster and in particular yield Very Low Rank Approximation (VLRA) at sublinear cost, using much fewer scalars and flops than the input matrix has entries. This is achieved at the price of deterioration of output accuracy, but according to our formal and empirical study subsequent oversampling improves accuracy to near-optimal level under the spectral norm for a large sub-class of matrices with fast decaying spectra of singular values.

Tables

Table 1. Table 6.1: Test results for Alg. 4.1 on Gravity, discretized single layer potential operator (SLP), synthetic matrices with fast decaying spectra (Fast Decay), and synthetic matrices with slow decaying spectra (Slow Decay).

Input Matrix	$\rho$
	$2r$	$3r$	$4r$	$5r$
Gravity ($r=45$)	1.000 $\pm$ 1.036E-05	1.000 $\pm$ 7.815E-06	1.000 $\pm$ 7.888E-06	1.000 $\pm$ 6.683E-06
SLP ($r=11$)	1.970 $\pm$ 4.083	1.000 $\pm$ 0.002	1.000 $\pm$ 1.018E-07	1.000 $\pm$ 1.362E-10
Fast Decay ($r=20$)	1.000 $\pm$ 3.086E-12	1.000 $\pm$ 3.590E-16	1.000 $\pm$ 2.470E-16	1.000 $\pm$ 2.617E-16
Slow Decay ($r=20$)	1.000 $\pm$ 3.604E-05	1.000 $\pm$ 1.244E-06	1.000 $\pm$ 1.543E-07	1.000 $\pm$ 3.754E-08

Table 2. Table 6.2: Test results for Alg. 5.1 on synthetic matrices with fast decaying spectra (Fast Decay), synthetic matrices with slow decaying spectra (Slow Decay), Shaw, Gravity, and discretized single layer potential operator (SLP).

Input Matrix	Multiplier Type	1st Itr.	2nd Itr.		3rd Itr.
			before	after	before	after
Fast Decay	Abr. Had.	3.1550E+00	2.9872E-11	1.0000E+00	2.4894E-11	1.0000E+00
	Gaussian	3.1202E+00	1.5322E-11	1.0000E+00	2.4075E-11	1.0000E+00
Slow Decay	Abr. Had.	5.0468E+00	3.3300E-02	1.0003E+00	3.1365E-02	1.0001E+00
	Gaussian	5.0755E+00	3.1995E-02	1.0002E+00	3.1017E-02	1.0001E+00
Shaw	Abr. Had.	2.8820E+01	5.3133E-01	1.0983E+00	5.4956E-01	1.1225E+00
	Gaussian	1.8235E+01	6.2568E-01	1.1517E+00	5.5411E-01	1.1189E+00
Gravity	Abr. Had.	1.5762E+01	7.6046E-03	1.0000E+00	8.4108E-03	1.0000E+00
	Gaussian	1.2917E+01	9.7073E-03	1.0000E+00	8.7926E-03	1.0000E+00
SLP	Abr. Had.	1.0931E+02	3.3201E-02	1.0014E+00	1.0484E-02	1.0000E+00
	Gaussian	5.2205E+00	4.1354E-03	1.0000E+00	4.6096E-03	1.0000E+00

Equations

$$M=\tilde{M}+E\mathrm{where}||E||\le ϵ||M||,\tilde{M}=AB,A\in {\mathbb{C}}^{m\times r},\mathrm{and}B\in {\mathbb{C}}^{r\times n}.$$

$$M=\tilde{M}+E,||E||\le ϵ,\tilde{M}=XTY,X\in {\mathbb{C}}^{m\times k},T\in {\mathbb{C}}^{k\times l},Y\in {\mathbb{C}}^{l\times n},$$

$$r=\mathrm{rank}(\tilde{M})\le k\le m,r\le l\le n.$$

$$|{\sigma }_j(M+E)-{\sigma }_j(M)|\le ||E||\mathrm{for}j=1,\dots ,\min \{m,n\}.$$

$${\nu }_{\rho }:=||M-M(\rho )||\ll ||M-M_r||$$

$$||E||\le ||M-M(\rho )||+||M(\rho )-(M(\rho ){)}_r||=||M(\rho )-(M(\rho ){)}_r||+{\nu }_{\rho }.$$

$$||M(\rho )-(M(\rho ){)}_r||\le ||M-M_r||+{\nu }_{\rho }.$$

$$\text{err}=\frac{||M-(M(\rho ){)}_r||}{||M-M_r||}$$

$$\text{err}=\frac{||M-X_i||}{||M-M_r||}$$

$$|||M-(AB{)}_{\rho }|||\le |||M-AB|||+|||AB-(AB{)}_{\rho }|||=|||M-AB|||+{\tau }_{\rho }(AB)$$

$$|||M-(AB{)}_{\rho }|||\le 2|||M-AB|||+{\tau }_{\rho }(M).$$

$$H_{d,0}:=I_{n/2^d},H_{d,i+1}:=\left(\begin{array}{cc}{H}_{d,i}& H_{d,i}\\ H_{d,i}& -H_{d,i}\end{array}\right)\mathrm{for}i=0,1,\dots ,d-1,$$

Full text

Key Words:

near-optimal algorithms, superfast (sublinear cost) algorithms, low rank approximation, iterative refinement

2020 Math. Subject Classification:

65F55, 65N75, 65Y20, 68Q25, 68W20

1 Introduction

1.1 Superfast algorithms by means of LRA:

benefits and an obstacle

In modern computations with matrices representing Big Data it is highly desirable to apply algorithms that run superfast, aka in sublinear time, that is, use much fewer memory cells and flops than an input matrix has entries.111Flop is a floating point arithmetic operation. Superfast algorithms (aka running in sublinear arithmetic time) are studied extensively in for Toeplitz, Hankel, Vandermonde, Cauchy, and other structured matrices having small displacement rank [22], [24], and we extend this concept to other matrices. Quite frequently such matrices admit their low rank approximations (LRAs), with which we can operate superfast. By virtue of Eckart-Young-Mirsky theorem, optimal rank-$r$ approximation of a matrix under both spectral and Frobenius matrix norms can be readily obtained from its SVD, but computation of SVD is costly and not superfast (see Sec. 2 for background). Can we still compute a close LRA superfast if a matrix admits one? Unfortunately, any superfast algorithm at best halves the large error of the trivial approximation by the matrix filled with 0s, even for an input family of $mn+1$ matrices having size $m\times n$ and rank at most 1 (see Appendix A or, say, [27, 28, 30, 29]). If randomization is applied, such an algorithm still at best halves the error of the trivial approximation with a probability not close to 0.

1.2 Recent and new progress

This is not the end of the story, however, because nontrivial superfast algorithms of [17, 2, 11, 13] compute near-optimal LRAs of matrices of some large and important classes. For surprisingly simple extension of this impressive progress to yet another class of matrices, we propose superfast algorithms that (i) first compute a crude initial LRA and (ii) then escalate it to near-optimal level.

1.3 Superfast computation of a crude LRA

The algorithm of [37] outputs a near-optimal LRA with a high probability (whp) under random sampling performed by means of pre- and/or post-multiplication of an input matrix $M$ by Gaussian (normal) random test matrices ${\bf\Omega}$ and ${\bf\Psi}$, respectively.222We call a random matrix Gaussian if it is filled with independent standard Gaussian (normal) random variables. Moreover, the algorithm is superfast except for this bottleneck multiplication stage.

We give up formal support for near-optimality to run the algorithm superfast; namely, we propose to replace its test matrices ${\bf\Omega}$ and ${\bf\Psi}$ with sparse multipliers. Then the output approximation errors tend to grow (see [12]) but still tend to stay reasonably low according to our study where we use abridged SRHT matrices for the sparse multipliers (see Secs. 3 and 6 and Appendix C).

1.4 Escalation to superfast near-optimal LRAs

Our superfast algorithms escalate crude LRAs to near-optimal LRAs in two ways – for two large classes of matrices, specified implicitly in Secs. 4 – 6, in terms of the spectra of their singular values and/or by the results of our numerical tests.

Further advance of our pioneering algorithms and of their formal and empirical support are natural challenges as well as further extension of our approach to tensor decomposition.

1.5 Related works

In our choice of abridged SRHT matrices for sparse multipliers we follow [27, 28, 29, 30], which proved some crude but reasonable probabilistic upper bounds on output errors of some superfast LRA algorithms under randomization in the input matrix space defined by a fixed probability structure in that space.333These papers continued the study in [32, 31, 34], which introduced a probability structure in the space of input matrices to motivate application of Gaussian elimination without pivoting.

In Sec. 5 we extend to LRA the classical iterative refinement techniques, which are popular for the solution of a linear system of equations (see [36, Secs. 3.3.4 and 4.2.5], [9, Ch. 12 and Sec. 20.5], [3], [7, Secs. 3.5.3 and 5.3.8], [1, Secs. 1.4.6, 2.1.4, and 2.3.7]) as well as for various other linear and nonlinear matrix computations (see [36, page 223 and the references on page 225]). Our extension is not straightforward – we incorporate the re-compression techniques of [19, 20, 25, 33, 26, 23] and [22, Secs. 6.3-6.5].

We know of no other attempts of devising iterative refinement algorithms for LRA except for the paper [14], whose algorithm uses no crude auxiliary LRA but instead bounds the Principal Angle Distances between singular spaces of $M$ and its LRAs and incorporates randomization based on the leverage scores of [4].

1.6 Organization

of our paper

In the next section we recall some background material. We specify sample superfast computation of a crude LRA in Sec. 3 and present our superfast LRA escalators in Secs. 4 and 5. We cover our test results in Sec. 6 and devote Sec. 7 to conclusions. In Appendix A we describe some small families of matrices for which any superfast LRA algorithm fails. In Appendix B we recall optimal and nearly optimal algorithms for LRA of a matrix based on computing its top SVD and rank-revealing QRP factorization, respectively. In Appendix C we generate abridged SRHT matrices, which we use for devising superfast LRA algorithms by means of random sampling.

2 Background

$||M||:=||M||_{2}$ denotes the spectral norm of a matrix $M$, to which we restrict our study.

Definition 2.1**.**

An $m\times n$ matrix $M$ has $\epsilon$-rank* at most $r$ for a fixed tolerance $\epsilon\geq 0$ if and only if $M$ admits its approximation within an error norm $\epsilon$ by a matrix $\tilde{M}$ of rank at most $r$ or equivalently if and only if there exist three matrices $A$, $B$ and $E$ such that*

[TABLE]

A matrix admits its LRA if and only if it has $\epsilon$-rank $r$ where $\epsilon$ and $r$ are small in context.

For a small positive $\epsilon$ one can regard $\epsilon$-rank as numerical rank (cf. [7, pages 275-276]), being aware that $\epsilon$-rank is numerically unstable where multiple singular values are close to $\epsilon$.

A 2-factor LRA $AB$ of $M$ of (2.1) can be generalized to a 3-factor LRA,

[TABLE]

[TABLE]

$\sigma_{j}(M)$ denotes the $j$th largest singular value of a matrix $M$, $\sigma_{1}(M)=||M||$.

For a matrix $M$ define its $r$-top SVD $M_{r}=U_{r}\Sigma_{r}V_{r}^{*}$ where

• •

$\Sigma_{r}=\operatorname{diag}(\sigma_{j}(M))_{j=1}^{r}$ is the diagonal matrix of the $r$ top (largest) singular values of $M$,

• •

$U_{r}$ and $V_{r}$ are two matrices of the associated top singular spaces of $M$, having orthonormal columns,

• •

we call the matrix $M_{r}$ the rank-$r$ truncation of the SVD of $M$ and can define that matrix by setting to 0 all singular values of $M$ except for the $r$ largest ones.

$M_{r}$ is optimal rank-$r$ approximation of $M$ under both spectral and Frobenius norms by virtue of the Eckart-Young-Mirsky theorem (cf. [35, 5, 16], [7, Thm. 2.4.8]); we only state it under the spectral norm:

Theorem 2.1**.**

$\sigma_{r+1}(M)=||M_{r}-M||=\min_{N:~{}\operatorname{rank}(N)\leq r}||M-N||$.

Theorem 2.2**.**

[7, Cor. 8.6.2].* For a pair of $m\times n$ matrices $M$ and $M+E$ it holds that*

[TABLE]

$M^{+}$ denotes the Moore–Penrose pseudo inverse of $M$.

Hereafter we write $a\ll b$ to show that $b$ exceeds $a$ greatly, in context, and let $m\geq n$; otherwise, if $m<n$, we can extend our study from $M$ to $M^{T}$, the transpose of $M$, because $(M_{r})^{T}=(M^{T})_{r}$.

For an $m\times n$ matrix $M$ and four integers $m,n,r$, and $\rho$ such that $1\leq r<\rho\leq n\leq m$, consider the following Assumptions:

(a) a fixed algorithm $\mathbb{A}$ computes a (crude but) reasonably close rank-$\rho$ approximation $M(\rho)$ of $M$ superfast (see Secs. 1.3 and 3) and

(b) $\rho\ll n$ and $\sigma_{r+1}(M)\gg\sigma_{\rho+1}(M),$ that is, the singular values $\sigma_{j+1}(M)$ decay strongly as $j$ increases from $r$ to $\rho$.

3 Superfast computation of a crude LRA

The following algorithm supports Assumption (a) by using abridged SRHT matrices of Appendix C for sparse multipliers.

Algorithm 3.1**.**

Input:

Three integers $m$, $n$, and $r$ satisfying $0<2r\leq n\leq m$; an $m\times n$ matrix $M$, given explicitly or implicitly.

Initialization:

1. 1.

Generate a pair of independent abridged SRHT matrices $F$ and $H$ of length 3 and sizes $2r\times m$ and $n\times r$, respectively (see Appendix C); 2. 2.

Compute the sketches $FM$ and $MH$.

Computations:

1. 1.

Compute matrix $Q$ as the Q factor of the thin QR factorization444Recall the uniqueness of the thin QR factorization of an $m\times k$ matrix where an $k\times k$ upper triangular factor $R$ has positive diagonal entries [7, Thm. 5.3.2].* of $EH$.* 2. 2.

Compute the matrices $U$ and $T$ of the thin QR factorization $UT$ of $FQ$. 3. 3.

Compute and output the matrix $\Delta=QT^{+}U^{T}(FM)$, an LRA of $M$.

4 LRA via a single escalation step

Algorithm 4.1: LRA via a single escalation step.

INPUT: An $m\times n$ matrix $M$ with $m\geq n$, two integers $r$ and $\rho$, and Algorithm $\mathbb{A}$, such that above Assumptions (a) and (b) hold.

COMPUTATIONS:

1. Apply Algorithm $\mathbb{A}$.

2. Compute and output a rank-$r$ truncation $(M(\rho))_{r}$ of its output matrix $M(\rho)$.

Assumption (b) implies that

[TABLE]

or, in words, even a crude rank-$\rho$ approximation $M(\rho)$ lies much closer to $M$ than its optimal rank-$r$ approximation $M_{r}$.

Algorithm $\mathbb{A}$ computes $M(\rho)$ superfast under Assumption (a). If rank $\rho$ is low, then Alg. 3.1 is superfast at stage 2 as well (see Appendix B).

Theorem 4.1**.**

For the error matrix $E:=M-(M(\rho))_{r}$ of Alg. 4.1 and $\nu_{\rho}$ of (4.1) it holds that $||E||\leq||M-M_{r}||+2\nu_{\rho}$ (cf. Thm. 2.1).

Proof.

Notice that $E=M-(M(\rho))_{r}=M-M(\rho)+(M(\rho)-(M(\rho))_{r},$ and so

[TABLE]

Eqn. (4.1) and Thm. 2.2 together imply that $\sigma_{j}(M(\rho))\leq\sigma_{j}(M)+\nu_{\rho}~{}{\rm for~{}all}~{}j.$ Hence

[TABLE]

Combine estimates (4.2) and (4.3). ∎

Thm. 4.1 implies that $||M-(M_{\rho})_{r}||\leq||M-M_{r}||+2\nu_{\rho}$ for $\nu_{\rho}$ of (4.1). Hence under the spectral norm $||\cdot||$ the rank-$r$ approximation $(M_{\rho})_{r}$ of $M$, output by Alg. 4.1, is optimal up to the term $2\nu_{\rho}:=2||M-M(\rho)||$, that is, near-optimal under (4.1).

5 Iterative refinement of an LRA with recursive re-compression

5.1 Introductory comments

In our alternative LRA algorithms we rely solely on Assumption (a) and assume no a priori information about the spectrum of the singular values of the matrix $M$, although in view of the Eckart-Young-Mirsky theorem the efficiency of the algorithm should depend on that spectrum.

We propose to apply Algorithm $\mathbb{A}$ recursively (for a certain policy of recursive selection of integers $\rho$) to the error matrices output at the previous recursive steps rather than to $M$. In other words, we escalate a crude LRA towards optimal LRA recursively – by means of extending to LRA the popular technique of iterative refinement.

To get progress in refinement we apply the customary mixed precision recipe, that is, increase (typically double) the computational precision at all subtraction stages. In our superfast iterative refinement of LRA we only compute sketches $FE$ and/or $EH$ for abridged SRHT multipliers $F$ and/or $H$ (or for some other tall, skinny, and sparse matrices $F^{T}$ and/or $H$) rather than all entries of error matrices $E$.

To extend the known refinement techniques to LRA, we also apply mixed compression recipe – by incorporating the re-compression techniques of [19, 20, 25, 33, 26, 23] and [22, Secs. 6.3-6.5]). Namely, in refinement iteration we approximate $M$ by the sums of all rank-$r$ approximations of the error matrices computed so far, beginning with the initial approximation of $M$ at the first iteration. In the $k$th iteration the sum can have rank as large as $kr$; its compression to optimal rank-$r$ approximation can be costly for large $k$, but we periodically compress the partial sums to their rank-$r$ approximations keeping bounded the number of rank-$r$ terms in every sum, e.g., just two rank-$r$ terms.

Using sketches and compression complicates formal support for iterative refinement LRA algorithms for LRA, but empirically they handle a wider class of inputs than Alg. 4.1.

5.2 Two refinement algorithms

Next we specify two variants of iterative refinement of LRA.

Algorithm 5.1: ITERATIVE REFINEMENT 1.

INPUT: an $m\times n$ matrix $M$, an integer $r$ such that $1\leq r\ll n\leq m$, and a black box algorithm $\mathbb{A}=\mathbb{A}_{\rho}$ that for a positive integer $\rho\leq n$ and an $m\times n$ matrix $E$ outputs a crude but reasonably close rank-$\rho$ approximation $X$ of $E$.

INITIALIZATION: Write $X_{0}:=0$ and $E_{0}:=M$; fix integers $k>0$ and $\rho_{i}\ll n$ for $i=1,\dots,k$ such that $\rho_{1}\geq r$ and $\rho_{i}>r$ for $i>1$.

COMPUTATIONS: Do for $i=1,2,\dots,k$:

1. $Y_{i}=\mathbb{A}_{\rho_{i}}(E_{i-1})$,

2. $X_{i}=(X_{i-1}+Y_{i})_{r}$,

3. $E_{i}=M-X_{i}$.

OUTPUT: $X:=(X_{k})_{r}$.

Next we update the error matrix as follows: $E_{i-1}-Y_{i}\mapsto E_{i}$. This should counter cancellation of significant digits in the subtraction $M-X_{i}$.

Algorithm 5.2: ITERATIVE REFINEMENT 2.

INPUT and INITIALIZATION: as in Alg. 5.1.

COMPUTATIONS: Do for $i=1,2,\dots,k$:

1. $Y_{i}=\mathbb{A}_{\rho_{i}}(E_{i-1})$,

2. $E_{i}=E_{i-1}-Y_{i}$.

OUTPUT: $X:=(Y_{1}+Y_{2}+\cdots+Y_{k})_{r}$.

Remark 5.1**.**

We avoid $k$-fold memory increase in the output summation of Alg. 5.2 by means of recursive truncation. We can just write $X_{1}:=Y_{1}$, $X_{i}:=(X_{i-1}+Y_{i})_{r}$ for $i=2,3,\dots,k$, $X:=X_{k}$ or, more generally, can recursively combine summation and periodic truncation with a fixed period $h$, $2\leq h\leq(k-1)/2$, as follows: $X_{0}:=Y_{1}$, $X_{i}:=X_{i-1}+\sum_{j=1}^{h-1}Y_{(i-1)(h-1)+j+1}$ for $i=1,\dots,q$, for $q:=\lceil(k-1)/(h-1)\rceil$ (in particular, $X_{1}:=(X_{0}+Y_{2}+\cdots+Y_{h})_{r}$), $X:=X_{q}$.

Remark 5.2**.**

Algs. 5.1 and 5.2 stop in $k$ iterations, where one can choose the integer parameter $k$ based on the previous record of convergence of the output matrices $X$ to $M_{r}$; our tests of Sec. 6 have consistently succeeded in two iterations. For an alternative heuristic stopping criterion, we can stop iterations where a fixed norm (say, 1-norm or $\infty$-norm) of the matrices $E_{i}H$ and $FE_{i}$ decreases below a fixed tolerance value for sufficiently many random choices of sparse multipliers $F$ and $H$. Optimization of stopping criteria for the refinement algorithms as well as their formal support and cost estimates are research challenges. For some test results of Alg. 5.1 see Sec. 6.

6 Numerical experiments

In this section we present the test results for our LRA algorithms on inputs made up of synthetic and real-world data with various spectra. We implemented our algorithm in Python, running it on a 64bit MacOS machine with 16GB Memory. We called scipy.linalg for numerical linear algebra routines such as QR factorization with pivoting, Moore-Penrose matrix inversion, and SVD.

6.1 Input Matrices

Real-world Input Matrices: We used three input matrices coming from different applications.

The matrix Shaw has target rank $\rho=20$ and comes from a one-dimensional image restoration model problem.

The matrix Gravity has target rank $\rho=45$ and comes from a one-dimensional gravity surveying model problem.

Both Shaw and Gravity are $1000\times 1000$ dense real matrices having low numerical rank; we padded them with 0s to increase their size to $1024\times 1024$. They represent discretization of Integral Equations, provided in the built-in problems of the Regularization Tools555See http://www.math.sjsu.edu/singular/matrices and http://www.imm.dtu.dk/$\sim$pcha/Regutools For more details see Ch. 4 of http://www.imm.dtu.dk/$\sim$pcha/Regutools/RTv4manual.pdf .

The third input matrix is generated from the discretization of a single layer potential (SLP) operator (see [10, Sec. 7.1] for more details). It has size $1024\times 1024$ and the target rank $\rho=11$.

We display the distribution of the top 50 singular values of these matrices in Fig. 1.

Synthetic Input Matrices: We used random synthetic $1024\times 1024$ input matrices of two kinds – with fast and slowly decaying spectra. In both cases we generated these matrices as the products $U\Sigma V^{T}$, where $U$ and $V$ were the matrices of the left and right singular vectors of a Gaussian random matrix, respectively.

By letting $\Sigma=\textrm{diag}(v)$, where $v_{i}=1$ for $i=1,2,3,\dots 20$, $v_{i}=\frac{1}{2^{i-20}}$ for $i=21,\dots,100$, and $v_{i}=0$ for $i>100$, we generated input matrices with fast decaying spectra and target rank $\rho=20$.

By letting $\Sigma=\textrm{diag}(u)$, where $u_{i}=1$ for $i=1,2,3,\dots 20$, and $u_{i}=\frac{1}{(1+i-20)^{2}}$ for $i>20$, we generated input matrices with slowly decaying spectra and target rank $\rho=20$.

6.2 Test results for Alg. 4.1

Table 6.1 displays the relative error

[TABLE]

for an input matrix $M$, a target rank $r$, and the $r$-truncation $M_{r}$ of $M$ and $(M(\rho))_{r}$ of $M(\rho))$ (cf. Alg. 4.1). We tested $m\times n$ matrices $M$ for $m=n=1024$ and let $\rho=2r$, $3r$, $4r$, and $5r$, to see the impact of increasing $\rho$.

We repeated the tests 100 times per test set, that is, per triple of $r$, $\rho$, and an input matrix type. Table 6.1 shows the mean and standard deviation of the 100 relative errors.

For our input classes, the relative error bound usually stayed near the optimal value 1.000 already where we choose $\rho=2r$, and the errors consistently decreased towards 1.000 as $\rho$ increased (see Table 6.1). Alg. 4.1, however, failed on the matrix class Shaw for $r=20$, which could be predicted because the spectrum of singular values $\sigma_{j}(M)$ of these matrices $M$ is flat for $j>20$ (decreasing roughly from $10^{-13}$ to between $10^{-15}$ and $10^{-16}$). The test results were much more favorable for $r=5$ and $r=15$.

6.3 Test results for iterative refinement of LRA

Our Table 6.2 displays the relative error

[TABLE]

for an input matrix $M$, the approximation $X_{i}$ output in the $i$-th iteration of Alg. 5.1, a target rank $r$, and the $r$-truncation $M_{r}$ of SVD of $M$.

We tested Alg. 5.1 for iterative refinement of an LRA where we incorporated Alg. 3.1 as Algorithm $\mathbb{A}$, applied recursively for superfast computation of crude LRA and using abridged SRHT multipliers $F\in\mathbb{R}^{2r\times 1024}$ and $H\in\mathbb{R}^{1024\times r}$.

For comparison, we have also tested Alg. 5.1 incorporating as Algorithm $\mathbb{A}$ the slower algorithm of [37] with Gaussian random multipliers of the same sizes.

To control rank growth of the auxiliary LRAs computed by Algorithm $\mathbb{A}$, we applied them for $\rho_{1}=r$ and $\rho_{i}=2r$ for all $i>1$.

For each input matrix and each type of multipliers, we run Alg. 5.1 one hundred (100) times with independent random choices of multipliers $F$ and $H$. We terminated each run in three refinement iterations, even though the approximations were always quite accurate already in two iterations. We recorded the mean relative error ratio of the approximation before and after the truncation for every iteration step in Table 6.2. With our choice of subalgorithm and multipliers, the rank of approximation before truncation could potentially be twice the target rank $r$ in the second and third iterations, and for this reason the relative error before truncation were much less than 1 in our tests.

Running LRA tests with both Gaussian random and abridged SRHT multipliers we observed similar output errors with only minor deterioration as the price for superfast acceleration.

Even in the tests with the matrices Shaw, the output relative error was close to 1, although not as close as in case of the other inputs, possibly because of the adverse impact of the clusters of trailing singular values in the spectra of matrices Shaw.

7 Conclusions

We proposed, tested, and briefly analyzed some simple but novel superfast LRA algorithms that are near-optimal for a large class of matrices.

Our study should motivate new effort for further advance in this direction. In particular, one can extend Alg. 3.1 by recalling some known algorithms that compute accurate LRA (e.g., those of [8, 18, 30, 14]) and accelerating them at the expense of some controlled deterioration of their output accuracy.

Another major research challenge is the extension of formal and experimental support for the known iterative refinement algorithms to new ones for LRA, in particular by combining our present algorithms and that of [14].

One can try to apply Alg. 4.1 to detection and verification of the gaps in the spectrum of the singular values of a matrix and to extend our whole approach to tensor decomposition.

Appendix

Appendix A Small families of matrices that are hard for superfast LRA

Any superfast algorithm fails to compute a close LRA of the following small families of matrices.

Example A.1**.**

Let $\Delta_{i,j}$ denote an $m\times n$ matrix of rank 1 filled with 0s except for its $(i,j)$th entry filled with 1. The $mn$ such matrices $\{\Delta_{i,j}\}_{i,j=1}^{m,n}$ form a family of $\delta$-matrices. We also include the $m\times n$ null matrix $O_{m,n}$ filled with 0s into this family. Now fix any superfast algorithm; it does not access the $(i,j)$th entry of its input matrices for some pair of $i$ and $j$. Therefore it outputs the same approximation of the matrices $\Delta_{i,j}$ and $O_{m,n}$, with an undetected error at least 1/2. Arrive at the same conclusion by applying the same argument to the set of $mn+1$ small-norm perturbations of the matrices of the above family and to the $mn+1$ sums of the latter matrices with any fixed $m\times n$ matrix of low rank. Furthermore, the same argument shows that superfast estimation of the norm of a matrix fails already on the same matrix family and that for both problems of LRA and norm estimation randomized superfast algorithms fail with error probability not close to 0.

Remark A.1**.**

How representative are the above matrix families of hard inputs? Clearly, they do not rule out superfast LRA computations if we restrict input matrices to be symmetric or diagonally dominant. Furthermore, the matrices of these families represent data singularities and are not relevant to a large class of input matrices representing regular processes or, say, smooth surfaces.

Appendix B Computation and approximation of $\rho$-top SVD of an LRA

B.1 Computation of $\rho$-top SVD

of an LRA

For completeness of our exposition we next recall the classical algorithm that computes $\rho$-top SVD of an LRA superfast (cf., e.g., [19]).

Algorithm B.1**.**

[Computation of $\rho$-top SVD of an LRA.]**

Input:

Four integers $\rho$, $k$, $m$, and $n$ such that $0<\rho\leq k\leq\min\{m,n\}$ and two matrices $A\in\mathbb{C}^{m\times k}$ and $B\in\mathbb{C}^{k\times n}$.

Output:

Three matrices $U\in\mathbb{C}^{m\times\rho}$ (unitary), $\Sigma\in\mathbb{R}^{\rho\times\rho}$ (diagonal), and $V\in\mathbb{C}^{n\times\rho}$ (unitary) such that $(AB)_{\rho}=U\Sigma V^{*}$ is a $\rho$-top SVD of $AB$.

Computations:

1. 1.

Compute SVDs $A=U_{A}\Sigma_{A}V^{*}_{A}$ and $B=U_{B}\Sigma_{B}V^{*}_{B}$ where $U_{B}\in\mathbb{C}^{m\times k}$, $V_{B}^{*}\in\mathbb{C}^{k\times n}$, and $\Sigma_{A},V^{*}_{A},U_{B},\Sigma_{B}\in\mathbb{C}^{k\times k}$. 2. 2.

Compute $k\times k$ matrices $W=\Sigma_{A}V^{*}_{A}U_{B}\Sigma_{B}$, $U_{W}$, $\Sigma_{W}$, and $V^{*}_{W}$ such that $W=U_{W}\Sigma_{W}V_{W}^{*}$ is SVD, $\Sigma_{W}=\operatorname{diag}(\Sigma,\Sigma^{\prime})$, and $\Sigma=\operatorname{diag}(\sigma_{j})_{j=1}^{\rho}$ and $\Sigma^{\prime}=\operatorname{diag}(\sigma_{j})_{j=\rho+1}^{k}$ are the matrices of the $\rho$ top (largest) and the $k-\rho$ trailing (smallest) singular values of the matrix $W$, respectively. Output the matrix $\Sigma$. 3. 3.

Compute and output the matrices $U$ and $V$ made up of the first $\rho$ columns of the matrices $U_{A}U_{W}$ and $V_{B}V_{W}$, respectively.

The algorithm involves $O((m+n)k^{2})$ flops; it is superfast if $k^{2}\ll n\leq m$. Its correctness follows from equations $AB=U_{A}WV^{*}_{B}$, $W=U_{W}\Sigma_{W}V^{*}_{W}$, and $\Sigma_{W}=\operatorname{diag}(\Sigma,\Sigma^{\prime})$.

For any matrix $M$ triangle inequality implies that

[TABLE]

where by following [10] we use unified notation $|||M|||$ for both spectral norm $||M||_{2}$ and Frobenius norm $||M||_{F}$ of a matrix $M$. Furthermore, $\tau_{\rho}(AB)\leq|||AB-M_{\rho}|||\leq|||AB-M|||+\tau_{\rho}(M)$, and so (cf. [37, Proposition 6.1])

[TABLE]

B.2 Approximation of $\rho$-top SVD of an LRA

We can readily extend Alg. B.1 to the computation of $\rho$-top SVD of a 3-factor LRA $M=XTY$, but next we quite closely approximate $\rho$-top SVD of an LRA $AB$ (and similarly of a 3-factor LRA) based on rank-revealing QRP factorization of the factors $A$ and $B$ rather than on their SVD.

Algorithm B.2**.**

[Approximation of $\rho$-top SVD of an LRA.]**

Input:

as in Alg. B.1.

Output:

Three matrices $U\in\mathbb{C}^{m\times\rho}$ (unitary), $\Sigma\in\mathbb{R}^{\rho\times\rho}$ (diagonal), and $V\in\mathbb{C}^{n\times\rho}$ (unitary) such that the matrix $U\Sigma V^{*}$ approximates a $\rho$-top SVD of $AB$.

Computations:

1. 1.

Compute Rank-revealing QRP and P*′LQ′** factorization $A=QRP$ and $B=P^{\prime}LQ;$ where $Q\in\mathbb{C}^{m\times k}$, $Q^{\prime}\in\mathbb{C}^{k\times n}$, $R$ and $L^{*}$ are $k\times k$ upper triangular matrices, and $P$ and $P^{\prime}$ are $k\times k$ permutation matrices.* 2. 2.

Define $\rho\times\rho$ leading principal submatrices $R_{\rho}$ of $R$ and $L^{\prime}_{\rho}$ of $L^{\prime}$. Define $\rho\times\rho$ matrix $P_{\rho}$ by first deleting the $k-\rho$ last rows of $P$ and then deleting the $k-\rho$ columns of the resulting matrix filled with 0s. Define $\rho\times\rho$ matrix $P_{\rho}^{\prime}$ by first deleting the $k-\rho$ last columns of $P^{\prime}$ and then deleting the $k-\rho$ rows of the resulting matrix filled with 0s. Compute the $\rho\times\rho$ nonsingular matrices $\bar{R}_{\rho}:=R_{\rho}P_{\rho}$, $\bar{L}^{\prime}_{\rho}:=P^{\prime}_{\rho}L^{\prime}_{\rho}$, and $A:=\bar{R}_{\rho}\bar{L}^{\prime}_{\rho}$. 3. 3.

Compute SVD $A=U_{A}\Sigma V^{*}_{A}$ where $U_{A}$, $\Sigma$, and $V^{*}_{A}$ are $\rho\times\rho$ matrices, $U_{A}$ and $V^{*}_{A}$ are unitary, and $\Sigma_{A}$ is a diagonal matrix filled with the singular values of $A$. Output the matrix $\Sigma$. 4. 4.

Compute and output the unitary matrices $U:=QU_{A}$ and $V:=V_{A}Q^{\prime}$.

The asymptotic complexity of the algorithm is $O((m+n)k^{2})$ flops, the same as for Alg. B.1, but the overhead constant in the ”O” notation is a little smaller because at stage 1 we substitute rank revealing QRP factorization for SVD. [6] computes QRP factorization that define LRA within the error bound $\sqrt{1+h^{2}(k-\rho)\rho}$ from optimal for any $h>1$, e.g., for $h=1.01$. Alternatively one can use rank-revealing $LUP$ factorization of [21] instead of $QRP$ factorization – staying within the same arithmetic cost bound but increasing the error bound to $1+h^{2}(k-\rho)\rho$.

Appendix C Generation of abridged SRHT matrices

In this section we specify the family of abridged SRHT matrices, used in Alg. 3.1 and defined by means of abridging the classical recursive processes of the generation of $n\times n$ SRHT matrices, obtained from the $n\times n$ dense matrices $H_{n}$ of Walsh-Hadamard transform for $n=2^{t}$ (cf. [15, Sec. 3.1]). The $n\times n$ matrices $H_{n}$ are obtained in $t=\log_{2}(n)$ recursive steps, but we only perform $d\ll t$ steps, and the resulting abridged matrix $H_{d,d}$ can be multiplied by a vector by using $2dn$ additions and subtractions. SRHT matrices are obtained from the matrices $H_{n}$ by means of random sampling and scaling, which we also apply to the $d$-abridged Hadamard transform matrices $H_{d,d}$. They turn into $H_{n}$ for $d=t$ but are sparse fo $d\ll t$. Namely, we write $H_{d,0}:=I_{n/2^{d}}$ and then specify the following recursive process:

[TABLE]

For any fixed pair of $d$ and $i$, each of the matrices $H_{d,i}$ is orthogonal up to scaling and has $2^{d}$ nonzero entries in every row and column; we can compute the product $MH$ for an $n\times k$ submatrix $H$ of $H_{d,d}$ by using less than $km2^{d}$ additions and subtractions.

Now define the $d$-Abridged Scaled and Permuted Hadamard matrices, $PDH_{d,d}$, where $P$ is a random sampling matrix and $D$ is the matrix of random integer diagonal scaling. Each random permutation or scaling contributes up to $n$ random parameters. We can involve more random parameters by applying random permutation and scaling also to some or all intermediate matrices $H_{d,i}$ for $i=0,1,\dots,d$.

The first $k$ columns of $H_{d,d}$ for $r\leq k\leq n$ form a $d$-Abridged Subsampled Randomized Hadamard Transform (SRHT) matrix $H$, which turns into a SRHT matrix for $d=t$, where $k=r+p$, $r$ is a target rank and $p$ is the oversampling parameter (cf. [10, Sec. 11]).