Low-rank approximation in the Frobenius norm by column and row subset selection

TL;DR

This paper introduces a fast, robust deterministic algorithm for low-rank matrix approximation using column and row subset selection, achieving error bounds close to the optimal Frobenius norm approximation, and extends the method to tensor approximations.

Contribution

It presents a novel, efficient deterministic algorithm for column and row subset selection with provable error bounds, improving upon previous randomized methods and extending to tensor approximations.

Findings

Algorithm guarantees error close to best Frobenius norm approximation.

Error bound within a factor of k+1 of optimal for cross approximation.

First polynomial-time deterministic algorithm with such error bounds.

Abstract

A CUR approximation of a matrix $A$ is a particular type of low-rank approximation $A \approx C U R$, where $C$ and $R$ consist of columns and rows of $A$, respectively. One way to obtain such an approximation is to apply column subset selection to $A$ and $A^T$. In this work, we describe a numerically robust and much faster variant of the column subset selection algorithm proposed by Deshpande and Rademacher, which guarantees an error close to the best approximation error in the Frobenius norm. For cross approximation, in which $U$ is required to be the inverse of a submatrix of $A$ described by the intersection of $C$ and $R$, we obtain a new algorithm with an error bound that stays within a factor $k + 1$ of the best rank-$k$ approximation error in the Frobenius norm. To the best of our knowledge, this is the first deterministic polynomial-time algorithm for which this factor is bounded by a polynomial in $k$. Our derivation and analysis of the algorithm is based on derandomizing a recent existence result by Zamarashkin and Osinsky. To illustrate the versatility of our new column subset selection algorithm, an extension to low multilinear rank approximations of tensors is provided as well.

Equations

$$\Vert A-A(:,I)A(:,I{)}^{+}A{\Vert }_F^2\le (k+1)\left({\sigma }_{k+1}(A{)}^2+\dots +{\sigma }_{\min \{m,n\}}(A{)}^2\right)$$

$$A\approx CUR,$$

$$U=C^{+}A{R}^{+}$$

$$A\approx A(:,J)A(I,J{)}^{-1}A(I,:),$$

$${\pi }_{i_1,\dots ,i_k}A:=A(:,I)\cdot A(:,I{)}^{+}\cdot A=Q{Q}^{T}A,$$

$$\mathbb{P}(X=I):=\frac{\mathrm{Vol}(A(:,I){)}^2}{{\sum }_{J\in \{1,\dots ,n{\}}^k}\mathrm{Vol}(A(:,J){)}^2}.$$

$$\mathbb{E}[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2]\le (k+1)\left({\sigma }_{k+1}^2+\dots +{\sigma }_m^2\right).$$

$$A=\left[\begin{array}{ccccc}a(1+\epsilon )& b& b& \dots & b\\ -b(1+\epsilon )& a& a& \dots & a\end{array}\right]$$

$$\mathbb{E}\left[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2\mid X_1=i_1\right]$$

$$\mathbb{E}\left[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2\mid X_1=i_1,\dots ,X_{t-1}=i_{t-1},X_t=i_t\right]$$

$$\mathbb{E}\left[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2\mid X_1=i_1,\dots ,X_t=i_t\right]=(k-t+1)\frac{c_{m-k+t-1}(B{B}^T)}{c_{m-k+t}(B{B}^T)},$$

$$(-\lambda {)}^m+c_{m-1}(-\lambda {)}^{m-1}+\dots +c_1(-\lambda )+c_0:=\det (B{B}^T-\lambda I).$$

$$\frac{c_{m-k+t-1}(B_i{B}_i^T)}{c_{m-k+t}(B_i{B}_i^T)}$$

$$B_i=B-{\pi }_{i}B=\left(I-\frac{b_i{b}_i^T}{\Vert b_i{\Vert }_2^2}\right)B$$

$$\det (\lambda I-B{B}^T)+\frac{1}{\Vert b_i{\Vert }_2^2}\sum _{j=1}^n{\sigma }_j^2(B)v_{ij}^2{g}_j(x).$$

$$A=\left[\begin{array}{cc}6.583644\cdot 10^{-7}& 8.113362\cdot 10^{-3}\\ 8.113362\cdot 10^{-3}& 100\end{array}\right],$$

$$B=U\Sigma V^T,U\in {\mathbb{R}}^{m\times m},\Sigma \in {\mathbb{R}}^{m\times m},V\in {\mathbb{R}}^{n\times m}.$$

$$U^T{B}_{i}V=(I-U^T{\pi }_{i}U)U^{T}BV=\left(I-\frac{U^T{b}_i{b}_i^{T}U}{\Vert b_i{\Vert }_2^2}\right)\Sigma =\left(I-q{q}^T\right)\Sigma ,$$

$$Q^T(I-q{q}^T)\Sigma W=(I-e_1{e}_1^T)(Q^T\Sigma W).$$

$$\mathbb{E}\left[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2\mid X_1=i_1,\dots ,X_{t-1}=i_{t-1}\right]\le (k+1)({\sigma }_{k+1}^2+\dots +{\sigma }_m^2)$$

$$\mathbb{E}\left[\Vert A-{\pi }_{X_1,\dots ,X_k}A{\Vert }_F^2\mid X_1=i_1,\dots ,X_t=i_t\right]\le (k+1)({\sigma }_{k+1}^2+\dots +{\sigma }_m^2)$$

$$\Vert A-a_1{a}_1^{+}A{\Vert }_F^2\le 2({\sigma }_2^2+\dots +{\sigma }_m^2).$$

$$A^{T}A=\left[\begin{array}{cc}1& b^T\\ b& A_2^T{A}_2\end{array}\right]=\left[\begin{array}{cc}1& b^T\\ b& B^{T}B+b{b}^T\end{array}\right]$$

