On maximum volume submatrices and cross approximation for symmetric semidefinite and diagonally dominant matrices

TL;DR

This paper investigates maximum volume submatrices in symmetric semidefinite and diagonally dominant matrices, providing bounds on approximation errors and enhancing cross approximation methods.

Contribution

It proves that maximum volume submatrices can be principal in certain matrices and extends error bounds for cross approximation methods.

Findings

Maximum volume submatrices can be principal in symmetric semidefinite and diagonally dominant matrices.

Error bounds for cross approximation are extended to more general matrices.

For doubly diagonally dominant matrices, the approximation error is within a modest factor of the best.

Abstract

The problem of finding a $k \times k$ submatrix of maximum volume of a matrix $A$ is of interest in a variety of applications. For example, it yields a quasi-best low-rank approximation constructed from the rows and columns of $A$. We show that such a submatrix can always be chosen to be a principal submatrix if $A$ is symmetric semidefinite or diagonally dominant. Then we analyze the low-rank approximation error returned by a greedy method for volume maximization, cross approximation with complete pivoting. Our bound for general matrices extends an existing result for symmetric semidefinite matrices and yields new error estimates for diagonally dominant matrices. In particular, for doubly diagonally dominant matrices the error is shown to remain within a modest factor of the best approximation error. We also illustrate how the application of our results to cross approximation for functions leads to new and better convergence results.

Equations

$$\Vert A-A(:,J)A(I,J{)}^{-1}A(I,:){\Vert }_{\max }\le (k+1){\sigma }_{k+1}(A),$$

$${\sigma }_i(A(I,I))={\sigma }_i\left(C(:,I{)}^{\ast }C(:,I)\right)={\sigma }_i(C(:,I){)}^2.$$

$$\det (A(I,J){)}^2$$

$$\le {\left({\Pi }_{i=1}^k{\sigma }_i(C(:,I)){\Pi }_{j=1}^k{\sigma }_j(C(:,J))\right)}^2$$

$$=\left({\Pi }_{i=1}^k{\sigma }_i(C(:,I{)}^{\ast }C(:,I))\right)\left({\Pi }_{j=1}^k{\sigma }_j(C(:,J{)}^{\ast }C(:,J))\right)$$

$$=\left({\Pi }_{i=1}^k{\sigma }_i(A(I,I))\right)\left({\Pi }_{j=1}^k{\sigma }_j(A(J,J))\right)$$

$$=\det (A(I,I))\cdot \det (A(J,J)),$$

$$\sum _{j=1,j\ne i}^n|a_{ij}|\le |a_{ii}|,i=1,\dots ,n.$$

$$\det (T(I,J))$$

$$\det (T(I,I))$$

$$|b_{i,1}|+|b_{i,2}|+\dots +|b_{i,k+1}|<1.$$

$$|\det(B)|=\Big{|}\sum_{j=1}^{k+1}(-1)^{i+j}b_{ij}\det(B_{ij})\Big{|}\leq\sum_{j=1}^{k+1}|b_{ij}|\,|\det(B_{ij})|\leq\sum_{j=1}^{k+1}|b_{ij}|<1.$$

$$|\det (A(I,J))|\le |\det (A(I,I))|.$$

$$A(I,I)=L(I,I)U(I,I),A(I,J)=L(I,I)U(I,J).$$

$$|\det (A(I,I))|=|\det (U(I,I))|,|\det (A(I,J)|=|\det (U(I,J))|.$$

$$A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right],A_{11}\in {\mathbb{R}}^{k\times k}.$$

$$A$$

$$=$$

$$A^{(k)}:=A_{22}-A_{21}{A}_{11}^{-1}{A}_{12}.$$

$$(i_{k+1},j_{k+1})=\operatorname*{arg\,max}\big{\{}\big{|}\det\!\big{(}A(I\cup\{i\},J\cup\{j\})\big{)}\big{|}:\,i\not\in I,j\not\in J\big{\}}.$$

$$\det\big{(}A(\tilde{I},\tilde{J})\big{)}=\det\big{(}A(I,J)\big{)}\cdot A^{(k)}(\tilde{i},\tilde{j}).$$

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

$$\rho_{k}:=\sup_{A}\big{\{}\|A^{(k)}\|_{\max}/\|A\|_{\max}\big{\}},$$

$${\rho }_k\le \sqrt{k+1}\cdot \sqrt{2\cdot 3^{1/2}\cdot 4^{1/3}\cdot \dots \cdot (k+1{)}^{1/k}}\le 2\sqrt{k+1}(k+1{)}^{\ln (k+1)/4}.$$

