Low Rank Approximation at Sublinear Cost by Means of Subspace Sampling

1 Introduction

LRA Background. Low rank approximation (LRA) of a matrix is a hot research area of Numerical Linear Algebra (NLA) and Computer Science (CS) with applications to fundamental matrix and tensor computations and data mining and analysis (see surveys [HMT11], [M11], [KS16], and [CLO16]). Matrices defining Big Data (e.g., unfolding matrices of multidimensional tensors) are frequently so immense that realistically one can access and process only a tiny fraction of their entries, although quite typically these matrices admit their LRA, that is, are close to low rank matrices or equivalently have low numerical rank. One can operate with low rank matrices at sub-linear computational cost, that is, by using much fewer arithmetic operations and memory cells than an input matrix has entries, but can we compute LRA at sub-linear cost? Yes and no. No, because every sub-linear cost LRA algorithm fails even on the small input families of Appendix C. Yes, because our sub-linear cost variations of a popular subspace sampling algorithm output accurate LRA for a large class of input. Let us provide some details.

Subspace sampling algorithms compute LRA of a matrix $M$ by using auxiliary matrices $F M$ , $M H$ or $F M H$ for random multipliers $F$ and $H$ , commonly called test matrices and having smaller sizes. Their output LRA are nearly optimal whp provided that $F$ and/or $H$ are Gaussian, Rademacher’s, SRHT or SRFT matrices;¹¹1

Here and hereafter “SRHT and SRFT” are the acronyms for “Subsample Random Hadamard and Fourier transforms”;

“Gaussian” stands for “standard Gaussian (normal) random”; Rademacher’s are the matrices filled with iid variables, each equal to 1 or

- 1

with probability 1/2. furthermore the algorithms consistently output accurate LRA in their worldwide application with these and some other random multipliers

F

and

H

, all of which are multiplied by

M

at super-linear cost (see [TYUC17, Section 3.9], [HMT11, Section 7.4], and the bibliography therein).

Our modifications of these algorithms use sparse orthogonal (e.g., sub-permutation) multipliers²²2We define sub-permutation matrices as full-rank submatrices of permutation matrices. $F$ and $H$ , run at sub-linear cost, and as we prove, whp output reasonably accurate dual LRA

, that is, LRA of a random input admitting LRA; we deduce our error estimates under three distinct models of random matrix computations in Sections

4.1 – 4.3.

How meaningful is our result? Our definitions of three classes of random matrices of low numerical rank are quite natural for various real world applications of LRA, but are odd for some other ones. This, however, applies to any definition of that kind.

Our approach enables new insight into the subject, and our formal study is in good accordance with our numerical tests for both synthetic and real world inputs, some from [HMT11].

Our upper bounds on the output error of LRA of an $m \times n$ matrix of numerical rank $r$ exceed the optimal error bound by a factor of $\sqrt{min {m, n} r}$ , but if this optimal bound is small enough we can apply iterative refinement of LRA running at sub-linear cost (see [PLa]).

As we have pointed out, any sub-linear cost LRA algorithm (and ours are no exception) fails on some families of hard inputs, but our analysis and tests show that the class of such inputs is narrow. We conjecture that it shrinks fast if we recursively apply the same algorithm with new multipliers; we propose some heuristic recipes for these recursive processes, and our numerical tests confirm their efficiency.

Impact of our study, its extensions and by-products:

(i) Our duality approach enables new insight into some fundamental matrix computations besides LRA: [PQY15], [PZ17a], and [PZ17b] provide formal support for empirical efficiency of dual Gaussian elimination with no pivoting, while [PLb] proposes a sub-linear cost modification of Sarlós’ algorithm of 2006 and then proves that whp it outputs nearly optimal solution of the highly important problem of Linear Least Squares Regression (LLSR) provided that its input is random. Then again this formal proof is in good accordance with the test results.

(ii) In [PLSZa] we proved that popular Cross-Approximation LRA algorithms running at sub-linear cost as well as our simplified sub-linear cost variations of these algorithms output accurate solution of dual LRA whp, and we also devised a sub-linear cost algorithm for transformation of any LRA into its special form of CUR LRA, which is particularly memory efficient.

(iii) The paper [PLa] has proposed and elaborated upon sub-linear cost refinement of a crude but reasonably close LRA.

Related Works. Huge bibliography on LRA can be partly accessed via [M11], [HMT11], [KS16], [PLSZ17], [TYUC17], [OZ18], and the references therein. [PLSZ16] and [PLSZ17] were the first papers that provided formal support for dual accurate randomized LRA computations performed at sub-linear cost (in these papers such computations are called superfast). The earlier papers [PQY15], [PLSZ16], [PZ17a], and [PZ17b] studied duality for other fundamental matrix computations besides LRA, while the paper [PLb] has extended our study to a sub-linear cost dual algorithm for the popular problem of Linear Least Squares Regression and confirmed accuracy of this solution by the results of numerical experiments.

Organization of the paper. In Section 2 we recall random sampling for LRA. In Sections 3 and 4 we prove deterministic and randomized error bounds, respectively, for our dual LRA algorithms running at sub-linear cost. In Section 5 we discuss multiplicative pre-processing and generation of multipliers for both pre-processing and sampling. In Section 6 we cover our numerical tests. Appendix A is devoted to background on matrix computations. In Appendix B we prove our error bounds for dual LRA. In Appendix C we specify some small families of hard inputs for sub-linear cost LRA.

Some definitions. The concepts “large”, “small”, “ill-” and “well-conditioned”, “near”, “close”, and “approximate” are usually quantified in the context. “ $≪$ ” and “ $≫$ ” mean “much less than” and “much greater than”, respectively. “Flop”