$$\Vert A^{T}A{\Vert }_2\le {\Vert \left[\begin{array}{cc}1& \Vert b{\Vert }_2\\ \Vert b{\Vert }_2& \Vert B^{T}B+b{b}^T{\Vert }_2\end{array}\right]\Vert }_2\le {\Vert \left[\begin{array}{cc}1& \Vert b{\Vert }_2\\ \Vert b{\Vert }_2& \Vert B{\Vert }_F^2+\Vert b{\Vert }_2^2\end{array}\right]\Vert }_2.$$

$$\Vert A-a_1{a}_1^{+}A{\Vert }_F^2=\Vert B{\Vert }_F^2\le 2(\Vert A{\Vert }_F^2-\Vert A{\Vert }_2^2),$$

$$\Vert A{\Vert }_2^2\le 1+\Vert b{\Vert }_2^2+\frac{\Vert B{\Vert }_F^2}{2}=:\gamma .$$

$$p(\lambda )=\left(\lambda -1\right)\left(\lambda -\Vert b{\Vert }_2^2-\Vert B{\Vert }_F^2\right)-\Vert b{\Vert }_2^2.$$

$$p(\gamma )=\Vert B{\Vert }_F^2/2\cdot \left(1-\Vert b{\Vert }_2^2-\Vert B{\Vert }_F^2/2\right)\ge 0,$$

$$A=\left[\begin{array}{ccc}1& 0& 10^{-b}\\ 0& 1& 10^{-b}\\ 0& 0& 10^{-2b}\end{array}\right]$$

$$\Vert A-CUR{\Vert }_F\le \sqrt{2k+2}\sqrt{{\sigma }_{k+1}(A{)}^2+\cdots +{\sigma }_m(A{)}^2}.$$

Full text

1 Introduction

Given an $m\times n$ matrix $A$ and an integer $k$, typically much smaller than $m$ and $n$, the column subset selection problem aims at determining an index set $I\subset\{1,\ldots,n\}$ of cardinality $k$ such that the corresponding $k$ columns $A(:,I)$ represent a good approximation of the range of $A$. This problem has broad applications in a diversity of disciplines, including scientific computing, model reduction, and statistical data analysis. While column subset selection is a classical problem in numerical linear algebra, closely connected to rank-revealing QR factorizations [9, 10, 23], new significant theoretical and algorithmic developments have been achieved during the last two decades within the model reduction and theory of algorithms communities. In particular, this concerns the interplay between column subset selection and interpolation [3, 11, 18] as well as the development and analysis of randomized algorithms, see [6, 17, 19, 36] for a few references representing this research direction.

This paper is concerned with algorithmic improvements and extensions of the seminal work by Deshpande, Rademacher and co-authors [14, 15] on column subset selection. In [15] the existence of an index set $I$ such that

[TABLE]

has been established. Here, $\|\cdot\|_{F}$ and $(\cdot)^{+}$ denote the Frobenius norm and the Moore–Penrose inverse of a matrix, respectively. We let $\sigma_{1}(A)\geq\sigma_{2}(A)\geq\cdots$ denote the singular values of $A$. Note that $A(:,I)A(:,I)^{+}$ is an orthogonal projector and the bound (1) measures how well all the columns of $A$ are approximated by the subset of column contained in $I$. The bound (1) is remarkable because the singular value decomposition (SVD) of $A$ implies that the best approximation error $\|A-QQ^{+}A\|_{F}^{2}$ attained by an arbitrary $m\times k$ matrix $Q$ is given by $\sigma_{k+1}(A)^{2}+\ldots+\sigma_{\min\{m,n\}}(A)^{2}$. The bound (1) is larger by a factor that is only linear in $k$. More generally, we will call any quasi-optimal bound with a factor that is at most polynomial in $k$ (and independent of $m,n$ or $A$) a polynomial bound. The proof of (1) proceeds by defining a suitable discrete probability distribution on index tuples such that the expected value of the error with respect to this distribution satisfies the bound. This then implies the existence of at least one index set satisfying the bound as well. We remark that the factor $k+1$ in (1) cannot be improved [15, Proposition 3.3]. In [14], a deterministic algorithm has been developed by derandomizing this approach using the method of conditional expectations. These conditional expectations are given in terms of coefficients of certain characteristic polynomials and the algorithm from [14] attains efficiency by cheaply updating these coefficients. However, it is well known that working with characteristic polynomials in finite precision arithmetic is prone to massive numerical cancellation [29] and, as we will see, the algorithm from [14] is also affected by numerical instability. Our first contribution, presented in Section 2, consists of deriving a formulation of the algorithm that updates singular values instead of coefficients of characteristic polynomials. While our new variant enjoys the same favorable complexity, numerical experiments with matrices of different singular value decay indicate that it is numerically robust, achieving (1) up to the level of roundoff error. Based on a minor extension of the theory from [14, 15], we will also present a modification of the column selection strategy that results in significant speed ups of the algorithm.

In Section 3, we extend the developments from [14] to the problem of determining a rank-$k$ approximation of the form

[TABLE]

where $C=A(:,J)$ and $R=A(I,:)$ contain $k$ selected columns and rows of $A$, respectively. There is a simple and well established strategy to derive such an approximation, see, e.g., [17, 31]: One first applies column subset selection to $A$ and $A^{T}$ in order to determine $C$ and $R$, respectively. Given $C$ and $R$, the choice

[TABLE]

then minimizes the Frobenius norm error. We will show that this strategy combined with (1) results in an error that is at most a factor $\sqrt{2k+2}$ larger than the best rank-$k$ approximation error. While this is clearly a favorable bound, the choice (2) comes with a disadvantage. It involves the full matrix $A$, but sometimes only partial information on $A$ is available and has been used to determine $I,J$. One example for such a situation is the Chebfun2 construction for approximating bivariate functions [34], which uses a coarse discretization to cheaply determine $I,J$ and then evaluates the full matrix $A$ only along the cross containing the rows and columns determined by $I$ and $J$, respectively. The choice $U=A(I,J)^{-1}$ then leads to a rank-$k$ approximation of the form

[TABLE]