$$\Vert A-A(:,J)A(I,J{)}^{-1}A(I,:){\Vert }_{\max }\le 4^m\cdot {\rho }_m\cdot {\sigma }_{m+1}(A).$$

$$\Vert L_{11}^{-1}\Vert \le 2^m\cdot \min \{|p_1|,\dots ,|p_{m+1}|{\}}^{-1},$$

$$\Vert A_{11}^{-1}\Vert =\Vert U_{11}^{-1}{L}_{11}^{-1}\Vert \le 4^m\cdot \min \{|p_1|,\dots ,|p_{m+1}|{\}}^{-1}.$$

$$\min \{|p_1|,\dots ,|p_{m+1}|\}\le 4^m\Vert A_{11}^{-1}{\Vert }^{-1}=4^m{\sigma }_{m+1}(A_{11})\le 4^m{\sigma }_{m+1}(A),$$

$$p_{m+1}=\Vert A^{(m)}{\Vert }_{\max }\le {\rho }_j\cdot \Vert A^{(m-j)}{\Vert }_{\max }={\rho }_j\cdot |p_{m-j+1}|\le {\rho }_m\cdot |p_{m-j+1}|$$

$$\Vert A-A(:,J)A(I,J{)}^{-1}A(I,:){\Vert }_{\max }$$

$$=\Vert A^{(m)}{\Vert }_{\max }=|p_{m+1}|\le {\rho }_m\min \{|p_1|,\dots ,|p_{m+1}|\}.$$

$${\gamma }_k(A):=\min \{\Vert E{\Vert }_{\max }:\text{rank}(A+E)\le k\},k=1,\dots ,n.$$

$${\gamma }_{n-1}(A)=\Vert A^{-1}{\Vert }_{\infty \to 1}^{-1},$$

$$\Vert B{\Vert }_{\alpha \to \beta }:=\underset{x\ne 0}{\sup }\Vert Bx{\Vert }_{\beta }/\Vert x{\Vert }_{\alpha }$$

Full text

1 Introduction

Given a matrix $A\in\mathbb{R}^{n\times n}$ and $1\leq k\leq n$, the volume of a submatrix $A(I,J)$ for two index sets $I,J\subset\{1,\ldots,n\}$ of cardinality $k$ is defined as the absolute value of the determinant. The problem of finding the submatrix of maximum volume is connected to a range of applications in discrete mathematics, engineering, and scientific computing; see, e.g., [1, 17, 36]. Our primary motivation is its connection to low-rank approximation. Specifically, if $A(I,J)$ is invertible then the matrix $A(:,J)A(I,J)^{-1}A(I,:)$ has rank $k$ where the colon is used to denote that all rows or columns are selected. Low-rank approximations that involve columns and rows of the original matrix come in various flavors, as pseudoskeleton approximation [16], cross approximation [4], CUR approximation [11], or (strong) rank-revealing LU factorizations [27, 29]. If $A(I,J)$ has maximum volume then, by a result of Goreinov and Tyrtyshnikov [15], we have

[TABLE]

where $\|\cdot\|_{\max}$ denotes the maximum absolute value of the entries of a matrix. The $(k+1)$th largest singular value of $A$, denoted by $\sigma_{k+1}(A)$, is well known to govern the best rank-$k$ approximation error of $A$ in the spectral norm; see, e.g., [22, Chapter 7.4]. In other words, the result (1) states that volume maximization yields a quasi-best low-rank approximation.

Finding the submatrix of maximum volume is a difficult problem. In fact, it is NP hard to determine the maximum volume submatrix of $A$ [7, 30]. In [9] it is shown that there exists a universal constant $c>1$ such that it is NP hard to approximate the maximum volume of a $k\times k$ submatrix of a matrix $A\in\mathbb{R}^{n\times k}$ within a factor $c^{k}$. By a trivial embedding, this implies that it is also NP hard to approximate the maximum volume of a $k\times k$ submatrix of an $n\times n$ matrix. However, it is important to emphasize that while a good approximation of the maximum volume submatrix yields a good low-rank approximation [15, Theorem 2.2], the converse is generally not true, see also Remark 11 below.