stands for “floating point arithmetic operation”; “

iid” for “independent identically distributed”. In context a “perturbation of a matrix” can mean a perturbation having a small relative norm.

R^{p \times q}

denotes the class of

p \times q

real matrices, We assume dealing with real matrices throughout, and so the Hermitian transpose of

M

turns into transpose,

M^{*} = M^{T}

, but most of our study can be readily extended to complex matrices; see some relevant results about complex Gaussian matrices in [E88], [CD05], [ES05], and [TYUC17].

2 LRA by means of subspace sampling

2.1 Four subspace sampling algorithms

Hereafter $| | \cdot | |$ and $| | \cdot | |_{F}$ denote the spectral and the Frobenius matrix norms, respectively; $| \cdot |$ can denote either of them. $M^{+}$ denotes the Moore – Penrose pseudo inverse of $M$ .

Algorithm 2.1.

Column Subspace Sampling or Range Finder (see Remark 2.1).

Input:

An $m \times n$ matrix $M$ and a target rank $r$ .

Output:

Two matrices $X \in R^{m \times l}$ and $Y \in R^{l \times m}$ defining an LRA $~ M = X Y$ of $M$ .

Initialization:

Fix an integer $l$ , $r \leq l \leq n$ , and an $n \times l$ matrix $H$ of full rank $l$ .

Computations:

Compute the $m \times l$ matrix $M H$ .
Fix a nonsingular $l \times l$ matrix $T^{- 1}$ and output the $m \times l$ matrix $X := M H T^{- 1}$ .
Output an $l \times n$ matrix $Y := {a r g m i n}_{V} | X V - M |$ .

Remark 2.1.

Let $rank (M H) = l$ . Then $Y = (M H)^{+} M$ and $X Y = M H (M H)^{+} M$ independently of the choice of $T^{- 1}$ , but its proper choice numerically stabilizes the computations of the algorithm. For $l > r \geq nrank (M H)$ the matrix $M H$ is ill-conditioned,³³3 $nrank (W)$ denotes numerical rank of $W$ (see Appendix A.1). but $M H T^{- 1}$ is orthogonal for $T = R$ , $X := Q = M H R^{- 1}$ and $Y := Q^{*} M$ where $Q$ and $R$ are the factors of the thin QR factorization of $M H$ (cf. [HMT11, Algorithm 4.1]). It is also orthogonal for $T = R Π$ and the factors $R$ and $Π$ in a rank-revealing QR $Π$ factorization $M H = Q R Π$ .

Column Subspace Sampling turns into Column Subset Selection in the case of a sub-permutation matrix $H$ and turns into Row Subspace Sampling if it is applied to the transpose $M^{T}$ .

Algorithm 2.2.

Row Subspace Sampling or Transposed Range Finder. See Remark 2.2.

Input:

As in Algorithm 2.1.

Output:

Two matrices $X \in R^{k \times n}$ and $Y \in R^{m \times k}$ defining an LRA $~ M = Y X$ of $M$ .

Initialization:

Fix an integer $k$ , $r \leq k \leq m$ , and a $k \times m$ matrix $F$ of full numerical rank $k$ .

Computations:

Compute the $k \times m$ matrix $F M$ .
Fix a nonsingular $k \times k$ matrix $S^{- 1}$ ; then output $k \times n$ matrix $X := S^{- 1} F M$ .
Output an $m \times k$ matrix $Y := {a r g m i n}_{V} | V X - M |$ .

Row Subspace Sampling turns into Row Subset Selection in case of a sub-permutation matrix $F$ .

Remark 2.2.

Let $rank (F M) = l$ . Then $Y = M (F M)^{+}$ and $Y X = M (F M)^{+} F M$ independently of the choice of $S^{- 1}$ , but a proper choice of $S$ numerically stabilizes the computations of the algorithm. For $k > r \geq nrank (F M)$ the matrix $F M$ is ill-conditioned, but $S^{- 1} F M$ is orthogonal if $S = L$ , $X := Q = L^{- 1} F M$ , $Y := Q^{*} M$ , and $L$ and $Q$ are the factors of the thin LQ factorization of $F M$ or if $S = Π L$ and $Π$ and $L$ are the factors in a rank-revealing $Π$ LQ factorization $F M = Π L Q$ .

The following algorithm combines row and column subspace sampling. In the case of the identity matrix

S

it turns into the algorithm of [CW09, Theorems 4.7 and 4.8] and [TYUC17, Section 1.4], whose origin can be traced back to [WLRT08].

Algorithm 2.3.

Row and Column Subspace Sampling. See Remark 2.3.

Input:: As in Algorithm 2.1.
Output:: Two matrices $X \in R^{m \times k}$ and $Y \in R^{k \times m}$ defining an LRA $~ M = X Y$ of $M$ .
Initialization:: Fix two integers $k$ and $l$ , $r \leq k \leq m$ and $r \leq l \leq n$ ; fix two matrices $F \in R^{k \times m}$ and $H \in R^{n \times l}$ of full numerical ranks and two nonsingular matrices $S \in R^{k \times k}$ and $T \in R^{l \times l}$ .
Computations:: 1. Output the matrix $X = M H T^{- 1} \in R^{m \times l}$ .

2. Compute the matrices $V = S^{- 1} F M \in R^{k \times n}$ and $W = S^{- 1} F X \in R^{m \times l}$ .

3. Output the $l \times n$ matrix $Y := {a r g m i n}_{V} | W^{+} V - S^{- 1} F M |$ .

Remark 2.3.