which is often called cross approximation. Choosing $I,J$ via column subset selection is not advisable in this setting; it may lead to (nearly) singular $A(I,J)$ and result in an unfavorable approximation error. On the other hand, Goreinov and Tyrtyshnikov [21] have established a polynomial bound for cross approximation in the maximum norm when choosing $I,J$ such that the volume of $A(I,J)$ is maximal. Recently, Zamarashkin and Osinsky [38] derived a polynomial bound in the Frobenius norm by extending the techniques from [15]. However, as far as we know, there is no polynomial time deterministic algorithm that guarantees a polynomial bound (in any norm); popular greedy algorithms lead to exponential bounds [12, 24] at best. One major contribution of this work is to derive such an algorithm via an extension of [14]; our algorithm guarantees a Frobenius norm error that is at most a factor $k+1$ larger than the best approximation error.

Section 4 contains an extension to tensors. In particular, we derive a deterministic algorithm that obtains a multilinear low-rank approximation that is constructed from the fibers of the tensor and satisfies a polynomial bound. Although our approach is a relatively straightforward extension of (2) and related approaches have been proposed in the literature [16, 22, 30], we are not aware that such an algorithm has been explicitly spelled out and analyzed.

2 Column subset selection

We start by providing more details on the approach from [14, 15] for the column subset selection problem. In the following we consider a matrix $A\in\mathbb{R}^{m\times n}$ with $m\leq n$ and rank at least $k$. We let $a_{i}$ denote the $i$th column of $A$ and $\pi_{i_{1},\ldots,i_{k}}A$ the orthogonal projection of $A$ on the subspace spanned by the columns $a_{i_{1}},\ldots,a_{i_{k}}$, that is,

[TABLE]

where $I=(i_{1},\ldots,i_{k})\in\{1,\ldots,n\}^{k}$ and $Q$ denotes an orthonormal basis of $A(:,I)$. Let us emphasize that $I$ is now a tuple. Although order is not important and we are ultimately interested in an index set, working with tuples simplifies the subsequent definition and manipulation of probability distributions. The volume of a rectangular matrix $B\in\mathbb{R}^{m\times k}$ with $k\leq m$ is defined as $\,\mathrm{Vol}(B):=\prod_{i=1}^{k}\sigma_{i}(B)^{2}$. Note that $\,\mathrm{Vol}(B)^{2}=\det(B^{T}B)$.

We now define a discrete probability distribution on integer tuples of the form $I\in\{1,\ldots,n\}^{k}$ corresponding to a selection of $k$ columns from $A$. For this purpose, let $X=(X_{1},\ldots,X_{k})$ be a $k$-tuple of random variables with values in $\{1,\ldots,n\}$ such that

[TABLE]

By convention, $\,\mathrm{Vol}\left(A(:,I)\right)=0$ whenever $i_{1},\ldots,i_{k}$ contain repeated indices. Then [15, Theorem 1.3] shows that

[TABLE]

In particular, this implies the existence of $I$ satisfying this bound.

In view of (3) and the prominent role played by maximum volume submatrices in low-rank approximation [21], it is tempting to expect that the $k$ columns of maximum volume satisfy (1). However, not only that it is NP hard to choose such columns [8], but they also fail to have this property. For instance, for $k=1$ consider the $2\times n$ matrix

[TABLE]

with $a^{2}+b^{2}=1$ and $\varepsilon>0$. The column of maximum volume (that is, of maximum Euclidean norm) is the first one. The approximation error obtained by this choice is given by $\|A-\pi_{1}A\|_{F}^{2}=n-1$, which is much larger than $2\sigma^{2}_{2}=2(1+\varepsilon)^{2}$ for $\varepsilon$ sufficiently small. Note that choosing any of the other columns yields the best approximation error $(1+\varepsilon)^{2}=\sigma_{2}^{2}$.

2.1 Algorithm by Deshpande and Rademacher

Deshpande and Rademacher [14] derived a deterministic algorithm for column subset selection by derandomizing (4) using the method of conditional expectations.

More specifically, the first step of the algorithm chooses an index $i_{1}$ such that

[TABLE]

is minimized. By construction, this quantity still satisfies the bound (4). More generally, having $t-1$ indices $i_{1},\ldots,i_{t-1}$ selected, step $t$ chooses an index $i_{t}$ such that

[TABLE]

is minimized. After $k$ steps we arrive at an index set $I$ of cardinality $k$ such that the desired bound (1) holds.

For the algorithm to be practical, it is crucial to compute the conditional expectations (5) efficiently. Lemma 21 in [14] shows that

[TABLE]

where the right-hand side involves the matrix $B=A-\pi_{i_{1},\ldots,i_{t}}A$ and coefficients $c_{j}\equiv c_{j}(BB^{T})$ of the characteristic polynomial

[TABLE]

It is therefore required to compute in every step for all values of $i$, the ratios

[TABLE]

where $B_{i}=A-\pi_{i_{1},\ldots,i_{t-1},i}A$.

In the following, we discuss the computation of (7) and show how the minimization problem (5) can be relaxed in order to accelerate the search for suitable indices.

2.2 Computation of characteristic polynomial coefficients

Assuming that the first $t-1$ indices have been selected, we set $B:=A-\pi_{i_{1},\ldots,i_{t-1}}A$. Then

[TABLE]

is a rank-1 modification of $B$. Deshpande and Rademacher [14] propose two methods to compute (7) for $i=1,\ldots,n$. In the following, we summarize them briefly.

Algorithm 2 in [14] computes $BB^{T}$ explicitly and then computes $B_{i}B_{i}^{T}$ as a rank-2 update of $BB^{T}$ for every $i=1,\ldots,n$. The characteristic polynomial of $B_{i}B_{i}^{T}$ is computed by establishing a similarity transformation to a matrix in Frobenius normal form [7, Section 16.6]. Fast matrix-matrix multiplication and inversion can be exploited so that the cost of this approach is $O(nm^{\omega}\log m)$, where $\omega\leq 2.373$ is the best exponent of matrix-matrix multiplication complexity. 2. 2.