Despite the difficulties associated with volume maximization, this concept has proven fruitful in the development of greedy and randomized algorithms that often yield reasonably good low-rank approximations. In particular, (adaptive) cross approximation [4], a greedy method for volume maximization, can now be regarded as the work horse for matrices $A$ that cannot be stored in memory, because, for example, $A$ has too many nonzero entries or it is too expensive to compute all the entries of $A$. This situation occurs frequently for discretized integral operators and cross approximation plays an important role in accelerating computations within the boundary element method [4] and uncertainty quantification [18]. Low-rank approximations of the form (1) also feature prominently in the Nyström method for kernel-based learning [2] and spectral clustering [13]. Let us also stress that cross approximation is equivalent to Gaussian elimination and (incomplete) LU factorization with complete pivoting; it primarily constitutes a different point of view with stronger emphasis on low-rank approximation.

In many of the applications mentioned above, the matrix $A$ carries additional structure. For example, if $A$ is the discretization of an integral operator with a positive semidefinite kernel then $A$ is symmetric positive semidefinite. In the first part of this work, we will show that the submatrix of maximum volume is always attained by a principal submatrix if $A$ is symmetric positive semidefinite (SPSD). This has a number of important consequences. For example, it allows us to draw a one-to-one correspondence to the column selection problem considered in [7]. In turn, the maximum volume problem remains NP hard when restricted to SPSD matrices. We also extend this result to diagonally dominant (DD) matrices. Somewhat surprisingly, we have not found such results for SPSD and DD matrices in the existing literature.

In the second part of this work, we derive a priori error bounds for the approximation returned by cross approximation. Although the literature on rank-revealing LU factorizations contains related results, see in particular [12, Corollary 5.3], the non-asymptotic bound of Theorem 6 appears to be new. Our result includes existing work by Harbrecht et al. [18] for SPSD matrices as a special case. For the particular case of doubly DD matrices, we show that the approximation error returned by cross approximation is at most $2(k+1)$ times larger than the right-hand side of (1). This class of matrices includes symmetric DD matrices, which play a prominent role in [24, 33]. Our result also allows us to obtain refined bounds for the convergence of cross approximation applied to functions [4, 35].

2 Maximum volume submatrices

In this section, we will prove for two classes of matrices that the submatrix of maximum volume can always be chosen to be a principal submatrix, that is, a submatrix of the form $A(I,I)$.

2.1 Symmetric positive semidefinite (SPSD) matrices

Theorem 1**.**

Let $A\in\mathbb{R}^{n\times n}$ be SPSD and let $1\leq k\leq n$. Then the maximum volume $k\times k$ submatrix of $A$ can be chosen to be a principal submatrix.

Proof.

Let $A(I,J)$ be any $k\times k$ submatrix of $A$. As $A$ is SPSD, it admits a Cholesky decomposition $A=C^{*}C$ and, in turn, $A(I,J)=C(:,I)^{*}C(:,J)$. The singular values of a principal submatrix satisfy

[TABLE]

Noting that the absolute value of the determinant equals the product of the singular values, we obtain

[TABLE]

where we used [21, Theorem 3.3.4] for the inequality. This implies that the volume of $A(I,J)$ is not larger than the maximum of the volumes of $A(I,I)$ and $A(J,J)$. In turn, $A(I,J)$ can be replaced by a principal submatrix without decreasing the volume. ∎

Trivially, the result of Theorem 1 extends to symmetric negative semidefinite matrices. On the other hand, it does not extend to the indefinite case; consider for example the $2k\times 2k$ matrix $A=\begin{bmatrix}0&I\\ I&0\end{bmatrix}$.

2.1.1 Connection to column selection

The volume of a general $n\times k$ matrix is defined as the product of its singular values. The column selection problem, which is also connected to low-rank approximation [8], is the following:

Given $B\in\mathbb{R}^{n\times m}$ and $1\leq k<m$, select the $n\times k$ submatrix of maximum volume.

In [7] it is shown that this problem is NP hard.

Theorem 1 allows us to relate the column selection problem to the classical maximum volume submatrix problem. Given $B\in\mathbb{R}^{n\times m}$, we consider the SPSD matrix $A=B^{*}B\in\mathbb{R}^{m\times m}$. As $A(I,I)=B(:,I)^{*}B(:,I)$ for any index set $I$, there is a one-to-one correspondence between the principal submatrices of $A$ and the subsets of $k$ columns of $B$. Moreover, as seen in the proof of Theorem 1, the volume of $A(I,I)$ is the square of the volume of $B(:,I)$. This shows that $B(:,I)$ has maximum volume if and only if $A(I,I)$ has maximum volume. In turn, this proves that the maximum volume submatrix problem remains NP hard when restricted to the subclass of SPSD matrices.