If the matrix $F M H$ has full rank $min {k, n}$ , then $Y X = M H (F H M)^{+} F M,$ independently of the choice of the matrices $S^{- 1}$ and $T^{- 1}$ , but a proper choice of $S$ numerically stabilizes the computations of the algorithm. For $min {k, l} > r \geq nrank (F X)$ the matrix $F M$ is ill-conditioned, but $S^{- 1} F M$ is orthogonal if $S = L$ , $X := Q = L^{- 1} F M$ , $Y := Q^{*} M$ , and $L$ and $Q$ are the factors of the thin LQ factorization of $F M$ or if $S = Π L$ and $Π$ and $L$ are the factors in a rank-revealing $Π$ LQ factorization $F M = Π L Q$ .

Remark 2.4.

By applying Algorithm 2.3 to the transpose matrix $M^{*}$ we obtain Algorithm 2.4. It begins with column subspace sampling followed by row subspace sampling. We only study Algorithms 2.1 and 2.3 for input $M$ , but they turn into Algorithms 2.2 and 2.4 for the input $M^{*}$ .

2.2 The known error bounds

Theorem 2.1.

(i) Let $2 \leq r \leq l - 2$ and apply Algorithm 2.1 with a Gaussian multiplier $H$ . Then (cf. [HMT11, Theorems 10.5 and 10.6]) ⁴⁴4[HMT11, Theorems 10.7 and 10.8] estimate the norm $| M - X Y |$ in probability.

E | | M - X Y | |_{F}^{2} \leq (1 + \frac{r}{l - r - 1}) σ_{F, r + 1}^{2} (M),

E | | M - X Y | | \leq (1 + \sqrt{\frac{r}{l - r - 1}}) σ_{r + 1} (M) + \frac{e \sqrt{l}}{l - r} σ_{F, r + 1} (M) .

(ii) Let $4 [\sqrt{r} + \sqrt{8 log (r n)}]^{2} log (r) \leq l \leq n$ and apply Algorithm 2.1 with an SRHT or SRFT multiplier $H$ . Then (cf. [T11], [HMT11, Theorem 11.2])

| M - X Y | \leq \sqrt{1 + 7 n / l} {¯ σ}_{r + 1} (M) w i t h a p r o b a b i l i t y i n 1 - O (1 / r) .

Clarkson and Woodruff prove in [CW09] that Algorithm 2.3 reaches the bound ${¯ σ}_{r + 1} (M)$ within a factor of $1 + ϵ$ whp if the multipliers $F \in G^{k \times m}$ and $H \in G^{n \times l}$ are Rademacher’s matrices and if $k$ and $l$ are sufficiently large, having order of $r / ϵ$ and $r / ϵ^{2}$ for small $ϵ$ , respectively.

Tropp et al. argue in [TYUC17, Section 1.7.3] that LRA is not practical if the numbers $k$ and $l$ of row and column samples are large; iterative refinement of LRA at sub-linear cost in [PLa] can be a partial remedy. [TYUC17, Theorem 4.3] shows that the output LRA $X Y$ of Algorithm 2.3 applied with Gaussian multipliers $F$ and $H$ satisfies⁵⁵5In words, the expected output error norm $E | | M - X Y | |_{F}$ is within a factor of $(\frac{k l}{(k - l) (l - r)})^{1 / 2}$ from its minimum value $σ_{F, r + 1} (M)$ ; this factor is just 2 for $k = 2 l = 4 r$ .

E | | M - X Y | |_{F}^{2} \leq \frac{k l}{(k - l) (l - r)} σ_{F, r + 1}^{2} (M) i f k > l > r .

(2.1)

3 Deterministic output error bounds for sampling algorithms

3.1 Deterministic error bounds of Range Finder

Theorem 3.1.

[HMT11, Theorem 9.1]. Suppose that Algorithm 2.1 has been applied to a matrix $M$ with a multiplier $H$ and let

C_{1} = V_{1}^{*} H, C_{2} = V_{2}^{*} H,

(3.1)

M = (\begin{matrix} U_{1} & Σ_{1} & V_{1}^{*} U_{2} & Σ_{2} & V_{2}^{*} \end{matrix}), M_{r} = U_{1} Σ_{1} V_{1}^{*}, a n d M - M_{r} = U_{2} Σ_{2} V_{2}^{*}

(3.2)

be SVDs of the matrices $M$ , its rank- $r$ truncation $M_{r}$ , and $M - M_{r}$ , respectively. [ $Σ_{2} = O$ and $X Y = M$ if $rank (M) = r$ . The columns of $V_{1}^{*}$ span the top right singular space of $M$ .] Then

| M - X Y |^{2} \leq | Σ_{2} |^{2} + | Σ_{2} C_{2} C_{1}^{+} |^{2} .

(3.3)

Notice that $| Σ_{2} | = {¯ σ}_{r + 1} (M)$ , $| C_{2} | \leq 1$ , and $| Σ_{2} C_{2} C_{1}^{+} | \leq | Σ_{2} | | C_{2} | | C_{1}^{+} |$ and obtain

| M - X Y | \leq (1 + | C_{1}^{+} |^{2})^{1 / 2} {¯ σ}_{r + 1} (M) f o r C_{1} = V_{1}^{*} H .

(3.4)

It follows that the output LRA is optimal up to a factor of $(1 + | C_{1}^{+} |^{2})^{1 / 2}$ .