Algorithm 3 in [14] computes the thin SVD of $B=U\Sigma V^{T}$, the characteristic polynomial of $BB^{T}$ from the squared singular values of $B$, and the auxiliary polynomials $g_{j}(x)=\prod_{\ell\neq j}\left(x-\sigma_{\ell}(B)^{2}\right)$ for $j=1,\ldots,m$. For $h=m-k+t$ and $h=m-k+t-1$, the coefficient $c_{h}(B_{i}B_{i}^{T})$ can then be computed as the coefficient of $x^{h}$ in

[TABLE]

The cost of this second approach is $O(m^{2}n)$.

The problem of computing the Frobenius normal form of a matrix is “numerically not viable” [28]. Also, updating directly the characteristic polynomial as in (8) is prone to numerical cancellation, leading to inaccurate results. For instance, consider the $2\times 2$ matrix

[TABLE]

and the column selection problem for $k=1$. Algorithm 4 in [14] using (8) selects the first column, giving an error $\|A-A(:,1)A(:,1)^{+}A\|_{F}\approx 1.2\cdot 10^{-6}\gg\sqrt{2}\sigma_{2}(A)=1.4\cdot 10^{-10}$.

Therefore, from now on we will avoid updating coefficients of characteristic polynomials and work with singular values instead. More specifically, we will compute the singular values of $B_{i}$ by updating the SVD of $B$ and then apply the Summation Algorithm [29, Algorithm 1] to compute the coefficients of the characteristic polynomial of $B_{i}B_{i}^{T}$ from its eigenvalues (that is, the squared singular values) with $O(m^{2})$ operations in a numerically forward stable manner. To describe the updating procedure, consider the (thin) SVD

[TABLE]

The (nonzero) singular values of $B_{i}$ and

[TABLE]

with $q=U^{T}b_{i}/\|b_{i}\|_{2}$, are identical. Using standard bulge chasing algorithms (see, e.g., [37, Algorithm 3.4] and [2]) it is possible to find orthogonal matrices $Q,W\in\mathbb{R}^{m\times m}$ such that $Q^{T}q=e_{1}$, where $e_{1}$ denotes the first unit vector, and $Q^{T}\Sigma W$ is upper bidiagonal. In turn, the singular values can be computed from the bidiagonal matrix

[TABLE]

The matrices $Q$ and $W$ are composed of $O(m^{2})$ Givens rotations [20, Section 5.1] and the computation of $Q^{T}\Sigma W$ requires to apply each of these rotations to at most $3$ vectors. In turn, the cost of computing this bidiagonal matrix is $O(m^{2})$, which is identical with the cost of computing its singular values [20, Section 8.6].

2.3 Overall algorithm

The described variation of the column subset selection algorithm by Deshpande and Rademacher is summarized in Algorithm 1. One execution of line 3 is $O(nm^{2})$, lines 6–9 are $O(m^{2})$, and lines 14–15 are $O(knm)$. In summary, the overall complexity of Algorithm 1 is $O(knm^{2})$. This is identical with the complexity of [14, Algorithm 4] combined with [14, Algorithm 3], and it is better than [14, Algorithm 4] combined with [14, Algorithm 2].

Note that instead of lines 14–15 we could have updated $B\leftarrow B-\pi_{i_{t}}B$. However, we noticed that recomputing $B$ in lines 18–19 tends to improve accuracy and it does not change the overall complexity.

2.4 Early stopping of column search

For each column index, Algorithm 1 needs to traverse $O(n)$ columns in order to find the one that minimizes the coefficient ratio or, equivalently, the conditional expectation. This column search can be shortened. To describe the idea, suppose that $i_{1},\ldots,i_{t-1}$ have already been selected such that

[TABLE]

holds. Now, we can choose any $i_{t}$ such that

[TABLE]

holds. The existence of $i_{t}$ is guaranteed by (9) but we do not need to find the one that minimizes the conditional expectation. It suffices to always choose in every step an index such that (10) is verified. By induction, the error bound (1) still holds.

The discussion above suggests to modify Algorithm 1 such that it computes

  $\,\mathrm{bound}=(k+1)\cdot\left(\sigma_{k+1}^{2}+\ldots+\sigma_{m}^{2}\right)$

in the beginning and substitute line 11 with

1:  **if ** $(k-t+1)\cdot\text{ratio}\leq\,\mathrm{bound}$ **then ** Set $i_{t}=i$ and break **end if **

To be able to stop the search early, it is important to test the columns in a suitable order. We found it beneficial to test the columns of $B$ in descending norm. For each step $t$, computing the norms of all columns of $B$ and sorting them has complexity $O(mn+n\log n)$. Although this choice is clearly heuristic, the following lemma provides some justification for it by showing that the column of largest norm is the right choice for $k=1$ provided that all other columns are sufficiently small.

Lemma 1**.**

Let $A=\begin{bmatrix}a_{1}&A_{2}\end{bmatrix}$. If $\|A_{2}\|_{F}\leq\|a_{1}\|_{2}$ then choosing the first column solves the column selection problem for $k=1$, that is,

[TABLE]

Note that the condition of the lemma is satisfied if the column norms of $A$ decay sufficiently fast, for instance if $\|a_{i}\|_{2}\leq\frac{\|a_{1}\|_{2}}{i}$ for $i=2,\ldots,n$.

Proof.

Without loss of generality we may assume that $\|a_{1}\|_{2}=1$. By setting $B=A_{2}-a_{1}a_{1}^{+}A_{2}=A_{2}-a_{1}a_{1}^{T}A_{2}$ and $b=A_{2}^{T}a_{1}$, we have

[TABLE]

and obtain

[TABLE]

Here, the first inequality is a norm-compression inequality [5, Section 9.10] and the second inequality follows from the fact that the involved matrices are positive.

We aim at proving

[TABLE]

which is equivalent to

[TABLE]

Thus, it remains to show that the larger eigenvalue of the symmetric positive definite $2\times 2$ matrix on the right-hand side of (11) is bounded by $\gamma$. For this purpose, we note that its characteristic polynomial is given by

[TABLE]