2.2 Diagonally dominant (DD) matrices

Definition 2**.**

A matrix $A\in\mathbb{R}^{n\times n}$ is called (row) diagonally dominant (DD) if

[TABLE]

If (2) holds with strict inequality for $i=1,\ldots,n$, we call $A$ strictly DD. A matrix $A$ is called doubly DD if both $A$ and $A^{*}$ are DD.

Lemma 3**.**

Let $T\in\mathbb{R}^{n\times n}$ be a strictly DD, upper triangular matrix. Then $|\det(T(I,J))|<|\det(T(I,I))|$ holds for every $I,J\subseteq\{1,\ldots,n\}$ with $|I|=|J|$ and $I\not=J$.

Proof.

Let $D$ be the diagonal matrix with $d_{ii}=t_{ii}$ and set $\tilde{T}=D^{-1}T$. Because of

[TABLE]

the statement of the lemma holds for $T$ if and only if it holds for $\tilde{T}$. In turn, this allows us to assume without loss of generality that $T$ has ones on the diagonal. In particular, $\det(T(I,I))=1$.

The statement of the lemma will be proven by induction on $k:=|I|=|J|$. The case $k=1$ follows immediately from the diagonal dominance of $T$. Suppose now that the statement of the lemma is true for fixed $k$. To prove the statement for $k+1$, we consider an arbitrary $(k+1)\times(k+1)$ submatrix $B:=T(I,J)$. If $I\not=J$ then there exists a row $B(i,:)$ that does not contain a diagonal element of $T$. By diagonal dominance of $T$,

[TABLE]

Denote by $B_{ij}$ the $k\times k$ submatrix of $T$ obtained from eliminating the $i$th row and $j$th column of $B$. By induction assumption, $|\det(B_{ij})|\leq 1$. Thus, combining (3) with the Laplace expansion gives

[TABLE]

In other words, $|\det(T(I,J))|<|\det(T(I,I))|$. ∎

Theorem 4**.**

Let $A\in\mathbb{R}^{n\times n}$ be a diagonally dominant matrix and $1\leq k\leq n$. Then the maximum volume $k\times k$ submatrix of $A$ can be chosen to be a principal submatrix.

Proof.

We prove the theorem in the case when $A$ is strictly DD; the DD case follows by a continuity argument, noting that volumes of submatrices are continuous in $A$. Let $A(I,J)$ be a $k\times k$ submatrix of $A$. Also, by applying a suitable permutation to the rows and columns of $A$, we may assume that $I=\{1,\ldots,k\}$ and $J=\{k-d+1,\ldots,2k-d\}$ with $d=|I\cap J|$. The result of the theorem follows if we can prove

[TABLE]

For this purpose, we note that the LU factorization $A=LU$ always exists with $U$ strictly DD; see Theorem 9.9 in [20]. We have that

[TABLE]

As $L(I,I)$ is lower triangular with ones on the diagonal, we obtain

[TABLE]

Thus, the inequality (4) follows from Lemma 3. ∎

For $k=n-1$, the result of Theorem 4 is covered in the proof of Theorem 2.5.12 in [21], while the result of Lemma 3 for $k=n-1$ follows from Proposition 2.1 in [31].

3 Cross Approximation

In the following, we summarize the idea behind Bebendorf’s cross approximation algorithm [4]. For this purpose, we first recall that an approximation of the form $A(:,J)A(I,J)^{-1}A(I,:)$ is closely connected to an incomplete LU decomposition of $A$. To see this, suppose that $A$ has been permuted such that $I=J=\{1,\ldots,k\}$ and partition

[TABLE]

Assume that $A_{11}$ is invertible and admits an LU decomposition $A_{11}=L_{11}U_{11}$, where $L_{11}$ is lower triangular and $U_{11}$ is upper triangular with ones on the diagonal. By setting $L_{21}=A_{21}U_{11}^{-1}$ and $U_{12}=L_{11}^{-1}A_{12}$, we obtain

[TABLE]

with the Schur complement

[TABLE]

This shows that the approximation error is governed by $A^{(k)}$. The factorized form (5) corresponds exactly to what is obtained after applying $k$ steps of the LU factorization to $A$, see, e.g., [14, Chapter 3.2].

Given index sets $I$ and $J$, one step of the greedy method for volume maximization consists of choosing indices such that

[TABLE]

