Given a matrix and , the volume of a submatrix for two index sets of cardinality 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 is invertible then the matrix has rank 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 , cross approximation , CUR approximation , or (strong) rank-revealing LU factorizations [27, 29]. If has maximum volume then, by a result of Goreinov and Tyrtyshnikov , we have
where denotes the maximum absolute value of the entries of a matrix. The
th largest singular value of, denoted by , is well known to govern the best rank- approximation error of 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 [7, 30]. In  it is shown that there exists a universal constant such that it is NP hard to approximate the maximum volume of a submatrix of a matrix within a factor . By a trivial embedding, this implies that it is also NP hard to approximate the maximum volume of a submatrix of an 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 , a greedy method for volume maximization, can now be regarded as the work horse for matrices that cannot be stored in memory, because, for example, has too many nonzero entries or it is too expensive to compute all the entries of . This situation occurs frequently for discretized integral operators and cross approximation plays an important role in accelerating computations within the boundary element method  and uncertainty quantification . Low-rank approximations of the form (1) also feature prominently in the Nyström method for kernel-based learning 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 carries additional structure. For example, if is the discretization of an integral operator with a positive semidefinite kernel then 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 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 . 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.  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 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 .
2.1 Symmetric positive semidefinite (SPSD) matrices
Let be SPSD and let . Then the maximum volume submatrix of can be chosen to be a principal submatrix.
Let be any submatrix of . As is SPSD, it admits a Cholesky decomposition and, in turn, . The singular values of a principal submatrix satisfy
Noting that the absolute value of the determinant equals the product of the singular values, we obtain
where we used [21, Theorem 3.3.4] for the inequality. This implies that the volume of is not larger than the maximum of the volumes of and . In turn, 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 matrix .
2.1.1 Connection to column selection
The volume of a general matrix is defined as the product of its singular values. The column selection problem, which is also connected to low-rank approximation , is the following:
Given and , select the submatrix of maximum volume.
In  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 , we consider the SPSD matrix . As for any index set , there is a one-to-one correspondence between the principal submatrices of and the subsets of columns of . Moreover, as seen in the proof of Theorem 1, the volume of is the square of the volume of . This shows that has maximum volume if and only if 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
A matrix is called (row) diagonally dominant (DD) if
If (2) holds with strict inequality for , we call strictly DD. A matrix is called doubly DD if both and are DD.
Let be a strictly DD, upper triangular matrix. Then holds for every with and .
Let be the diagonal matrix with and set . Because of
the statement of the lemma holds for if and only if it holds for . In turn, this allows us to assume without loss of generality that has ones on the diagonal. In particular, .
The statement of the lemma will be proven by induction on . The case follows immediately from the diagonal dominance of . Suppose now that the statement of the lemma is true for fixed . To prove the statement for , we consider an arbitrary submatrix . If then there exists a row that does not contain a diagonal element of . By diagonal dominance of ,
Denote by the submatrix of obtained from eliminating the th row and th column of . By induction assumption, . Thus, combining (3) with the Laplace expansion gives
In other words, . ∎
Let be a diagonally dominant matrix and . Then the maximum volume submatrix of can be chosen to be a principal submatrix.
We prove the theorem in the case when is strictly DD; the DD case follows by a continuity argument, noting that volumes of submatrices are continuous in . Let be a submatrix of . Also, by applying a suitable permutation to the rows and columns of , we may assume that and with . The result of the theorem follows if we can prove
For this purpose, we note that the LU factorization always exists with strictly DD; see Theorem 9.9 in . We have that
As is lower triangular with ones on the diagonal, we obtain
3 Cross Approximation
In the following, we summarize the idea behind Bebendorf’s cross approximation algorithm . For this purpose, we first recall that an approximation of the form is closely connected to an incomplete LU decomposition of . To see this, suppose that has been permuted such that and partition
Assume that is invertible and admits an LU decomposition , where is lower triangular and is upper triangular with ones on the diagonal. By setting and , we obtain
with the Schur complement
This shows that the approximation error is governed by . The factorized form (5) corresponds exactly to what is obtained after applying steps of the LU factorization to , see, e.g., [14, Chapter 3.2].
Given index sets and , one step of the greedy method for volume maximization consists of choosing indices such that
Again, let us assume that and set , for some . Then (6) implies
In other words, the local optimization problem (7) is solved by searching the entry of that has maximum modulus. This choice leads to Algorithm 1, which is equivalent to applying LU factorization with complete pivoting to .
Because of (5), the remainder term of Algorithm 1 satisfies 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 that need to be evaluated when running Algorithm 1. It also implies that Algorithm 1 returns , which aligns nicely with the results from Section 2. Notice that if is an SPSD matrix then the cross approximation
3.1 Error analysis for general matrices
Although not desirable, it may happen that the pivots in Algorithm 1 grow. Upper bounds on the growth factor play an important role in the error analysis of Gaussian elimination (see e.g. ). In the setting of complete pivoting, we can define
but it is known that this bound cannot be attained for . For matrices occurring in practice, it is rare to see any significant growth and it is not unreasonable to consider ; we refer to [20, Section 9.4] for more details. Extending the proof of [18, Theorem 3.2], we obtain the following result.
Let have rank at least . Then the index sets returned by Algorithm 1 satisfy
Without loss of generality, we may assume . We perform one more step of Algorithm 1 and consider the relation from (5) for . Because of complete pivoting, the element of largest modulus in the th column of is on the diagonal and equals . For such triangular matrices, Theorem 6.1 in  gives
where denotes the spectral norm of a matrix. Analogously, using that the element of largest modulus in every row of is on the diagonal and equals , we obtain . Hence,
where we used interlacing properties of singular values, see [22, Corollary 7.3.6].
Because of the factor , 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 and rarely assume the exponential growth estimates used in the proof of Theorem 6. In turn, the factor 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 . This will be useful later on, in Section 3.5, when considering approximation of functions. Let us define the approximation numbers
Because of , we have . If is invertible then . 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
with denoting the matrix norm induced by the - and -norms. More generally, we set
for vector norms, . Note that .
Under the assumptions of Theorem 6 and with the notation introduced above, we have
Along the lines of the proof of Theorem 6, we first note that
which can be shown by induction on . Combined with , see [19, Theorem 6.1], we obtain
with the last inequality being a direct consequence of the definition of . 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 , one obtains . In turn, the following result due to Harbrecht et al.  is a corollary of Theorem 6.
For an SPSD matrix of rank at least , the bound of Theorem 6 improves to
3.3 Error analysis for DD matrices
For a DD matrix of rank at least , the bound of Theorem 6 improves to
It is well known that the factor in the LU decomposition of a DD matrix is again DD; see . In particular, this implies that the unit upper triangular matrix in the proof of Theorem 6 is DD. Then, for every entry of we have by [31, Prop. 2.1]. Therefore,
This shows that the factor can be reduced to in the bound of Theorem 6. Combined with , this establishes the desired result. ∎
For a doubly DD matrix of rank at least , the bound of Corollary 9 improves to
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 and consider , where . Then Algorithm 1 does not perform any pivoting during its steps and thus the submatrix is selected. Its volume is , while the volume of is exponentially larger, it grows like .
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 steps of Algorithm 1 yields . As in the proof of Theorem 6, we exploit the relation (5) for to obtain the factorization
From now on, let be a DD matrix. In this case, the pivot growth factor does not exceed and we have that . In turn, the ratio between and the best rank- approximation error satisfies
Inheriting the diagonal dominance from , the matrix is well conditioned; see (13). Therefore, large require to become large.
The quantity also plays a prominent role in the stability analysis of LDU decompositions, see  and the references therein. In particular, the potential rapid growth of under complete pivoting has motivated the search for alternative pivoting strategies . 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 matrix
where is even and each block has size . When applying complete pivoting to this matrix, no interchanges are performed and the LDU factorization satisfies
Note that, for this example, the right-hand side of (14) overestimates the error. This example attains quadratic growth: . This is still far away from the exponential growth estimated in Corollary 9, but closer than the example from [3, Sec. 3], which yields .
For a doubly DD matrix, one obtains linear growth in (14) by considering the lower bidiagonal matrix having on the diagonal and on the first subdiagonal. In this case, , and hence , showing that can grow at least linearly with . 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 by a sum of separable functions:
In the context of cross approximation, the factors are restricted to functions of the form and , where and are fixed elements of . In particular, Micchelli and Pinkus 
considered interpolating approximations of the following form:
for some . Townsend and Trefethen  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 . 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 and consider the points and obtained by the first steps of Algorithm 2. Consider what happens when applying Algorithm 1 to the matrix obtained by interpolating in the points mentioned above:
The first chosen pivot will be because it is the largest entry of the matrix. Now observe that the Schur complement obtained after the first step, is the matrix that interpolates the function in the points and . At this point, the second pivot chosen by Algorithm 1 will be because of how Algorithm 2 chose in line . After steps of Algorithm 1 we will be left with only one nonzero entry in position and this will be . This allows us to estimate via Theorem 7:
The last thing we need is an estimate on that is uniform in . This will follow from analyticity assumptions on the functions for .
The Bernstein ellipse of radius is the ellipse with foci in and and with sum of the semi-axes equal to .
Let be such that admits an analytic extension - which we will denote by - in the Bernstein ellipse of radius for each . Let and
After steps of Algorithm 2 the error function satisfies
Fix and let be the vector-valued function defined by
The analyticity hypothesis allows us to apply standard polynomial approximation results (see e.g. Corollary 2.2 in ) and conclude that there exists an approximation given by
where are constant vectors and are polynomials, such that
for any . We can clearly bound .
To get convergence of the error function to zero as , in Corollary 13 it is sufficient that the function admits an analytic extension to the Bernstein ellipse with for each , because the factor has subexponential growth. Our result compares favorably to Theorem 8.1 in , which requires an analytic extension to the region consisting of all points at a distance from .