Setting $\gamma=1+\|b\|_{2}^{2}+\|B\|_{F}^{2}/2$, we obtain

[TABLE]

where we used that $\|b\|_{2}^{2}+\|B\|_{F}^{2}=\|A_{2}\|^{2}_{F}\leq\|a_{1}\|^{2}_{2}=1$. Because $p$ is a parabola with vertex $(1+\|b\|_{2}^{2}+\|B\|_{F}^{2})/2\leq\gamma$, it follows that the larger root of $p$ is bounded by $\gamma$, which completes the proof. ∎

It is important to not draw too many conclusions from Lemma 1. Consider, for example, the matrix

[TABLE]

for some integer $b$, say $b=16$. For $k=1$, the optimal choice is the third column, which is the one of smallest norm. This matrix also nicely illustrates that the obvious greedy approach (in order to get $k$ columns of $A$, one first chooses the best column, then the best column in the orthogonal complement, and so on) comes with no guarantees and may, in fact, utterly fail. For $k=2$ the optimal choice consists of the first two columns. On the other hand, the greedy approach for $k=2$ first selects the third column and then the first column, resulting in the arbitrarily bad error ratio $\frac{\text{error greedy}}{\text{error best}}\approx 10^{b}$.

2.5 Numerical experiments

Both variants of Algorithm 1, without and with early stopping, have been implemented in Matlab version R2019a. As the bulge chasing algorithm in line 7 would perform poorly in Matlab, this part has been implemented in C++ and is called via a MEX interface. All numerical experiments in this work have been run on an eight-core Intel Core i7-8650U 1.90 GHz CPU, 256 KB of level 2 Cache and 16 GB of RAM. Multithreading has been turned off in order to not distort the findings.

We have applied the algorithm to the following three matrices:

the Hilbert matrix $A_{1}\in\mathbb{R}^{200\times 200}$ given by $A_{1}(i,j)=\frac{1}{i+j-1}$; 2. 2.

$A_{2}\in\mathbb{R}^{100\times 200}$ given by $A_{2}(i,j)=\exp(-0.3\cdot|i-j|/200)$; 3. 3.

$A_{3}\in\mathbb{R}^{100\times 200}$ given by $A_{3}(i,j)=\left(\left(\frac{i}{200}\right)^{20}+\left(\frac{j}{200}\right)^{20}\right)^{1/20}$.

The obtained results are shown in Figures 1, 2, and 3 respectively. Each left plot contains, for different values of $k$, the approximation error $\|A_{i}-A_{i}(:,S)A_{i}(:,S)^{+}\|_{F}$ returned by Algorithm 1, without and with early stopping. We compare with the best rank-$k$ approximation error $\sqrt{\sigma_{k+1}^{2}+\ldots+\sigma_{m}^{2}}$ and the upper bound (1), that is, $\sqrt{(k+1)(\sigma_{k+1}^{2}+\ldots+\sigma_{m}^{2})}$. It can be seen that both variants of our algorithm stay below the upper bound, until it reaches the level of roundoff error. Interestingly, for the matrix $A_{2}$, which features the slowest singular value decay, the observed approximation error is much closer to the best approximation error than to the upper bound. The right plots of the figures show, for different values of $k$, the ratio between the execution times of Algorithm 1 without early stopping and with early stopping. For the variant with early stopping, we also plot the number of columns that were examined. In the most optimistic scenario, only $k$ columns need to be examined, which means that in every step of the algorithm already the first verifies the desired criterion. The plots reveal that our algorithm actually stays pretty close to this ideal situation, at least for the matrices considered. Note that for values of $k$ larger than the numerical rank of the matrix, Algorithm 1 starts computing ratios (6) from singular values of the order of machine precision. In turn, the computations are severely affected by roundoff error and it may, in fact, happen that the early stopping criterion is never satisfied. This leads to meaningless results and we therefore truncate the plots before this happens. A proper implementation of Algorithm 1 needs to detect such a situation and reduce $k$ accordingly.

3 Matrix approximation

In this section, we extend the developments from Section 2 on column subset selection to compute certain low-rank matrix approximations of a matrix $A\in\mathbb{R}^{m\times n}$ with $m\leq n$. As already discussed in the introduction, we will pursue two ways. First, in Section 3.1, we discuss a general CUR approximation obtained from applying column subset selection to the columns and rows of the matrix. Second, in Section 3.2, we present a novel approach to cross approximation, a specific type of CUR approximation, with guaranteed error bounds.

3.1 CUR approximation induced by column subset selection

Suppose that $C\in\mathbb{R}^{m\times k}$ and $R\in\mathbb{R}^{k\times n}$ have been chosen. Then the matrix $U\in\mathbb{R}^{k\times k}$ that minimizes $\|A-CUR\|_{F}$ is given by the projection $U=C^{+}AR^{+}$, see [33, p. 320]. The following corollary provides an error bound for the case when $C$ and $R$ are determined by the techniques from Section 2, leading to Algorithm 2. Closely related results can be found in the literature; see, for example, [17, Theorem 4], [30, Corollary 3.5], and [32, Theorem 4.1].

Corollary 2**.**

Let $A\in\mathbb{R}^{m\times n}$, with $1\leq k\leq m\leq n$. Then the CUR approximation returned by Algorithm 2 satisfies

[TABLE]

Proof.

Using the inequality (1) twice and the fact that $CC^{+}$ is an orthogonal projection, we obtain

[TABLE]

∎

3.1.1 Numerical experiments

We have tested a Matlab implementation of Algorithm 2 in the setting and for the matrices $A_{1},A_{2},A_{3}$ described in Section 2.5. Figure 4 displays the obtained approximation errors $\|A_{i}-CUR\|_{F}$ for different values of $k$. Again, we have tested both variants of Algorithm 1, without and with early stopping, within Algorithm 2. The speedups obtained from early stopping are very similar to the ones reported Section 2.5 and, therefore, we refrain from providing details.

We also consider, for $0<\alpha<1$, the $n\times n$ matrix

[TABLE]