Again, let us assume that $I=J=\{1,\ldots,k\}$ and set $\tilde{I}=I\cup\{k+\tilde{i}\}$, $\tilde{J}=J\cup\{k+\tilde{j}\}$ for some $\tilde{i},\tilde{j}\in\{1,\ldots,n-k\}$. Then (6) implies

[TABLE]

In other words, the local optimization problem (7) is solved by searching the entry of $A^{(k)}$ that has maximum modulus. This choice leads to Algorithm 1, which is equivalent to applying LU factorization with complete pivoting to $A$.

Remark 5**.**

Because of (5), the remainder term of Algorithm 1 satisfies $R_{k}=\begin{bmatrix}0&0\\ 0&A^{(k)}\end{bmatrix}$ after a suitable permutation of the indices. Both for SPSD and DD matrices, the element of maximum modulus is on the diagonal. Positive definiteness and diagonal dominance are preserved by taking Schur complements; see, e.g., [38, Chapter 4]. In turn, the search for the pivot element in Step 3 can be restricted to the diagonal for such matrices. This significantly reduces the number of entries of $A$ that need to be evaluated when running Algorithm 1. It also implies that Algorithm 1 returns $I=J$, which aligns nicely with the results from Section 2. Notice that if $A$ is an SPSD matrix then the cross approximation

[TABLE]

obtained by Algorithm 1 is SPSD. In contrast, diagonal dominance is generally not preserved by the low-rank approximation returned by Algorithm 1.

3.1 Error analysis for general matrices

Although not desirable, it may happen that the pivots $p_{k}$ in Algorithm 1 grow. Upper bounds on the growth factor $\|A^{(k)}\|_{\max}/\|A\|_{\max}$ play an important role in the error analysis of Gaussian elimination (see e.g. [37]). In the setting of complete pivoting, we can define

[TABLE]

where the supremum is taken over all matrices of rank at least $k$. This condition ensures that there is no breakdown in the first $k$ steps of Algorithm 1. By definition, $1\leq\rho_{1}\leq\rho_{2}\leq\ldots\leq\rho_{k}$. Wilkinson [37] proved that

[TABLE]

but it is known that this bound cannot be attained for $k\geq 3$. For matrices occurring in practice, it is rare to see any significant growth and it is not unreasonable to consider $\rho_{k}=O(1)$; we refer to [20, Section 9.4] for more details. Extending the proof of [18, Theorem 3.2], we obtain the following result.

Theorem 6**.**

Let $A\in\mathbb{R}^{n\times n}$ have rank at least $m<n$. Then the index sets $I,J$ returned by Algorithm 1 satisfy

[TABLE]

Proof.

Without loss of generality, we may assume $I=J=\{1,\ldots,m\}$. We perform one more step of Algorithm 1 and consider the relation $A_{11}=L_{11}U_{11}$ from (5) for $k=m+1$. Because of complete pivoting, the element of largest modulus in the $j$th column of $L_{11}$ is on the diagonal and equals $p_{j}$. For such triangular matrices, Theorem 6.1 in [19] gives

[TABLE]

where $\|\cdot\|$ denotes the spectral norm of a matrix. Analogously, using that the element of largest modulus in every row of $U_{11}$ is on the diagonal and equals $1$, we obtain $\|U_{11}^{-1}\|\leq 2^{m}$. Hence,

[TABLE]

This implies

[TABLE]

where we used interlacing properties of singular values, see [22, Corollary 7.3.6].

On the other hand, as $A^{(m)}$ is the matrix obtained after $j$ steps of Algorithm 1 applied to the matrix $A^{(m-j)}$, the definition (8) gives the inequalities

[TABLE]

for $j=1,2,\ldots,m$. We therefore obtain

[TABLE]

Combined with (10), this shows the result of the theorem. ∎

Because of the factor $4^{m}$, Theorem 6 only guarantees good low-rank approximations when the singular values are strongly decaying. This limitation does not correspond to the typical behavior observed in practice; the quantities $\|{L_{11}^{-1}}\|$ and $\|{U_{11}^{-1}}\|$ rarely assume the exponential growth estimates used in the proof of Theorem 6. In turn, the factor $4^{m}$ usually severely overestimates the error. Nevertheless, there are examples for which the error estimate of Theorem 6 is asymptotically tight; see Section 3.2 below.

The matrix norms on the two sides of the estimate (9) do not match. In the following we develop a variant of Theorem 6 in which the best approximation error is also measured in terms of $\|\cdot\|_{\max}$. This will be useful later on, in Section 3.5, when considering approximation of functions. Let us define the approximation numbers

[TABLE]