Next we deduce an upper bound on the norm $| C_{1}^{+} |$ in terms of $| | ((M H)_{r})^{+} | |$ , $| | M | |$ , and $η := 2 σ_{r + 1} (M) | | ((M H)_{r})^{+} | |$ . Given $M H$ we can compute the norm $| | ((M H)_{r})^{+} | |$ at sub-linear cost if $l ≪ m$ , and in some applications reasonable upper estimates for $| | M | |$ and $σ_{r + 1}$ are available.

Corollary 3.1.

Under the assumptions of Theorem 3.1 let the matrix $M_{r} H$ have full rank $r$ . Then

| (M_{r} H)^{+} | / | M_{r}^{+} | \leq | C_{1}^{+} | \leq | (M_{r} H)^{+} | | M_{r} | \leq | (M_{r} H)^{+} | | M | .

Proof.

Deduce from (3.1) and (3.2) that $M_{r} H = U_{1} Σ_{1} C_{1}$ . Hence $C_{1} = Σ_{1}^{- 1} U_{1}^{*} M_{r} H$ .

Recall that the matrix $M_{r} H$ has full rank $r$ , apply Lemma A.1, recall that $U_{1}$

is an orthogonal matrix, and obtain

| (M_{r} H)^{+} | / | Σ_{1}^{- 1} | \leq | C_{1}^{+} | \leq | (M_{r} H)^{+} | | Σ_{1} |

Substitute $| Σ_{1} | = | M_{r} |$ and $| Σ_{1}^{- 1} | = | M_{r}^{+} |$ and obtain the corollary. ∎

Corollary 3.2.