where $Q\in\mathbb{R}^{n\times n}$ is determined as the orthogonal factor from the QR decomposition of

[TABLE]

This is known to be a challenging example for the CUR approximation induced by DEIM (discrete interpolation method); see [32, Section 4.2], which determines the row and column indices by greedily choosing a maximum volume submatrix of $U_{k}$ and $V_{k}$ containing the first $k$ left and right singular vectors of $A$, respectively. For the example above, the DEIM induced CUR approximation always chooses $1,\ldots,k$ for the column and row indices. For $\alpha=0.1$, $n=6$, $k=5$, the error resulting from this choice is given by

[TABLE]

which is a magnitude larger than the upper bound $\sqrt{2(k+1)}\sigma_{6}\approx 3.5\cdot 10^{-10}$ guaranteed by Algorithm 2. Note that the latter algorithm selects the last $5$ rows and columns for this example, leading to an error of $\approx 1.3\cdot 10^{-10}$.

3.2 Cross approximation

We now consider cross approximations, which take the form

[TABLE]

for row/column index tuples

[TABLE]

The different choice of the middle matrix makes a fundamental difference. In particular, as the following example shows, choosing the indices $I,J$ as in Algorithm 2 may lead to poor approximation error.

Example 3**.**

Consider $A=\begin{bmatrix}2\varepsilon&1\\ 1&\varepsilon\end{bmatrix}$ for $\varepsilon>0$ and $k=1$. Clearly, the first column and row satisfy the bound (1) for $k=1$ with respect to $A$ and $A^{T}$, respectively. However, the error of the corresponding cross approximation, $\|A-A(:,1)A(1,1)^{-1}A(1,:)\|_{F}=\frac{1}{2\varepsilon}-\varepsilon$, becomes arbitrarily large as $\varepsilon\to 0$.

Zamarashkin and Osinsky [38] have shown the existence of a cross approximation that satisfies a polynomial error bound in the Frobenius norm. To summarize their result, let

[TABLE]

be a $(2k)$-tuple of random variables with values in $\Omega$ such that

[TABLE]

Note that $\,\mathrm{Vol}\left(A(I,J)\right)=0$ whenever $i_{1},\ldots,i_{k}$ or $j_{1},\ldots,j_{k}$ contain repeated indices. Then [38, Theorem 1] shows that

[TABLE]

In particular, this implies that there exists $(I,J)\in\Omega$ such that

[TABLE]

In analogy to Section 2.1 and [14], we will now derandomize this result producing a polynomial-time deterministic algorithm that returns a cross approximation satisfying (15). The key for doing so is to find an expression for the conditional expectations that is easy to work with.

3.2.1 Conditional expectations

Lemma 4**.**

Let $1\leq t\leq k$ and $(i_{1},\ldots,i_{t},j_{1},\ldots,j_{t})$ be such that

[TABLE]

for a random $(2k)$-tuple $(X,Y)$ with the probability distribution defined by (13). Consider

[TABLE]

the remainder of cross approximation after choosing row indices $i_{1},\ldots,i_{t}$ and column indices $j_{1},\ldots,j_{t}$. Then

[TABLE]

with the coefficients $c_{m-k+t},c_{m-k+t-1}$ defined as in (6).

Proof.

To simplify notation, we let $I_{1}=(i_{1},\ldots,i_{t})$, $I_{2}=(i_{t+1},\ldots,i_{k})$, $I=(I_{1},I_{2})=(i_{1},\ldots,i_{k})$ and define $J_{1},J_{2},J$ analogously. In the following, we always use the convention that row and column summation indices range from $1$ to $m$ and from $1$ to $n$, respectively. We have that

[TABLE]

with

[TABLE]

For establishing the equality in (16) we used from [38, Lemma 1] that

[TABLE]

and, from (13), that

[TABLE]

We now aim at simplifying the expression (16). For this purpose, we assume without loss of generality that $i_{1}=1,\ldots,i_{t}=t$ and $j_{1}=1,\ldots,j_{t}=t$. This allows us to partition

[TABLE]

where $B(I_{2},J_{2})=A(I_{2},J_{2})-A(I_{2},J_{1})A(I_{1},J_{1})^{-1}A(I_{1},J_{2})$ by the definition of $B$. By the relation between determinants and Schur complements [25, Equation (0.8.5.1)], $\,\mathrm{Vol}(A(I,J))=\,\mathrm{Vol}(A(I_{1},J_{1}))\cdot\,\mathrm{Vol}(B(I_{2},J_{2}))$. Therefore,

[TABLE]

Analogously, one shows

[TABLE]

Inserting these expressions into (16) yields

[TABLE]

By [26, Theorem 7] this ratio is equal to

[TABLE]

3.2.2 Derandomized cross approximation algorithm

With Lemma 4 at hand, we can proceed analogously to Section 2.1 and sequentially find $k$ pairs of row/column indices such that (15) is satisfied. Suppose that $t-1$ index pairs $(i_{1},j_{1}),\ldots,(i_{t-1},j_{t-1})$ have been determined. Then the $t$th step of the algorithm proceeds by choosing $(i_{t},j_{t})$ such that

[TABLE]

is minimized. We will show in Theorem 5 below that this choice of index pairs leads to a cross approximation satisfying the desired error bound (15). In view of Lemma 4, the minimization of (17) means that in each step of the algorithm we need to compute the ratios

[TABLE]

where

[TABLE]

Parallelizing the developments in Section 2.2, we now show how the coefficients in (18) can be computed via updating the singular values of $C_{ij}$. Let us denote the remainder from the previous step by

[TABLE]

Then it follows that

[TABLE]

see, e.g., [4]. We compute a thin SVD $B=U\Sigma V^{T}$ such that $U\in\mathbb{R}^{m\times m}$, $V\in\mathbb{R}^{m\times n}$ have orthonormal columns and $\Sigma\in\mathbb{R}^{m\times m}$ is diagonal. Note that

[TABLE]

Inserted into (19), this shows that the nonzero singular values of $C_{ij}$ match the singular values of