Because of $\|E\|/n\leq\|E\|_{\max}\leq\|E\|$, we have $\sigma_{k}(A)/n\leq\gamma_{k-1}(A)\leq\sigma_{k}(A)$. If $A$ is invertible then $\sigma_{n}(A)=\|A^{-1}\|^{-1}$. This result, relating the distance to singularity to the norm of the inverse, extends to general subordinate matrix norms; see, e.g., [20, Theorem 6.5]. In particular, we have

[TABLE]

with $\|\cdot\|_{\infty\to 1}$ denoting the matrix norm induced by the $1$- and $\infty$-norms. More generally, we set

[TABLE]

for vector norms $\|\cdot\|_{\alpha}$, $\|\cdot\|_{\beta}$. Note that $\|B\|_{1\to\infty}=\|B\|_{\max}$.

Theorem 7**.**

Under the assumptions of Theorem 6 and with the notation introduced above, we have

[TABLE]

Proof.

Along the lines of the proof of Theorem 6, we first note that

[TABLE]

which can be shown by induction on $m$. Combined with $\|U_{11}^{-1}\|_{1\to 1}\leq 2^{m}$, see [19, Theorem 6.1], we obtain

[TABLE]

where we used submultiplicativity [20, Eqn (6.7)]. Using (12), this implies

[TABLE]

with the last inequality being a direct consequence of the definition of $\gamma_{m}$. The rest of the proof is identical with the proof of Theorem 6. ∎

3.2 Error analysis for SPSD matrices

In the SPSD case, the pivot elements of Algorithm 1 are always non-increasing. Thus, when restricting the supremum in (8) to SPSD matrices of rank at least $k$, one obtains $\rho_{k}=1$. In turn, the following result due to Harbrecht et al. [18] is a corollary of Theorem 6.

Corollary 8**.**

For an SPSD matrix $A$ of rank at least $m$, the bound of Theorem 6 improves to

[TABLE]

The bound of Corollary 8 is asymptotically tight, see [18, Remark 3.3] and [23, p. 791]. As the growth factor $\rho_{m}$ which comes into play in Theorem 6 is small compared to the $4^{m}$ factor, this also proves that the bound of Theorem 6 is almost tight.

3.3 Error analysis for DD matrices

When restricting the supremum in (8) to DD matrices of rank at least $k$, one obtains $\rho_{k}\leq 2$; see Theorem 13.8 in [20].

Corollary 9**.**

For a DD matrix $A$ of rank at least $m$, the bound of Theorem 6 improves to

[TABLE]

Proof.

It is well known that the factor $U$ in the LU decomposition of a DD matrix is again DD; see [6]. In particular, this implies that the $(m+1)\times(m+1)$ unit upper triangular matrix $U_{11}$ in the proof of Theorem 6 is DD. Then, for every entry of $U_{11}^{-1}$ we have $|(U_{11}^{-1})_{ij}|\leq 1$ by [31, Prop. 2.1]. Therefore,

[TABLE]

This shows that the factor $4^{m}$ can be reduced to $(m+1)2^{m}$ in the bound of Theorem 6. Combined with $\rho_{m}\leq 2$, this establishes the desired result. ∎

Corollary 10**.**

For a doubly DD matrix $A$ of rank at least $m$, the bound of Corollary 9 improves to

[TABLE]

Proof.

Trivially, $A^{*}$ is DD. By the same arguments as in the proof of Corollary 9 this implies that not only $U_{11}$ but also $L_{11}^{*}$ is DD. Proceeding as in the derivation of (13), we get

[TABLE]

This shows that the factor $(m+1)\cdot 2^{m+1}$ of Corollary 9 can be improved to $2(m+1)^{2}$. ∎

Remark 11**.**

The fact that Algorithm 1 gives a polynomially good low-rank approximation of a doubly DD matrix does not imply that it also gives a polynomially good approximation of the maximum volume submatrix. For instance, let $n=2m$ and consider $A=\begin{bmatrix}I_{m}&0\\ 0&B_{m}\end{bmatrix}$, where $B_{m}=\,\mathrm{tridiag}[\frac{1}{2},1,-\frac{1}{2}]$. Then Algorithm 1 does not perform any pivoting during its $m$ steps and thus the submatrix $I_{m}$ is selected. Its volume is $1$, while the volume of $B_{m}$ is exponentially larger, it grows like $\left(\frac{1+\sqrt{2}}{2}\right)^{m}$.

3.4 Tightness of estimates for DD matrices