Under the assumptions of Corollary 3.1 let $η := 2 σ_{r + 1} (M) | | ((M H)_{r})^{+} | | < 1$ and $η^{'} := \frac{2 σ_{r + 1} (M)}{1 - η} | | ((M H)_{r})^{+} | | < 1.$ Then

\frac{1 - η^{'}}{| | M_{r}^{+} | |} | | ((M H)_{r})^{+} | | \leq | | C_{1}^{+} | | \leq \frac{| | M | |}{1 - η} | | ((M H)_{r})^{+} | | .

Proof.

Lemma A.3 implies that $max {| M_{r} H - M H |, | M H - (M H)_{r} |} \leq {¯ σ}_{r + 1} (M)$ .

Consequently $| M_{r} H - (M H)_{r} | \leq 2 {¯ σ}_{r + 1} (M)$ , and so $| | (M_{r} H)^{+} | | \leq \frac{1}{1 - η} | | ((M H)_{r})^{+} | |$ by virtue of Lemma A.2 if $η = 2 σ_{r + 1} (M) | | ((M H)_{r})^{+} | | < 1$ .

If in addition $η^{'} = \frac{2 σ_{r + 1} (M)}{1 - η} | | ((M H)_{r})^{+} | | < 1$ , then $2 σ_{r + 1} (M) | | (M_{r} H)^{+} | | < 1$ and therefore $| | ((M H)_{r})^{+} | | \leq \frac{1}{1 - η^{'}} | | (M_{r} H)^{+} | |$ by virtue of Lemma A.2.

Combine these bounds and obtain $(1 - η^{'}) | | ((M H)_{r})^{+} | | \leq | (M_{r} H)^{+} | \leq \frac{1}{1 - η} | | ((M H)_{r})^{+} | |$ .

Together with Corollary 3.1 this implies Corollary 3.2. ∎

For a given matrix $M H$ we can compute the norm $| | ((M H)_{r})^{+} | |$ at sub-linear cost if $l ≪ m$ . If also some reasonable upper bounds on $| | M | |$ and $σ_{r + 1} (M)$ are known, then Corollary 3.2 implies a posteriori estimates for the output errors of Algorithm 2.1.

3.2 Impact of pre-multiplication on the errors of LRA (deterministic estimates)

Lemma 3.1.

[The impact of pre-multiplication on LRA errors.] Suppose that Algorithm 2.3 outputs a matrix $X Y$ for $Y = (F X)^{+} F M$ and that $m \geq k \geq l = rank (X)$ . Then

M - X Y = W (M - X X^{+} M) f o r W = I_{m} - X (F X)^{+} F,

(3.5)

| M - X Y | \leq | W | | M - X X^{+} M |, | W | \leq | I_{m} | + | X | | F | | (X F)^{+} | .

(3.6)

Proof.

Recall that $Y = (F X)^{+} F M$ and notice that $(F X)^{+} F X = I_{l}$ if $k \geq l = rank (F X)$ . Therefore $Y = X^{+} M + (F X)^{+} F (M - X X^{+} M)$ . Consequently (3.5) and (3.6) hold. ∎

We bounded the norm $| M - X X^{+} M |$ in the previous subsection; next we bound the norms $| (F X)^{+} |$ and $| W |$ at sub-linear cost for $k l ≪ n$ , a fixed orthogonal $X$ , and proper choice of sparse $F$ .

Theorem 3.2.

[P00, Algorithm 1] for a real $h > 1$ applied to an $m \times l$ orthogonal matrix $X$ performs $O (m l^{2})$ flops and outputs an $l \times m$ sub-permutation matrix $F$ such that $| | (F X)^{+} | | \leq \sqrt{(m - l) l h^{2} + 1}$ , and $| | W | | \leq 1 + \sqrt{(m - l) l h^{2} + 1}$ , for $W = I_{m} + X (F X)^{+} F$ of (3.5) and any fixed $h > 1$ ; $| | W | | \approx \sqrt{m l}$ for $m ≫ l$ and $h \approx 1$ .

[P00, Algorithm 1] outputs $l \times m$ matrix $F$ . One can strengthen deterministic bounds on the norm $| W |$ by computing $k \times m$ sub-permutation matrices $F$ for $k$ of at least order $l^{2}$ .

Theorem 3.3.

For $k$ of at least order $l^{2}$ and a fixed orthogonal multiplier $X$ compute a $k \times m$ sub-permutation multiplier $F$ by means of deterministic algorithms by Osinsky, running at sub-linear cost and supporting [O18, equation (1)]. Then $| | W | | \leq 1 + | | (F X)^{+} | | = O (l)$ for $W$ of (3.5).

4 Accuracy of sub-linear cost dual LRA algorithms

Next we estimate the output errors of Algorithm 2.1 for a fixed orthogonal matrix $H$ and two classes of random inputs of low numerical rank, in particular for perturbed factor-Gaussian inputs of Definition A.1. These estimates formally support the observed accuracy of Range Finder with various dense multipliers (see [HMT11, Section 7.4], and the bibliography therein), but also with sparse multipliers, with which Algorithms 2.3 and 2.4 run at sub-linear cost. By applying the results of the previous section we extend these upper estimates for output accuracy to variations of Algorithm 2.3 that run at sub-linear cost; then we extend them to Algorithm 2.4 by means of transposition of an input matrix.

Our estimates involve the norms of a Gaussian matrix and its pseudo inverse (cf. Appendix A.4).

Definition 4.1.

A matrix is Gaussian if its entries are iid Gaussian variables. $G^{p \times q}$ is the class of $p \times q$ Gaussian matrices. $ν_{p, q} = | G |$ , $ν_{s p, p, q} = | | G | |$ , $ν_{F, p, q} = | | G | |_{F}$ , $ν_{p, q}^{+} = | G^{+} |$ , $ν_{s p, p, q}^{+} = | | G^{+} | |$ , and $ν_{F, p, q}^{+} = | | G^{+} | |_{F}$ for a matrix $G \in G^{p \times q}$ . [ $ν_{p, q} = ν_{q, p}$ and $ν_{p, q}^{+} = ν_{q, p}^{+}$ , for all pairs of $p$ and $q$ .]

Theorem 4.1.

[Non-degeneration of a Gaussian Matrix.] Let $F \in G^{r \times p}$ , $H \in G^{q \times r}$ , $M \in R^{p \times q}$ and $r \leq rank (M)$ . Then the matrices $F$ , $H$ , $F M$ , and $M H$ have full rank $r$ with probability 1.

Assumption 4.1.

We simplify the statements of our results by assuming that a Gaussian matrix has full rank and ignoring the probability 0 of its degeneration.

4.1 Output errors of Range Finder for a perturbed factor-Gaussian input

Theorem 4.2.

[Errors of Range Finder for a perturbed factor-Gaussian matrix.] Apply Algorithm 2.1 to a perturbation $M = ~ M + E$ of a right $m \times n$ factor-Gaussian matrix $~ M$ of rank $r$ such that

α := | | E | |_{F} / (σ_{r} (M) - σ_{r + 1} (M)) \leq 0.2 a n d ξ := 4 α ϕ < 1 f o r ϕ = ν_{s p, r, n} ν_{s p, r, l}^{+} | | H | |_{F} | | H^{+} | | .

(4.1)

(i) Then

| | M - X Y | |^{2} \leq (1 + (\frac{ϕ}{1 - ξ})^{2}) σ_{r + 1}^{2} (M) .

(ii) Let $ξ \leq 1 / 2$ with a probability close to 1 and let the integer $l - r$ be at least moderately large. Then with a probability close to 1

E | | M - X Y | | \leq θ σ_{r + 1} (M) f o r θ^{2} \approx 1 + ((2 e | | H | |_{F} | | H^{+} | | (\sqrt{n l} + \sqrt{r l}) / (l - r))^{2},

which is close to $1 + (2 e | | H | |_{F} | | H^{+} | |)^{2} n / l$ if $r ≪ l$ . Here and hereafter $e := 2.71828182 \dots$

Proof.

We prove claim (i) in Appendix B. Let us deduce claim (ii). Recall from Theorems A.4 and A.5

that the random variables

ν_{s p, r, n}

and

ν_{s p, r, l}^{+}

are strongly concentrated about their expected values

E (ν_{s p, r, n}) = \sqrt{n} + \sqrt{r}

and

E (ν_{s p, r, l}^{+}) = \frac{e \sqrt{l}}{l - r}

, respectively. Substitute these equations into the bound of claim (i), apply Jensen’s inequality, and deduce claim (ii) of the theorem. ∎

4.2 Output errors of Range Finder near a matrix with a random singular space

Next we prove similar estimates under an alternative randomization model for dual LRA.

Theorem 4.3.

[Errors of Range Finder for an input with a perturbed random singular space.] Let the matrix $V_{1}$ in Theorem 3.1 be the $n \times r$ Q factor in a QR factorization of a normalized $n \times r$ Gaussian matrix $G_{n, r}$ and let the multiplier $H$ be any $n \times l$ matrix of full rank $l \geq r$ .

(i) Then for $ν_{r, n}$ and $ν_{r, l}^{+}$ of Definition 4.1 it holds that

| M - X Y | / {¯ σ}_{r + 1} (M) \leq ϕ_{r, l, n} := (1 + (ν_{n, r} ν_{r, l}^{+} | H^{+} |)^{2})^{1 / 2} .

(ii) For a large or reasonably large integer $l - r$ , the random variable $ϕ_{r, l, n}$ is strongly concentrated about its expected values

E (ϕ_{s p, r, l, n}) = (1 + (\frac{\sqrt{n l} + \sqrt{r l}) e | | H^{+} | |}{l - r})^{2})^{1 / 2} a n d E (ϕ_{F, r, l, n}) = (1 + (\frac{n r^{2} | | H^{+} | |_{F}}{l - r - 1})^{2})^{1 / 2},

which turn into

E (ϕ_{s p, r, l, n}) \approx (1 + \frac{(e | | H^{+} | |)^{2} n}{l})^{1 / 2} a n d E (ϕ_{F, r, l, n}) \approx (1 + \frac{(e | | H^{+} | |)^{2} n}{l})^{1 / 2}

if $r ≪ l$ . Here $| | H^{+} | | = 1$ and $| | H^{+} | |_{F} = l$ if the matrix $H$ is orthogonal.

Proof.

Write

G_{n, r} = V_{1} R, V_{1} = G_{n, r} R^{- 1}, R = V_{1}^{*} G_{n, r},

(4.2)

and so $V_{1}^{*} = (R^{*})^{- 1} G_{r, n}$ and $V_{1}^{*} H = (R^{*})^{- 1} G_{r, n} H$ .

Let $H = U_{H} Σ_{H} V_{H}^{*}$ be SVD. Then $G_{r, n} U_{H} = G_{r, l}$ (cf. Lemma A.4), and hence

V_{1}^{*} H = (R^{*})^{- 1} G_{r, l} Σ_{H} V_{H}^{*} .

Therefore $| (V_{1}^{*} H)^{+} | \leq | R | ν_{r, l}^{+} | Σ_{H}^{- 1} |$ (apply Lemma A.1 and recall that $V_{H}$ is an orthogonal matrix).

Substitute $| Σ_{H}^{- 1} | = | H^{+} |$ and obtain $| (V_{1}^{*} H)^{+} | \leq ν_{r, l}^{+} | H^{+} | | R |$ .

Deduce from (4.2) that the matrices $R$ and $G_{n, r}$

share all their singular values. Therefore

| R | = ν_{n, r}

, and so

| (V_{1}^{*} H)^{+} | \leq ν_{n, r} ν_{r, l}^{+} | H^{+} |

By combining this bound with (3.4) prove claim (i) of the theorem.

Already for a reasonably large ratio $l / r$ the random variable $ϕ_{r, l, n}^{2} = 1 + (ν_{n, r} ν_{r, l}^{+} | H^{+} |)^{2}$ is strongly concentrated about its expected values

E (ϕ_{s p, r, l, n}^{2}) = 1 + (\frac{(\sqrt{n l} + \sqrt{r} l) e | | H^{+} | |}{l - r})^{2} a n d E (ϕ_{F, r, l, n}^{2}) = 1 + (\frac{n r^{2}}{l - r - 1} | | H^{+} | |_{F})^{2},

respectively (cf. Theorems A.4 and A.5), which are close to

1 + (e | | H^{+} | |)^{2} n / l a n d 1 + (n r^{2} | | H^{+} | |_{F} / l)^{2},

respectively, if $r ≪ l \leq n$ . Apply Jensen’s inequality and deduce claim (ii) of the theorem. ∎

Bound the output errors of Algorithms 2.3–2.4 by combining the estimates of this section and Section 3.2 and by transposing the input matrix $M$ .

4.3 Impact of pre-multiplication in the case of Gaussian noise

Next we deduce randomized estimates for the impact of pre-multiplication in the case where an input matrix $M$ includes considerable additive white Gaussian noise,⁶⁶6

Additive white Gaussian noise is statistical noise having a probability density function (PDF) equal to that of the Gaussian (normal) distribution. Additive white Gaussian noise is widely adopted in information theory and used in signal and image processing; in many cases it properly represents the errors of measurement and rounding (cf.

[SST06]). which is a classical representation of natural noise in information theory, is widely adopted in signal and image processing, and in many cases properly represents errors of measurement and rounding (cf. [SST06]).

Theorem 4.4.

Suppose that the multipliers $F$ and $H$ are orthogonal and let the input matrix $M$ be the sum of a fixed matrix $A$ and a scaled Gaussian matrix $E$ ,

M = A + E, \frac{1}{λ_{E}} E = G_{m, n} \in G^{m \times n}

(4.3)

for a constant $λ_{E}$ proportional to the norm $| | E | |$ . Then for any pair of orthogonal multipliers $F$ and $H$ it holds that

| W | \leq | I_{m} | + λ_{E} | M | min {ν_{n - l, l}^{+}, ν_{l, l}^{+}}, | W | \leq | I_{m} | + λ_{E} | M | ν_{n - l, l}^{+} f o r n \geq 2 l,

(4.4)

(4.5)

Proof.

Assumption (4.3) and Lemma A.4 together imply that $F E H$ is a scaled Gaussian matrix: $\frac{1}{λ_{E}} F E H \in G^{k \times l}$ . Hence $F M H = F A H + \frac{1}{λ_{E}} G_{k, l}$ . Apply Theorem A.3 and obtain that $| (F M H)^{+} | \leq λ_{E} min {ν_{n - l, l}^{+}, ν_{l, l}^{+}}$ . Recall from (3.5) that $| W | \leq | I_{m} | + | (F M H)^{+} | | M |$ since the multipliers $F$ and $H$ are orthogonal. Substitute the above bound on $| (F M H)^{+} |$ and obtain (4.4). Substitute equations $| | I_{m} | |_{F} = \sqrt{m}$ and $| | I_{m} | | = 1$ and claim (iii) of Theorem A.5 and obtain (4.5). ∎

Remark 4.1.

For $k = l = ρ$ , $S = T = I_{k}$ , sub-permutation matrices $F$ and $H$ , and a nonsingular matrix $F M H$ , Algorithms 2.3 and 2.4 output LRA in the form $C U R$ where $C \in R^{m \times ρ}$ and $R \in R^{ρ \times n}$ are two submatrices made up of $ρ$ columns and $ρ$ rows of $M$ and $U = (F M H)^{- 1}$ . [PLSZa] extends our current study to devising and analyzing algorithms for the computation of such CUR LRA in the case where $k$ and $l$ are arbitrary integers not exceeded by $ρ$ .

5 Multiplicative pre-processing and generation of multipliers

5.1 Multiplicative pre-processing for LRA

We proved that sub-linear cost variations of Algorithms 2.3 and 2.4 tend to output accurate LRA of random input whp. In the real world computations input matrices are not random, but we can randomize them by multiplying them by random matrices.

Algorithms 2.1 – 2.4 output accurate LRA whp if the multipliers are Gaussian, SRHT, SRFT or Rademacher’s (cf. [HMT11, Sections 10 and 11], [T11], [CW09]), but multiplication by these matrices run at super-linear cost. Our heuristic recipe is to apply these algorithms with a small variety of sparse multipliers $F_{i}$ and/or $H_{i}$ , $i = 1, 2, \dots$ , with which computational cost becomes sub-linear and then to monitor the accuracy of the output LRA by applying the criteria of the previous section, [PLa], and/or [PLSZa].

Various families of sparse multipliers have been proposed in [PLSZ16] and [PLSZ17]. One can readily complement these families with sub-permutation matrices and, say, sparse quasi Radmacher’s multipliers (see [PLSZa]) and then combine these basic multipliers together by using orthogonalized sums, products or other lower degree polynomials of these matrices as multipliers (cf. [HMT11, Remark 4.6]).

Next we specify a particular family of sparse multipliers, which was highly efficient in our tests when we applied them both themselves and in combination with other sparse multipliers.

5.2 Generation of abridged Hadamard and Fourier multipliers

We define multipliers of this family by means of abridging the classical recursive processes of the generation of $n \times n$ SRHT and SRFT matrices for $n = 2^{t}$ . These matrices are obtained from the $n \times n$ dense matrices $H_{n}$ of Walsh-Hadamard transform (cf. [M11, Section 3.1]) and $F_{n}$ of discrete Fourier transform (DFT) at $n$ points (cf. [P01, Section 2.3]), respectively. Recursive representation in $t$ recursive steps enables multiplication of the matrices $H_{n}$ and $F_{n}$

by a vector in

2 t n

additions and subtractions and

O (t n)

flops, respectively.

We end these processes in $d$ recursive steps for a fixed recursion depth $d$ , $1 \leq d \leq t$ , and obtain the $d$ -abridged Hadamard (AH) and Fourier (AF) matrices $H_{d, d}$ and $F_{d, d}$ , respectively, such that $H_{t, t} = H_{n}$ and $F_{t, t} = F_{n}$ . Namely write $H_{d, 0} = F_{d, 0} = I_{n / 2^{d}}$ , $i = \sqrt{- 1}$ , and $ω_{s} = exp (2 π i / s)$ denoting a primitive $s$ -th root of 1, and then specify two recursive processes:

H_{d, 0} = I_{n / 2^{d}}, H_{d, i + 1} = (\begin{matrix} H_{d, i} & H_{d, i} H_{d, i} & - H_{d, i} \end{matrix}) f o r i = 0, 1, \dots, d - 1,

(5.1)

F_{d, i + 1} = {ˆ P}_{i + 1} (\begin{matrix} F_{d, i} & F_{d, i} F_{d, i} {ˆ D}_{i + 1} & - F_{d, i} {ˆ D}_{i + 1} \end{matrix}), {ˆ D}_{i + 1} = diag (ω_{2^{i + 1}}^{j})_{j = 0}^{2^{i} - 1}, i = 0, 1, \dots, d - 1,

(5.2)

where ${ˆ P}_{i}$ denotes the $2^{i} \times 2^{i}$ matrix of odd/even permutations such that ${ˆ P}_{i} u = v$ , $u = (u_{j})_{j = 0}^{2^{i} - 1}$ , $v = (v_{j})_{j = 0}^{2^{i} - 1}$ , $v_{j} = u_{2 j}$ , $v_{j + 2^{i - 1}} = u_{2 j + 1}$ , $j = 0, 1, \dots, 2^{i - 1} - 1$ .⁷⁷7For $d = t$ this is a decimation in frequency (DIF) radix-2 representation of FFT. Transposition turns it into the decimation in time (DIT) radix-2 representation of FFT.

For any fixed pair of $d$ and $i$ , each of the matrices $H_{d, i}$ (resp. $F_{d, i}$ ) is orthogonal (resp. unitary) up to scaling and has $2^{d}$ nonzero entries in every row and column. Now make up multipliers $F$ and $H$ of $k \times m$ and $n \times l$ submatrices of $F_{d, d}$ and $H_{d, d}$ , respectively. Then in view of sparseness of $F_{d, d}$ or $H_{d, d}$ , we can compute the products $F M$ and $M H$ by using $O (k n 2^{d})$ and $O (l m 2^{d})$ flops, respectively, and they are just additions or subtractions in the case of submatrices of $H_{d, d}$ .

By combining random permutation with either Rademacher’s diagonal scaling for AH matrices $H_{d, d}$ or or random unitary diagonal scaling for AF matrices $F_{d, d}$ , we obtain the $d$ –Abridged Scaled and Permuted Hadamard (ASPH) matrices, $P D H_{n}$ , and $d$ –Abridged Scaled and Permuted Fourier (ASPF) $n \times n$ matrices, $P D F_{n}$ , where $P$ and $D$ are two matrices of permutation and diagonal scaling. Likewise define the families of ASH, ASF, APH, and APF matrices, $D H_{n, d}$ , $D F_{n, d}$ , $H_{n, d} P$ , and $F_{n, d} P$ , respectively. Each random permutation or scaling contributes up to $n$ random parameters. We can involve more random parameters by applying random permutation and scaling to the intermediate matrices $H_{d, i}$ and $F_{d, i}$ for $i = 0, 1, \dots, d$ .

Now the first $k$ rows for $r \leq k \leq n$ or first $l$ columns for $r \leq l \leq n$ of $H_{d, d}$ and $F_{d, d}$ form a $d$ -abridged Hadamard or Fourier multiplier, which turns into a SRHT or SRFT matrix, respectively, for $d = t$ . For $k$ and $l$ of order $r log (r)$ Algorithm 2.1 with a SRHT or SRFT multiplier outputs whp accurate LRA of any matrix $M$ admitting LRA (see [HMT11, Section 11]), but in our tests the output was consistently accurate even with sparse abridged SRHT or SRFT multipliers computed just in three recursive steps.

6 Numerical tests

In this section we cover our tests of dual sub-linear cost variants of Algorithm 2.1. The tests for Tables 6.1–6.4 have been performed by using MatLab on a Dell computer with the Intel Core 2 2.50 GHz processor and 4G memory running Windows 7; the standard normal distribution function randn of MATLAB has been applied in order to generate Gaussian matrices. The MATLAB function ”svd()” has been applied in order to calculate the $ξ$ -rank, i.e., the number of singular values exceeding $ξ$ for $ξ = 10^{- 5}$ in Sections 6.2 and 6.3 and $ξ = 10^{- 6}$ in Section 6.4. The tests for Tables 6.5–6.7 have been performed on a 64-bit Windows machine with an Intel i5 dual-core 1.70 GHz processor by using custom programmed software in $C^{+ +}$ and compiled with LAPACK version 3.6.0 libraries.

6.1 Input matrices for LRA

We generated the following classes of input matrices $M$ for testing LRA algorithms.

Class I: Perturbed $n \times n$ factor-Gaussian matrices with expected rank $r$ , $M = G_{1} * G_{2} + 10^{- 10} G_{3},$ for three Gaussian matrices $G_{1} \in G^{n \times r}$ , $G_{2} \in G^{r \times n}$ , and $G_{3} \in G^{n \times n}$ .

Class II: $M = U_{M} Σ_{M} V_{M}^{*}$ , for $U_{M}$ and $V_{M}$ being the Q factors of the thin QR orthogonalization of $n \times n$ Gaussian matrices and $Σ_{M} = diag (σ_{j})_{j = 1}^{n}$ ; $σ_{j} = 1 / j, j = 1, \dots, r$ , $σ_{j} = 10^{- 10}, j = r + 1, \dots, n$ (cf. [H02, Section 28.3]), and $n = 256, 512, 1024$ . (Hence $| | M | | = 1$ and $κ (M) = 10^{10}$ .)

Class III: (i) The matrices $M$ of the discretized single-layer Laplacian operator of [HMT11, Section 7.1]: $[S σ] (x) = c \int_{Γ_{1}} log | x - y | σ (y) d y, x \in Γ_{2}$ , for two circles $Γ_{1} = C (0, 1)$ and $Γ_{2} = C (0, 2)$ on the complex plane. We arrived at a matrix $M = (m_{i j})_{i, j = 1}^{n}$ , $m_{i, j} = c \int_{Γ_{1, j}} log | 2$

Low Rank Approximation at Sublinear Cost by Means of Subspace Sampling

Authors

Related Research

Refinement of Low Rank Approximation of a Matrix at Sub-linear Cost

Low Rank Approximation Directed by Leverage Scores and Computed at Sub-linear Cost

Low Rank Approximation of a Matrix at Sub-linear Cost

Rank Determination for Low-Rank Data Completion

Randomized Numerical Linear Algebra: Foundations Algorithms

Tensor-Krylov method for computing eigenvalues of parameter-dependent matrices

Theory meets Practice at the Median: a worst case comparison of relative error quantile algorithms

This week in AI

1 Introduction

2 LRA by means of subspace sampling

2.1 Four subspace sampling algorithms

Algorithm 2.1.

Remark 2.1.

Algorithm 2.2.

Remark 2.2.

Algorithm 2.3.

Remark 2.3.

Remark 2.4.

2.2 The known error bounds

Theorem 2.1.

3 Deterministic output error bounds for sampling algorithms

3.1 Deterministic error bounds of Range Finder

Theorem 3.1.

Corollary 3.1.

Proof.

Corollary 3.2.

Proof.

3.2 Impact of pre-multiplication on the errors of LRA (deterministic estimates)

Lemma 3.1.

Proof.

Theorem 3.2.

Theorem 3.3.

4 Accuracy of sub-linear cost dual LRA algorithms

Definition 4.1.

Theorem 4.1.

Assumption 4.1.

4.1 Output errors of Range Finder for a perturbed factor-Gaussian input

Theorem 4.2.

Proof.

4.2 Output errors of Range Finder near a matrix with a random singular space

Theorem 4.3.

Proof.

4.3 Impact of pre-multiplication in the case of Gaussian noise

Theorem 4.4.

Proof.

Remark 4.1.

5 Multiplicative pre-processing and generation of multipliers

5.1 Multiplicative pre-processing for LRA

5.2 Generation of abridged Hadamard and Fourier multipliers

6 Numerical tests

6.1 Input matrices for LRA