[TABLE]

where $x=\Sigma V(j,:)^{T}$ and $y=\frac{1}{B(i,j)}\Sigma U(i,:)^{T}$ are vectors of length $m$ and can be computed with $O(m^{2})$ operations.

Similarly as in Section 2.2, we transform $\Sigma-xy^{T}$ into bidiagonal form, after which its singular values can be computed with $O(m^{2})$ operations. This transformation proceeds in three steps:

We compute orthogonal matrices $Q$ and $W$ such that $Q^{T}\Sigma W$ is upper bidiagonal and $Q^{T}x=\pm\|x\|_{2}\cdot e_{1}$ using, for example, [37, Algorithm 3.4]. In turn, the matrix

[TABLE]

is bidiagonal with an additional nonzero first row; see the first plot in Figure 5 for an illustration. 2. 2.

By a bulge chasing algorithm, we transform $D_{1}$ to an upper banded matrix $D_{2}$ with two superdiagonals using $O(m^{2})$ Givens rotations. We refrain from giving a detailed description of the algorithm and refer to Figure 5 for an illustration. 3. 3.

The banded matrix $D_{2}$ is reduced to a bidiagonal matrix $D_{3}$ using the LAPACK [1] routine dgbbrd.

The overall procedure described above can be implemended by means of $O(m^{2})$ Givens rotations, each of which is applied to a small matrix of size independent of $m,n$. Hence, it has complexity $O(m^{2})$.

Algorithm 3 summarizes our newly proposed method for cross approximation. The SVD needed at the beginning of each outer loop is of complexity $O(m^{2}n)$ and each of the $mn$ inner loops costs $O(m^{2})$ operation; the total complexity of Algorithm 3 is therefore $O(km^{3}n)$.

Theorem 5**.**

For a matrix $A$ of rank at least $k$, Algorithm 3 returns index sets $I$ and $J$ such that (15) is satisfied.

Proof.

Let $B_{\{X,Y\}}=A-A(:,Y)A(X,Y)^{-1}A(X,:)$. For $t=1,\ldots,k$ we have that

[TABLE]

Therefore, as (14) holds, the choice (17) inductively ensures that

[TABLE]

Therefore, the index sets $I$ and $J$ computed by Algorithm 3 satisfy the bound (15). ∎

In analogy to the discussion in Section 2.4, let us emphasize that it is not necessary to select the pair $(i_{t},j_{t})$ that minimizes the ratio $r$. Any pair $(i,j)$ for which the inequality

[TABLE]

holds will lead to index sets $I$ and $J$ such that (15) is satisfied. Inspired by adaptive cross approximation with full pivoting [4], we traverse the entries of $B$ from the largest to the smallest (in magnitude) and stop the search once we have found an index pair $(i_{t},j_{t}$) satisfying (21).

3.2.3 Numerical experiments

We have implemented both variants of Algorithm 3, without and with early stopping, in Matlab. Again, the two inner loops have been implemented in a C++ function that is called via a MEX interface. The computational environment is the one described in Section 2.5 but the test matrices are smaller because Algorithm 3 without early stopping is significantly slower. We choose $A_{1}$ to be $100\times 100$, $A_{2}$ to be $50\times 100$, and the matrix $A_{3}\in\mathbb{R}^{50\times 100}$ is given by $A_{3}(i,j)=\left(\left(\frac{i}{100}\right)^{10}+\left(\frac{j}{100}\right)^{10}\right)^{1/10}$.

The left plots of Figures 6, 7, 8 display the approximation error $\|A-A(:,J)A(I,J)^{-1}A(I,:)\|_{F}$ for the index sets returned by both variants of Algorithm 3. It can be seen that the approximation errors often stay close to the best rank-$k$ approximation error $\sqrt{\sigma_{k+1}^{2}+\ldots+\sigma_{m}^{2}}$ and do not exceed the upper bound (15), modulo roundoff error. However, for larger values of $k$, Algorithm 3 without early stopping appears to encounter stability issues; the approximation error is distorted well above the level of roundoff error. This appears to be due to the fact that $A(I,J)$ becomes almost singular. For instance, Algorithm 3 without early stopping applied to $A_{2}$ with $k=48$ yields a matrix $A(I,J)$ with condition number $\approx 1.3\cdot 10^{18}$. The variant with early stopping appears to lead to lower condition numbers and does not exhibit numerical instability for the matrices considered. The right plots of the figures display the ratios between the execution time of Algorithm 3 without and with early stopping, as well as the total number of index pairs that needed to be tested in Algorithm 3 with early stopping. It can be observed that early stopping dramatically accelerates the computation and is thus the preferred variant.

We also consider the $n\times n$ matrix $A=LDL^{T}$, where

[TABLE]

with $s=\sin(\theta),c=\cos(\theta)$ for some $0<\theta<\pi$. This is known to be a challenging example for greedy cross approximation [24]: When $k=n-1$ the greedy algorithm selects the leading $k\times k$ submatrix and returns an approximation error that is exponentially larger than the best approximation error. In contrast, Algorithm 3, with and without early stopping, makes the correct choice by selecting the last $n-1$ rows and columns. For instance, for $n=6$ and $\theta=0.1$, we obtain the error

[TABLE]

Selecting the first $5$ rows and columns results in an error of $9.8\cdot 10^{-11}$.

Finally, we would like to point out an interesting observation concerning the preservation of structure. In joint work with Massei [12], we have shown that for a symmetric positive definite matrix $A$ there is always a symmetric choice of indices, $J=I$, leading to a symmetric cross approximation such that the favorable error bound of Goreinov and Tyrtyshnikov [21] is attained. For cross approximation in the Frobenius norm, the situation appears to be more complicated; it is generally not true that a symmetric choice of indices achieves the error bound (15) even when $A$ is symmetric positive definite. For instance, for $n=3$ and $k=1$ consider

[TABLE]

The best symmetric choice is $I=J=\{3\}$ but this leads to an error $\approx 0.1911>2\sqrt{\sigma_{2}^{2}+\sigma_{3}^{2}}\approx 0.1821$.