To study the tightness of the estimates from Section 3.3, it is useful to connect Algorithm 1 to LDU decompositions. Suppose that the application of $m=n-1$ steps of Algorithm 1 yields $I=J=\{1,\ldots,n-1\}$. As in the proof of Theorem 6, we exploit the relation (5) for $k=m+1=n$ to obtain the factorization

[TABLE]

where $L:=L_{11}D^{-1}$ and $U:=U_{11}$ are lower and upper unit triangular matrices, respectively. We recall from (11) that the error of the approximation returned by Algorithm 1 is governed by $|p_{n}|$.

From now on, let $A$ be a DD matrix. In this case, the pivot growth factor does not exceed $2$ and we have that $|p_{n}|\leq 2\|D^{-1}\|^{-1}\leq 2|p_{n}|$. In turn, the ratio between $|p_{n}|$ and the best rank-$(n-1)$ approximation error satisfies

[TABLE]

Inheriting the diagonal dominance from $A$, the matrix $U$ is well conditioned; see (13). Therefore, large $r_{m}$ require $\|L^{-1}\|$ to become large.

The quantity $\|L^{-1}\|$ also plays a prominent role in the stability analysis of LDU decompositions, see [10] and the references therein. In particular, the potential rapid growth of $\|L^{-1}\|$ under complete pivoting has motivated the search for alternative pivoting strategies [31]. However, the existing literature is scarce on examples actually exhibiting such rapid growth. The worst example we could find is by Barreras and Peña [3, Sec. 3], which exhibits linear growth. A more rapid growth is attained by the $n\times n$ matrix

[TABLE]

where $n$ is even and each block has size $n/2\times n/2$. When applying complete pivoting to this matrix, no interchanges are performed and the LDU factorization satisfies

[TABLE]

Note that, for this example, the right-hand side of (14) overestimates the error. This example attains quadratic growth: $r_{m}=\|A^{-1}\|=\Theta(m^{2})$. This is still far away from the exponential growth estimated in Corollary 9, but closer than the example from [3, Sec. 3], which yields $r_{m}=\Theta(m\sqrt{m})$.

For a doubly DD matrix, one obtains linear growth in (14) by considering the $n\times n$ lower bidiagonal matrix $B$ having $1$ on the diagonal and $-1$ on the first subdiagonal. In this case, $L=B$, $D=U=I_{n}$ and hence $\|B_{n}^{-1}\|=\Theta(m)$, showing that $r_{m}$ can grow at least linearly with $m$. We have not found an example exhibiting the quadratic growth estimated by Corollary 10.

3.5 Cross approximation for functions

Let us consider the approximation of a function $f:[-1,1]^{2}\to\mathbb{R}$ by a sum of separable functions:

[TABLE]

In the context of cross approximation, the factors are restricted to functions $f_{i}^{(1)}$ of the form $f_{i}^{(1)}=f(x,\bar{y}_{i})$ and $f_{i}^{(2)}=f(\bar{x}_{i},y)$, where $\bar{x}_{i}$ and $\bar{y}_{i}$ are fixed elements of $[-1,1]$. In particular, Micchelli and Pinkus [26] considered interpolating approximations of the following form:

[TABLE]

for some $x_{1},\ldots,x_{m},y_{1},\ldots,y_{m}\in[-1,1]$. Townsend and Trefethen [35] use a strategy for choosing the interpolation points which is basically equivalent to Algorithm 1 and they prove a convergence result under some analyticity hypotheses on the function $f$. There also exist error analyses for cross approximation of functions when using different pivoting strategies, see, e.g., [5, 32].

Algorithm 2 summarizes cross approximation of functions with complete pivoting.

We now explain the connection to Algorithm 1. Fix $(x,y)\in[-1,1]^{2}$ and consider the points $x_{1},\ldots,x_{m}$ and $y_{1},\ldots,y_{m}$ obtained by the first $m$ steps of Algorithm 2. Consider what happens when applying Algorithm 1 to the $(m+1)\times(m+1)$ matrix obtained by interpolating $f$ in the points mentioned above:

[TABLE]

The first chosen pivot will be $p_{1}=f(x_{1},y_{1})$ because it is the largest entry of the matrix. Now observe that the Schur complement $A^{(1)}$ obtained after the first step, is the matrix that interpolates the function $e_{1}$ in the points $x_{2},\ldots,x_{m},x$ and $y_{2},\ldots,y_{m},y$. At this point, the second pivot chosen by Algorithm 1 will be $e_{1}(x_{2},y_{2})$ because of how Algorithm 2 chose $(x_{2},y_{2})$ in line $5$. After $m$ steps of Algorithm 1 we will be left with only one nonzero entry in position $(m+1,m+1)$ and this will be $e_{m}(x,y)$. This allows us to estimate $|e_{m}(x,y)|$ via Theorem 7:

[TABLE]

The last thing we need is an estimate on $\gamma_{m}(A_{(x,y)})$ that is uniform in $(x,y)\in[-1,1]^{2}$. This will follow from analyticity assumptions on the functions $f(\cdot,y)$ for $y\in[-1,1]$.

Definition 12**.**

The Bernstein ellipse $\mathcal{E}_{r}$ of radius $r>1$ is the ellipse with foci in $-1$ and $1$ and with sum of the semi-axes equal to $r$.

Corollary 13**.**

Let $f:[-1,1]^{2}\to\mathbb{R}$ be such that $f(\cdot,y)$ admits an analytic extension - which we will denote by $\tilde{f}$ - in the Bernstein ellipse $\mathcal{E}_{r_{0}}$ of radius $r_{0}$ for each $y\in[-1,1]$. Let $1<r<r_{0}$ and

[TABLE]

After $m$ steps of Algorithm 2 the error function satisfies

[TABLE]

Proof.

Fix $(x,y)\in[-1,1]^{2}$ and let $b:[-1,1]\to\mathbb{R}^{m+1}$ be the vector-valued function defined by

[TABLE]

The analyticity hypothesis allows us to apply standard polynomial approximation results (see e.g. Corollary 2.2 in [25]) and conclude that there exists an approximation $\hat{b}:[-1,1]\to\mathbb{R}^{m+1}$ given by

[TABLE]

where $v_{k}\in\mathbb{R}^{m+1}$ are constant vectors and $p_{k}:[-1,1]\to\mathbb{R}$ are polynomials, such that

[TABLE]

for any $1<r<r_{0}$. We can clearly bound $\max_{\alpha\in\mathcal{E}_{r}}\|b(\alpha)\|_{\max}\leq M$.

The matrix $A_{(x,y)}$ is obtained by sampling $b$ in the points $x_{1},\ldots,x_{m},x$, i.e.

[TABLE]

Let us define, analogously,

[TABLE]

Notice that $\hat{A}_{(x,y)}$ has rank as most $m$ because by (16) each of the $m+1$ columns of $\hat{A}_{(x,y)}$ is a linear combination of the $m$ vectors $v_{1},\ldots,v_{m}$, so

[TABLE]

The result then follows from Equation (15). ∎

To get convergence of the error function to zero as $m\to\infty$, in Corollary 13 it is sufficient that the function $f(\cdot,y)$ admits an analytic extension to the Bernstein ellipse $\mathcal{E}_{r_{0}}$ with $r_{0}>4$ for each $y$, because the factor $\rho_{m}$ has subexponential growth. Our result compares favorably to Theorem 8.1 in [35], which requires an analytic extension to the region $K$ consisting of all points at a distance $\leq 4$ from $[-1,1]$.

Figure 1 compares the two domains and it is evident that the requirement from [35] is significantly more restrictive.

For a positive semidefinite kernel function $f$, the matrix $A_{(x,y)}$ in (15) is positive semidefinite and hence the bound of Corollary 13 holds with $\rho_{m}\equiv 1$. This matches an asymptotic result given in [35, Theorem 9.1].

4 Conclusions

The fact that the search for the maximum volume submatrix can be restricted to principal submatrices for SPSD and DD matrices appears intuitive and is sometimes used without proof, see, e.g., [28, Theorem 1]. As far as we know, Theorems 1 and 4 are the first results providing a mathematical justification to this intuition.

For cross approximation, Theorem 6 appears to be the first non-asymptotic error bound that holds for general matrices. Except for [18], previous results for cross approximation applied to matrices or functions [5, 35] are based on a step-by-step analysis of the error. In contrast, our technique takes a more global view and can, in turn, leverage existing results on the pivot growth in Gaussian elimination. As illustrated in Section 3.5, this can yield significant advantages.

A number of fundamental questions remain open. Most importantly, there is a mismatch between the derived error bounds and the known worst-case examples for cross approximation applied to DD and doubly DD matrices. Especially for DD matrices, this problem appears to be difficult to overcome and was encountered previously in the context of the error analysis of LDU factorizations [3, 10].

5 Acknowledgments

We thank Froilan Dopico for an insightful discussion on DD matrices.