4 Tensor approximation

As shown, e.g., in [16, 22, 30], column subset selection can be used to approximate tensors as well. In the following, we demonstrate the use of the algorithm from Section 2 for to obtain approximations of low multilinear rank constructed from the fibers of a third-order tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$.

First, we briefly recall some basic definitions for tensors and refer to [27] for more details. Generalizing the notion of rows and columns of a matrix, the vectors obtained from $\mathcal{A}$ by fixing all indices but the $\mu$th one are called $\mu$-mode fibers. The matrix $A^{(\mu)}\in\mathbb{R}^{n_{\mu}\times(n_{1}n_{2}n_{3})/n_{\mu}}$ containing all $\mu$-mode fibers as columns is called the $\mu$-mode matricization of $\mathcal{A}$. The $\mu$-mode product of a matrix $B\in\mathbb{R}^{m\times n_{\mu}}$ with $\mathcal{A}$ is denoted by $B\times_{\mu}\mathcal{A}$ and it is the tensor such that its $\mu$-mode matricization is given by $B\cdot A^{(\mu)}$. We use the Frobenius norm of a tensor defined by

[TABLE]

and recall that $\|\mathcal{A}\|_{F}=\|A^{(\mu)}\|_{F}$ for $\mu=1,2,3$. The tuple $(k_{1},k_{2},k_{3})$ defined by $k_{\mu}=\,\mathrm{rank}(A^{(\mu)})$ is called the multilinear rank of $\mathcal{A}$ and we can decompose $\mathcal{A}$ as

[TABLE]

for coefficient matrices $B_{\mu}\in\mathbb{R}^{n_{\mu}\times k_{\mu}}$ for $\mu=1,2,3$ and a so called core tensor $\mathcal{C}\in\mathbb{R}^{k_{1}\times k_{2}\times k_{3}}$. This so called Tucker decomposition is particularly beneficial when the multilinear rank is much smaller than the size of a tensor.

Algorithm 4 produces an approximate Tucker decomposition for a given tensor such each coefficient matrix $B_{\mu}$ is composed of $\mu$-mode fibers. The following result shows that the obtained approximation error remains close to the best approximation error.

Corollary 6**.**

Consider $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and integers $k_{1},k_{2},k_{3}$ such that $1\leq k_{\mu}\leq n_{\mu}$ for $\mu=1,2,3$. Then the output of Algorithm 4 satisfies

[TABLE]

where $\mathcal{A}_{\mathrm{best}}$ is the best Tucker approximation of $\mathcal{A}$ of multilinear rank at most $(k_{1},k_{2},k_{3})$.

Proof.

The proof is similar to existing proofs on the quasi-optimality of the Higher-Order SVD [13] and related results in [16, 22, 30].

Using (1) and setting $\pi_{\mu}=B_{\mu}B_{\mu}^{+}$, the result of Algorithm 1 applied to $A^{(\mu)}$ satisfies

[TABLE]

where the second inequality follows from the fact that $A_{\mathrm{best}}^{(\mu)}$, the $\mu$-mode matricization of $\mathcal{A}_{\mathrm{best}}$, has rank at most $k_{\mu}$. Using the orthogonality of the projections $\pi_{\mu}$, we obtain

[TABLE]

where the second equality follows from [35, Theorem 5.1]. ∎

Remark 7**.**

Algorithm 4 easily generalizes to tensors of arbitrary order. Given a tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times\ldots\times n_{d}}$ and integers $k_{1},\ldots,k_{d}$, this generalization constructs subsets of fibers $B_{1},\ldots,B_{d}$ and a core tensor $\mathcal{C}$ such that

[TABLE]

4.1 Numerical experiments

We have implemented Algorithm 4 in Matlab and tested it on two $50\times 50\times 50$ tensors, given by $\mathcal{A}_{1}(i,j,h)=\frac{1}{i+j+h-1}$ and $\mathcal{A}_{2}(i,j,h)=\left(i^{10}+j^{10}+h^{10}\right)^{1/10}/50$. We choose $k_{1}=k_{2}=k_{3}=k$ and report in Figure 9 the obtained approximation errors $\|\mathcal{A}_{i}-B_{1}\times B_{2}\times B_{3}\times\mathcal{C}\|_{F}$ for different values of $k$, where $B_{1},B_{2},B_{3}$, $\mathcal{C}$ are returned by Algorithm 4, with and without early stopping in the column selection part. We compare with the quantity

[TABLE]

which provides a (tight) upper bound on the best approximation error. It can be seen that the errors obtained from Algorithm 4 remain close to this quasi-best approximation error.

5 Conclusions

In this work, we have proposed several improvements to the column selection algorithm by Deshpande and Rademacher [14]. The numerical experiments indicate that updating singular values (instead of characteristic polynomials) leads to numerical robustness, in the sense that the approximation error obtained in finite precision arithmetic is not affected unduly by roundoff error. We have also developed an extension of [14] to produce cross approximations of matrices and, to the best of our knowledge, this extension constitutes the first deterministic polynomial time algorithm that yields a cross approximation with a guaranteed polynomial error bound. We have introduced a mechanism for stopping early the search for indices in column subset selection or cross approximation. Although relatively simple, this mechanism tremendously reduces the execution time for all examples tested.

A number of issues remain for future study, such as the numerical stability analysis of our algorithms. In particular, it would be desirable to study the numerical robustness of the cross approximation returned by Algorithm 3 with early stopping. Also, by combining early stopping with a more aggressive reuse of the SVD might lead to further complexity reduction, but a rigorous complexity analysis would require deeper understanding of early stopping, well beyond the limited scope of Lemma 1. Finally, we would like to stress that the algorithms presented in this work are intented for small to medium sized matrices and tensors. For large-scale data, the algorithms presented in this paper need to be combined with other, possibly heuristic dimensionality reduction techniques.

Acknowledgements.

The authors thank Sergey Dolgov for helpful discussions on topics related to the work presented in this paper.