Abstract
The problem is to select \(k\) columns with maximal volume from a real \(m\times n\) matrix \(X\). We show that the logarithm of the volume is a submodular function on columns of \(\mx\), and for full column-rank matrices \(\mx\) with sufficiently large singular values, it is a non-negative non-decreasing function. As a consequence, traditional Businger-Golub QR with column pivoting is a greedy algorithm, with a relative error of at most 37 percent. In contrast, Gu-Eisenstat strong rank-revealing QR is a 1-interchange algorithm, with a relative error of at most 50 percent. The higher accuracy, under this metric, of the simple QR with column pivoting confirms its well known effectiveness in practice. For general matrices, we show that Businger-Golub QR is a greedy algorithm, and Gu-Eisenstat QR a 1-interchange algorithm for maximizing the trace of the upper triangular matrix in a QR factorization of \(\mx\), with Businger-Golub again having a better error bound. This is extended to finding \(k\times k\) submatrices of maximal volume in symmetric positive-definite matrices.
Submodular functions, greedy algorithm, interchange algorithm, volume, trace, pivoted QR factorization, pivoted Cholesky factorization
15A12, 15A15, 15A18, 15A23, 65F20, 90B80, 90C59
1 Introduction↩︎
Many problems in combinatorial optimization [1], machine learning [2], [3], management science [4], determinantal point processes [5], [6], matroids [7], sensor placement [8], [9], and D-optimal design [10] can be expressed as the maximization of a submodular function under cardinality constraints. Such maximization problems can be NP-complete, such as the account location problem [4], and heuristics have been developed, including greedy and interchange algorithms that come with relative error bounds [1].
To take advantage of these error bounds in the context of selecting \(k\) ‘good’ columns from a matrix \(\mx\in\rmn\), we express column subset selection in terms of volume, and show that the log-volume (logarithm of the volume) is a submodular function on column submatrices of \(\mx\).
For matrices with full column-rank and sufficiently large singular values, we show that traditional Businger-Golub QR with column pivoting [11], GoVL13? is a greedy algorithm for log-volume maximization, with a relative error of at most 37 percent. In contrast, Gu-Eisenstat strong rank-revealing QR [12] is a 1-interchange algorithm for log-volume maximization, with a relative error of at most 50 percent –however, only in theory, when run without a relaxation factor. The higher accuracy, under this metric, of the simple QR with column pivoting confirms its well known empirical effectiveness, and competitiveness with more sophisticated pivoting strategies.
For general, possibly rank-deficient matrices, however the log-volume can be negative and decreasing, hence the error bounds in [1] do not apply. As an alternative, we consider the trace of the triangular matrix in a pivoted QR factorization3. We show that Businger-QR with column pivoting approximates it with a relative error of at most 37 percent, while Gu-Eisenstat strong-rank revealing QR has a relative error of at most 50 percent. So, this alternative metric as well suggests that QR with column pivoting can be more accurate.
1.1 Contributions↩︎
We express existing algorithms for column subset selection and submatrix selection as greedy and 1-interchange algorithms for submodular functions, and derive error bounds on their accuracy.
We introduce variants of 1-interchange algorithms for submodular functions with additive and multiplicative relaxation factors, and derive error bounds for their accuracy (Theorem 3).
We show that the log-volume of a general matrix is a submodular function (Theorem 4), and introduce a new class of submodular functions acting on columns of a general matrix (Definition 5 and Theorem 5).
We show that Businger-Golub QR with column pivoting [11], GoVL13? is a greedy algorithm for log-volume maximization for matrices that have full column-rank (Theorem 6 and Corollary 2), with a relative error of at most 37 percent.
We show that Gu-Eisenstat strong rank-revealing QR [12] is a 1-interchange algorithm for log-volume maximization for matrices that have full column rank (Theorem 8, Corollary 3 and Theorem 9), with a relative error of at most 50 percent.
According to this metric, the traditional Businger-Golub QR with column pivoting can be more accurate than the strong rank-revealing Gu-Eisenstat QR factorization. This explains the accuracy of QR with column pivoting frequently observed in practice.
For general matrices we show that Businger-Golub QR with column pivoting [11], GoVL13? is a greedy algorithm (Theorem 7), and that Gu-Eisenstat strong rank-revealing QR [12] is a 1-interchange algorithm (Theorem 11) for maximizing the trace of the upper triangular matrix in a QR factorization; and again Businger-Golub can be more accurate.
For symmetric positive definite matrices, we show that pivoted Cholesky factorizations [13], [14], GoVL13? and strong rank-revealing Cholesky factorization [13] are, respectively, greedy and 1-interchange algorithms for log-volume maximization of principal submatrices (Theorems 13 and 16).
1.2 Overview↩︎
After reviewing submodular functions and their algorithms (Section 2), we discuss submodular functions for column subset selection (Section 3). Then we express and analyze the Businger-Golub QR factorization as a greedy algorithm (Section 4); the Gu-Eisenstat QR factorization as a 1-interchange algorithm (Section 5); and pivoted Cholesky factorizations as greedy and 1-interchange algorithms for submatrix selection (Section 6); followed by possible directions of future work (Section 7). Each section contains a brief review of relevant literature.
1.3 Notation↩︎
Let \(\mx\in\rmn\) be matrix with singular values \(\sigma_1(\mx)\geq \cdots\geq\sigma_{\min\{m,n\}}(\mx)\geq 0\). For some \(0<k\leq \rank(\mx)\), let \(\mc\in\real^{m\times k}\) be a submatrix of \(k\) columns of \(\mx\), with singular values \[\label{e95so} \sigma_1(\mc)\geq \cdots\geq\sigma_k(\mc)\geq 0.\tag{1}\] Singular value interlacing [15] implies \[\label{e95inter} \sigma_{\min\{m,n\}-k+j}(\ma)\leq \sigma_j(\mc)\leq \sigma_j(\ma),\qquad 1\leq j\leq k.\tag{2}\] The two-norm condition number with respect to left inversion of \(\mx\in\rmn\) with \(\rank(\mx)=n\) is \(\kappa_2(\mx) \equiv \sigma_1(\mx)/\sigma_n(\mx)\).
The eigenvalues of a symmetric positive semi-definite matrix \(\ma\in\rnn\) are \(\lambda_1(\ma)\geq \cdots\geq \lambda_n(\ma)\geq 0\). For symmetric positive semi-definite matrices \(\ma, \mb\in\rmm\), the Loewner inequality \(\ma\preceq \mb\) says that \(\mb-\ma\) is positive semi-definite.
We denote by \(\mi_n\equiv\begin{bmatrix} \ve_1& \cdots &\ve_n\end{bmatrix}\in\rnn\) the identity matrix, and by \(\ve_j\in\rn\) its columns for \(1 \le j \le n\). The norm \(\|\cdot\|\) is the Euclidean two-norm.
Definition 1. For a matrix \(\mz\in\real^{m\times k}\) we denote by \(\mathcal{R}(\mz)\in\real^{k\times k}\) the upper triangular matrix with non-negative diagonal elements in the thin QR decomposition \(\mz=\mathcal{Q}(\mz)\,\mathcal{R}(\mz)\), where \(\mathcal{Q}(\mz)\in\real^{m\times k}\) has orthonormal columns.
Definition 2. Let \(\mx\in\rmn\) and \(\mathcal{S}=\{s_1,\ldots,s_k\} \subset\{1,\ldots, n\}\) for some \(1\leq k\leq n\). Then \(\mx_{\mathcal{S}}=\begin{bmatrix}\vx_{s_1} &\ldots & \vx_{s_k}\end{bmatrix}\in\real^{m\times k}\) is a submatrix of \(\mx\) with columns indexed by \(\mathcal{S}\).
If \(\ma\in\rnn\), then \(\ma[\mathcal{S}]\in\real^{k\times k}\) is the principal submatrix of \(\ma\) indexed by rows and columns in \(\mathcal{S}\).
2 Submodular functions↩︎
After presenting the problem of maximization under cardinality constraints, and the definition of submodular functions (Definition 3), we review greedy algorithms and interchange algorithms (Section 2.1). We introduce and analyze a relaxed version of the interchange algorithms (Section 2.2), and review existing work (Section 2.3).
2.0.0.1 Maximization under cardinality constraints
Given a finite set \(\mathcal{X}\), a real valued function \(f\) on the subsets of \(\mathcal{X}\), and an integer \(0<k<|\mathcal{X}|\), we want to maximize \(f\) over all subsets \(\mathcal{S}\subset\mathcal{X}\) of cardinality \(k\), that is, \[\max_{\mathcal{S}\subset\mathcal{X}}\{f(\mathcal{S}):\;|\mathcal{S}|=k\}.\]
The notion of submodularity, defined below, is ‘in some sense a combinatorial analogue of concavity’ [1].
Definition 3 ((1.4), (1.5) and Proposition 2.1 in [1]). Given a finite set \(\mathcal{X}\), a real-valued function \(f\) on the subsets of \(\mathcal{X}\) is a submodular function* if \[f(\mathcal{T}\cup \{x\}) - f(\mathcal{T})\leq f(\mathcal{S}\cup\{x\})-f(\mathcal{S})\quad \text{for all} \quad \mathcal{S}\subset\mathcal{T}\subset\mathcal{X},\quad x\in\mathcal{X}-\mathcal{T}.\] Furthermore, a submodular function \(f\) is non-decreasing if \[f(\mathcal{T}\cup \{x\}) -f(\mathcal{T})\geq 0\quad \text{for all}\quad \mathcal{T}\subset\mathcal{X}, \quad x\in\mathcal{X}-\mathcal{T}.\]*
That is, adding a new element \(x\) to a smaller set \(\mathcal{S}\) is more beneficial than adding \(x\) to a larger set \(\mathcal{T}\); and adding a new element \(x\) can possibly increase the function, but cannot decrease it. For instance, linear functions are submodular.
2.1 Algorithms↩︎
We review two maximization algorithms: a greedy algorithm (Algorithm 1) and its accuracy (Theorem 1); and a 1-interchange algorithm (Algorithm 2) and its accuracy (Theorem 2).
The greedy algorithm in [1] is more general than Algorithm 1 because it applies to general submodular functions \(f\) which are not necessarily non-decreasing. However, here we need a non-negative non-decreasing function \(f\) so that we can set \(f(\emptyset)=0\).
![Figure 1: Greedy heuristic for non-negative non-decreasing submodular functions f[2], [1]](/article/2607.13823/_pdflatex/fwremskn.png "fig:")
Theorem 1 (Theorem 1 in [4], Theorem 5 in [2], Proposition 4.3 in [1]). Let \(f\) be a non-negative non-decreasing submodular function on a finite set \(\mathcal{X}\) with \(f(\emptyset)=0\); \(0<k<|\mathcal{X}|\) an integer; and \(\mathcal{S}_k^*\) a set that maximizes \(f(\mathcal{S})\) over all subsets \(\mathcal{S}\subset\mathcal{X}\) with \(|\mathcal{S}|=k\).
Then the set \(\mathcal{S}_k\) returned by Algorithm 1 satisfies \[0\leq f(\mathcal{S}_k^*)-f(\mathcal{S}_k) \leq (1-\tfrac{1}{k})^k\,f(\mathcal{S}_k^*),\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e} <.37\).
The upper bound in Theorem 1 is tight for the uncapacitated location problem [4], [1]. For certain classes of problems, [16] shows that the upper bound cannot be improved by an algorithm with a polynomial number of function evaluations.
![Figure 2: 1-interchange algorithm [1]](/article/2607.13823/_pdflatex/nxwtagrq.png "fig:")
Theorem 2 (Theorem 5 in [4], Theorem 5.1 in [1]). Let \(f\) be a non-negative non-decreasing submodular function on a finite set \(\mathcal{X}\) with \(f(\emptyset)=0\); \(0<k<|\mathcal{X}|\) an integer; and \(\mathcal{S}_k^*\) a set that maximizes \(f(\mathcal{S})\) over all subsets \(\mathcal{S}\subset\mathcal{X}\) with \(|\mathcal{S}|=k\).
Then the set \(\mathcal{S}_k\) returned by Algorithm 2 satisfies \[0\leq f(\mathcal{S}_k^*)-f(\mathcal{S}_k) \leq \frac{k-1}{2k-1}\,f(\mathcal{S}_k^*),\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\).
The upper bound is tight for the uncapacitated location problem [1].
2.2 Algorithm 2 with relaxation factors↩︎
We adapt the analysis of [1] to the case where interchanges take place only if the objective function improves by an additive relaxation factor \(\theta_a \ge 0\), or by a multiplicative relaxation factor \(\theta_m \ge 1\). For an additive factor, the condition in the while loop in line 2 of Algorithm 2 becomes \[\label{e95add} f( (\mathcal{S}^{(i)} - \{x_i\}) \cup \{y_j\} ) > f(\mathcal{S}^{(i)}) + \theta_a.\tag{3}\] The output \(\mathcal{S}_k\) of Algorithm 2 is locally optimal because for every \(x \in \mathcal{T}_k\) and \(y \in \{1,\ldots,n\} - \mathcal{T}_k\), \[f( (\mathcal{T}_k - \{ x \}) \cup \{y\}) \le f(\mathcal{S}_k) + \theta_a.\] For a multiplicative factor, the condition in the while loop in line 2 of Algorithm 2 becomes \[\label{e95mult} f( (\mathcal{S}^{(i)} - \{x_i\}) \cup \{y_j\} ) > \theta_m f(\mathcal{S}^{(i)}).\tag{4}\] Again, the output \(\mathcal{S}_k\) from Algorithm 2 is locally optimal.
Theorem 3. Let \(f\) be a non-negative non-decreasing submodular function on a finite set \(\mathcal{X}\) with \(f(\emptyset)=0\); \(0<k<|\mathcal{X}|\) an integer; and \(\mathcal{S}_k^*\) a set that maximizes \(f(\mathcal{S})\) over all subsets \(\mathcal{S}\subset\mathcal{X}\) with \(|\mathcal{S}|=k\). If \(\mathcal{S}_k\) is the output of the 1-interchange Algorithm 2 with additive relaxation factor \(\theta_a \ge 0\), then \[\begin{align} 0\leq f(\mathcal{S}_k^*)-f(\mathcal{S}_k) \leq\frac{k-1}{2k-1}f(\mathcal{S}_k^*)+\frac{k^2}{2k-1}\theta_a, \end{align}\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\) and \(\frac{k^2}{2k-1}\leq k\).
If \(\mathcal{S}_k\) is the output of the 1-interchange Algorithm 2 with multiplicative relaxation factor \(\theta_m\geq 1\), then \[\begin{align} 0\leq f(\mathcal{S}_k^*)-f(\mathcal{S}_k) \leq \frac{k-1 +k^2(\theta_m-1)}{2k-1 +k^2(\theta_m-1)} f(\mathcal{S}_k^*). \end{align}\]
In the special case \(\theta_a = 0\) or \(\theta_m = 1\), Theorem 3 reduces to Theorem 2.
Proof. We proceed as in the proof of [1]. Denote the elements of the optimal solution by \(\mathcal{S}_k^* = \{s_1^*,\dots,s_k^*\}\).
Apply Algorithm 1 to the output \(\mathcal{S}_k\) of Algorithm 2 to order the elements \(\mathcal{S}_k = \{s_1,\dots,s_k\}\), and define the subsets \(\mathcal{S}_{j} \equiv \{s_1,\dots,s_j\}\), \(1\le j \le k\) and \(\mathcal{S}_0 \equiv \emptyset\).
Define the increments \(\zeta_j \equiv f(\mathcal{S}_s) - f(\mathcal{S}_{s-1})\), \(1 \le j \le k\). Then \[\begin{align} \zeta_{j} &= f(\mathcal{S}_j) - f(\mathcal{S}_{j-1}) = f(\mathcal{S}_{j-1} \cup \{s_j\}) - f(\mathcal{S}_{j-1}) \le \> f(\mathcal{S}_{j-2} \cup \{s_j\}) - f(\mathcal{S}_{j-2}) \\ &\le \> f(\mathcal{S}_{j-2} \cup \{s_{j-1}\}) - f(\mathcal{S}_{j-2}) = f(\mathcal{S}_{j-1}) - f(\mathcal{S}_{j-2}) = \zeta_{j-1}, \qquad 2\leq j\leq k, \end{align}\] where the first inequality follows from submodularity, and the second from greedy ordering. Thus, the increments are non-increasing \(\zeta_1 \ge \dots \ge \zeta_{k}\). With \(f(\emptyset)=0\) this implies \[\label{eqn:decomp} f(\mathcal{S}_k) = f(\emptyset) + \sum_{j=1}^{k}\zeta_j \ge k \zeta_{k}.\tag{5}\] From [1] or [1] follows \[\label{eqn:decomp2} f(\mathcal{S}_k^*) \le f(\mathcal{S}_{k-1}) + \sum_{s_j^* \in \mathcal{S}_k^* - \mathcal{S}_{k-1}} (f(\mathcal{S}_{k-1} \cup \{s_j^*\} ) - f(\mathcal{S}_{k-1})).\tag{6}\] We consider additive and multiplicative relaxation factors separately.
2.2.0.1 Additive relaxation
The interchange property (3 ) ensures \[f(\mathcal{S}_{k-1} \cup \{s_j^*\} ) \le f(\mathcal{S}_k) + \theta_a , \qquad s_j^* \in \mathcal{S}_k^* - \mathcal{S}_{k-1}.\] Subtract \(f(\mathcal{S}_{k-1})\) on both sides, \[f(\mathcal{S}_{k-1} \cup \{s_j^*\} ) - f(\mathcal{S}_{k-1}) \le f(\mathcal{S}_k) - f(\mathcal{S}_{k-1}) +\theta_a = \zeta_{k} +\theta_a , \qquad s_j^* \in \mathcal{S}_k^* - \mathcal{S}_{k-1}.\] Insert this inequality into 6 \[\begin{align} f(\mathcal{S}_k^*) \le & \> f(\mathcal{S}_{k-1}) + | \mathcal{S}_k^*-\mathcal{S}_{k-1} |(\zeta_{k} + \theta_a) \\ \le & \> f(\mathcal{S}_{k-1}) + k\zeta_k + k\theta_a = f(\mathcal{S}_k) + (k-1)\zeta_{k} + k\theta_a \end{align}\] where the last inequality follows from 5 and we have used \(1 \le |\mathcal{S}_k^*-\mathcal{S}_{k-1} |\leq k\). From 5 follows \(\zeta_k \le kf(\mathcal{S}_k)\), so that \[f(\mathcal{S}_k^*) \le \frac{2k-1}{k}f(\mathcal{S}_k) +k\theta_a.\] Dividing by \(\frac{k}{2k-1}=1-\frac{k-1}{2k-1}\) gives \[\begin{align} f(\mathcal{S}_k^*)\leq f(\mathcal{S}_k) +\frac{k-1}{2k-1}f(\mathcal{S}_k^*)+\frac{k^2}{2k-1}\theta_a. \end{align}\]
2.2.0.2 Multiplicative relaxation
The interchange property (4 ) ensures \[f(\mathcal{S}_{k-1} \cup \{s_j^*\}) \le \theta_m f(\mathcal{S}_k), \qquad 1 \le j \le k.\] Subtract \(f(\mathcal{S}_{k-1})\) on both sides, and use the expression for \(\zeta_k\), \[f(\mathcal{S}_{k-1} \cup \{s_j^*\} ) - f(\mathcal{S}_{k-1}) \le \theta_m f(\mathcal{S}_k) - f(\mathcal{S}_{k-1}) = \zeta_{k} + (\theta_m-1) f(\mathcal{S}_{k}) , \qquad 1\le j \le k.\] Insert this inequality into 6 , \[\begin{align} f(\mathcal{S}_k^*) &\le f(\mathcal{S}_{k-1}) +k (\theta_m-1) f(\mathcal{T}_{k}) + k \zeta_{k} \leq f(\mathcal{S}_k)\left( 1 + k(\theta_m-1) + \frac{k-1}{k} \right)\\ &= \frac{f(\mathcal{S}_k)}{k}\left(2k-1+k^2(\theta_m-1)\right), \end{align}\] where the last inequality follows from 5 . ◻
Theorem 3 implies that the solution from the 1-interchange Algorithm 2 tends to differ more from the optimal solution, as the relaxation factors \(\theta_a\) and \(\theta_m\) and the desired number \(k\) of columns increase.
2.3 Existing Work↩︎
A survey of submodular functions, their applications and extensions is given in [17]. Distributed greedy algorithms for maximizing submodular functions are analyzed in [3]. In the context of determinantal point processes, there are greedy algorithms for log-submodular functions [6] and lazy greedy algorithms [18]. Extensions to non-monotone submodular functions are analyzed in [19]–[22]. Adaptive submodularity is introduced and surveyed in [23], [24]. Randomized greedy algorithms for sensor selection are presented in [8].
3 Submodular functions in subset selection↩︎
After defining the volume of a matrix (Definition 4), we review work on column subset selection based on volume maximization (Section 3.1), show that the log-volume is a submodular function (Section 3.2), and introduce a generalized trace function (Section 3.3).
Unlike [25] and [26], the definition below is strict in that it allows a positive volume only for matrices with non-zero singular values.
Definition 4. The volume of a matrix \(\mx\in\real^{m\times n}\) with \(\rank(\mx)=n\) equals \[\vol(\mx)\equiv\sqrt{\det(\mx^T\mx)}= \prod_{i=1}^n{\sigma_i(\mx)}.\] The volume of a matrix \(\mx\in\real^{m\times n}\) with \(\rank(\mx)=m\) equals \[\vol(\mx)\equiv\sqrt{\det(\mx\mx^T)}=\prod_{i=1}^m{\sigma_i(\mx)}.\] The volume of a nonsingular matrix \(\mx\in\rmm\) equals \[\vol(\mx)=|\det(\mx)|=\prod_{i=1}^m{\sigma_i(\mx)}.\]
3.1 Existing work↩︎
Volume maximization is NP-hard. Specifically, given \(\mx\in\rmn\) with \(\rank(\mx)\geq k\), finding a \(m\times k\) column submatrix of \(\mx\) with maximal volume is a NP-hard problem [27], [28]. Moreover [29], given \(\mx\in\real^{k\times n}\) with \(\rank(\mx)=k\), there exists a universal constant such \(c>0\) so that finding a \(k\times k\) maximal-volume submatrix of \(\ma\) within a factor \(c^k\) is also NP-hard.
The Businger-Golub algorithm [11], GoVL13? computes the pivoted QR factorization of a matrix \(\mx\in\rmn\), and it corresponds to Khachiyan’s algorithm for finding the maximum subdeterminant [29]. If \(\mx_k^*\in\real^{m\times k}\) is a submatrix of \(\mx\) with maximal volume, then the output \(\mx_k\in\real^{m\times k}\) from the Businger-Golub algorithm satisfies [27] \[\vol(\mx_k)\geq \frac{1}{k!} \vol(\mx_k^*).\]
Bounds for the volume can also be derived from bounds for individual singular values.
Lemma 1. Let \(\mx \in \rmn\) and \(k\) a fixed integer with \(0<k \le \rank(\mx)\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\). Then the output \(\mx_k\in\real^{m\times k}\) from the Businger-Golub algorithm satisfies \[\frac{\vol(\mx_k^*)}{\prod_{j=1}^k 2^j \sqrt{n-j}} \le \vol(\mx_k) \le \vol(\mx_k^*).\]
Proof. The upper bound follows from the definition of \(\mx_k^*\). The bounds for individual singular values from [12] imply \[\frac{\sigma_j(\mx)}{2^j \sqrt{n-j}} \le \sigma_j(\mx_k) \le \sigma_j(\mx), \qquad 1 \le j \le k.\] Multiplication gives \[\prod_{j=1}^k \frac{\sigma_j(\mx)}{2^j\sqrt{n-j}} \le {\prod_{j=1}^k\sigma_j(\mx_k)} = {\vol(\mx_k)} .\] Singular value interlacing 2 implies \(\vol (\mx_k^*) = \prod_{j=1}^k \sigma_j(\mx_k^*) \le \prod_{j=1}^k \sigma_j(\mx)\). Insert this into the previous inequality to obtain the lower bound. ◻
Additional bounds for the case \(m\le k \le n\) are given in [30].
Next are bounds from the Gu-Eisenstat strong rank revealing QR algorithm [12].
Lemma 2. Let \(\mx \in \rmn\) and \(k\) a fixed integer with \(0<k \le \rank(\mx)\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\). Then the output \(\mx_k\in\real^{m\times k}\) from the Gu-Eisenstat [12] with relaxation factor \(f = \eta\) satisfies \[\frac{\vol(\mx_k^*)}{(1+\eta^2k(n-k))^{k/2}} \le {\vol(\mx_k)} \le \vol(\mx_k^*).\]
Algorithms for volume maximization are presented in [30], [31] when \(k\) exceeds the number of rows, with bounds on the norm of the pseudo-inverse.
The notion of submodularity for subset selection has been explored in the context of regression [2], graph Laplacian column selection [32], and minimization of the matrix residual \(\min_{\mc}\|\ma - \mc\mc^\dagger\ma\|_F\) [33]. The submodularity of the log-volume is discussed in [34], [35].
3.2 Log-volume↩︎
After presenting two examples (Examples 1 and 2), and expressing submodular functions in terms of matrices (Remark [r951]), we show that the log-volume is a submodular function (Theorem 4) and establish conditions when it is non-decreasing (Lemma 4).
Given a matrix \(\mx\in\rmn\). With the notation in Definition 2, let \(\mathcal{S}\subset\mathcal{T}\subset\{1,\ldots, n\}\), so that \(\mx_{\mathcal{S}}\) is a submatrix of \(\mx_{\mathcal{T}}\) which, in turn, is a submatrix of \(\mx\). Let \(t\in\{1,\ldots,n\}-\mathcal{T}\).
A function \(f\) on column subsets of \(\mx\) is submodular if for all submatrices \[f(\mx_{\mathcal{S}\cup\{t\}})-f(\mx_{\mathcal{S}})\geq f(\mx_{\mathcal{T}\cup\{t\}})-f(\mx_{\mathcal{T}}).\] That is, adding a column \(\mx_{\{t\}}\) to a smaller submatrix \(\mx_{\mathcal{S}}\) is more beneficial than adding \(\mx_{\{t\}}\) to a larger submatrix \(\mx_{\mathcal{T}}\).
Example 1. The volume itself is not* a submodular function on column submatrices.*
With the notation in Definition 2, let \[\mx=\begin{bmatrix} 2 & 0 & 0&1\\ 0&3 &0 &0\\ 0&0&4&0\end{bmatrix}\in\real^{3\times 4},\] \(\mathcal{S}=\{1\}\), and \(\mathcal{T}=\{1,2\}\). Then \[\begin{align} \vol(\mx_{\mathcal{S}})=2, \quad \vol(\mx_{\mathcal{S}\cup\{3\}})=8,\quad \vol(\mx_{\mathcal{T}})=6, \quad \vol(\mx_{\mathcal{T}\cup\{3\}})=24, \end{align}\] and \[\vol(\mx_{\mathcal{S}\cup\{3\}})-\vol(\mx_{\mathcal{S}})=6\not\geq 18= \vol(\mx_{\mathcal{T}\cup\{3\}})-\vol(\mx_{\mathcal{T}}).\] However, since the matrices are diagonal, \[\log{\vol(\mx_{\mathcal{S}\cup\{3\}})}-\log{\vol(\mx_{\mathcal{S}})}=\log(4)= \log{\vol(\mx_{\mathcal{T}\cup\{3\}})}- \log{\vol(\mx_{\mathcal{T}})}.\]
The log-volume of diagonal matrices is easy to determine.
Example 2. Let \(\mx\in\rmn\) be a diagonal matrix with \(\rank(\mx)=n\). Then \(\log\vol(\mx)\) is a modular* function on column submatrices of \(\mx\).*
A submodular function that satisfies the condition in Definition 3 with equality is called modular* [3], \[f(\mathcal{T}\cup \{x\}) - f(\mathcal{T})= f(\mathcal{S}\cup\{x\})-f(\mathcal{S})\] for all \(\mathcal{S}, \mathcal{T}\subset\mathcal{X}\) and \(x\in\mathcal{X}-(\mathcal{S}\cup\mathcal{T})\).*
Although the submodularity of the log-volume was proved in [34], [35] and [9] we present a clean proof for general matrices below, starting with a well known result about determinants.
Lemma 3. Let \(\mx\in\rmn\), \(\mathcal{S}\subset\{1,\ldots,n\}\), \(t\in\{1,\ldots,n\} -\mathcal{S}\), and \(\mP_{\mathcal{S}}\equiv \mx_{\mathcal{S}}(\mx_{\mathcal{S}}^T\mx_{\mathcal{S}})^{\dagger}\mx_{\mathcal{S}}\) the orthogonal projector onto \(\range(\mx_{\mathcal{S}})\). Then \[\begin{align} \vol(\mx_{\mathcal{S}\cup\{t\}}) = \vol(\mx_{\mathcal{S}})\| (\mi-\mP_{\mathcal{S}} )\mx_{\{t\}}\|. \end{align}\]
Proof. Applying Schur’s determinant formula for the generalized Schur complement [36] to \[\begin{align} \mx_{\mathcal{S}\cup\{t\}}^T\mx_{\mathcal{S}\cup\{t\}}= \begin{bmatrix}\mx_{\mathcal{S}}^T\mx_{\mathcal{S}}& \mx_{\mathcal{S}}^T\mx_{\{t\}}\\ \mx_{\{t\}}^T\mx_{\mathcal{S}} &\mx_{\{t\}}^T\mx_{\{t\}} \end{bmatrix} \end{align}\] gives \[\begin{align} \det\left(\mx_{\mathcal{S}\cup\{t\}}^T\mx_{\mathcal{S}\cup\{t\}}\right)= \det\left(\mx_{\mathcal{S}}^T\mx_{\mathcal{S}}\right) \det\left(\mx_{\{t\}}^T(\mi-\mP_{\mathcal{S}})\mx_{\{t\}}\right). \end{align}\] Hence \(\vol(\mx_{\mathcal{S}\cup\{t\}}) = \vol(\mx_{\mathcal{S}})\| (\mi-\mP_{\mathcal{S}} )\mx_{\{t\}}\|\). ◻
Theorem 4. For \(\mx\in\rmn\), \(\log\vol(\mx)\) is a submodular function on column submatrices of \(\mx\).
Proof. Since \(\log{\vol(\mx)}=\tfrac{1}{2}\log{\det(\mx^T\mx)}\), it suffices to show the submodularity of \(\log{\det(\mx^T\mx)}\).
Let \(\mathcal{S}\subset\mathcal{T}\subset\{1,\ldots,n\}\), \(t\in\{1,\ldots,n\}-\mathcal{T}\), and define orthogonal projectors \(\mP_{\mathcal{S}}\equiv \mx_{\mathcal{S}}(\mx_{\mathcal{S}}^T\mx_{\mathcal{S}})^{\dagger}\mx_{\mathcal{S}}^T\) and \(\mP_{\mathcal{T}}\equiv \mx_{\mathcal{T}}(\mx_{\mathcal{T}}^T\mx_{\mathcal{T}})^{\dagger}\mx_{\mathcal{T}}^T\). Lemma 3 implies \[\begin{align} \vol(\mx_{\mathcal{S}\cup\{t\}}) &= \vol(\mx_{\mathcal{S}})\| (\mi-\mP_{\mathcal{S}} )\mx_{\{t\}}\|\\ \vol(\mx_{\mathcal{T}\cup\{t\}}) & = \vol(\mx_{\mathcal{T}})\| (\mi-\mP_{\mathcal{T}} )\mx_{\{t\}}\|. \end{align}\] The Loewner inequality \(\mP_{\mathcal{S}} \preceq \mP_{\mathcal{T}}\) implies \[\| (\mi-\mP_{\mathcal{S}} )\mx_{\{t\}}\|^2= \mx_{\{t\}}^T (\mi - \mP_{\mathcal{S}})\mx_{\{t\} } \ge\mx_{\{t\}}^T (\mi - \mP_{\mathcal{T}})\mx_{\{t\}}= \| (\mi-\mP_{\mathcal{T}} )\mx_{\{t\}}\|^2.\] Combining everything gives \[\begin{align} \vol(\mx_{\mathcal{S}\cup\{t\}}) \geq\vol(\mx_{\mathcal{S}})\| (\mi-\mP_{\mathcal{T}}) \mx_{\{t\}}\| =\vol(\mx_{\mathcal{S}})\vol(\mx_{\mathcal{T}\cup\{t\}})/ \vol(\mx_{\mathcal{T}}). \end{align}\] At last, take logarithms and interpret \(0/0\) as \(0\). ◻
For full column-rank matrices with all singular values greater or equal to one, log-volume is a non-negative non-decreasing function.
Lemma 4. Let \(\mx\in\real^{m\times n}\) have \(\rank(\mx)=n\) and singular values \(\sigma_1(\mx)\geq \cdots \geq \sigma_n(\mx)\geq 1\). Then \(\log\vol(\mx)\geq 0\).
Furthermore, if \(\mathcal{S} \subset \mathcal{T}\subset\{1,\ldots,n\}\), then \(0 \le \log\vol(\mx_{\mathcal{S}}) \le \log\vol(\mx_{\mathcal{T}})\).
Proof. From \(\det(\mx^T\mx)=\prod_{i=1}^n{\sigma_i(\mx)^2}\geq 1\) follows \(\log(\vol(\mx))\geq 0\).
Singular value interlacing (1 ) implies \(1\leq \vol(\mx_{\mathcal{S}})\leq \vol(\mx_{\mathcal{T}})\). Therefore \(0\leq\log{\vol(\mx_{\mathcal{S}})}\leq \log{\vol(\mx_{\mathcal{T}})}\). ◻
3.3 A generalized trace function↩︎
Since log-volume is a non-decreasing function only for matrices with full column-rank and with singular values at least \(1\), we need an alternative criterion that applies to all matrices. To this end, we define a generalized trace function (Definition 5) that is non-decreasing and non-negative for general matrices, and show that it is submodular (Theorem 5, Corollary 1).
Let \(\mx\in\rmn\) and a (partial) pivoted QR decomposition GoVL13? \[\label{eqn:xsqr} \mx\mP=\begin{bmatrix} \mx_1 & \mx_2\end{bmatrix}=\begin{bmatrix}\mq_1 & \mq_2\end{bmatrix} \begin{bmatrix} \mr_{11} & \mr_{12}\\ \vzero & \mr_{22}\end{bmatrix},\tag{7}\] where \(1\leq k\leq \rank(\mx)\), \(\mP\in\rnn\) is a permutation matrix, \(\mx_1\in\real^{m\times k}\), \(\mq_1\in\real^{m\times k}\) has orthonormal columns, and \(\mr_{11}\in\real^{k\times k}\) is nonsingular upper triangular with positive diagonal elements. We will show that \(\trace(\mr_{11})=\sum_{j=1}^k(\mr_{11})_{jj}\) is a non-negative, non-deceasing, submodular function.
3.3.0.1 Preparation
Since the diagonal elements of \(\mr_{11}\) depend on each other, we need to identify this dependence clearly and fix subsets of columns of \(\mx\) ahead of time, so that we can define a function on subsets of columns of \(\mx\).
This is done based on the modified Gram-Schmidt algorithm GoVL13?. Assume that pivoting is performed ahead of time, and distinguish the columns \(\mx_1=\begin{bmatrix}\vx_1& \cdots &\vx_k\end{bmatrix}\). Focusing on the computation of the diagonal elements of \(\mr_{11}\), and initializing \(\mP_0\equiv \vzero_{m\times m}\) gives \[\begin{align} \vy_j\equiv (\mi-\mP_{j-1})\vx_j, \quad \mP_j\equiv \mP_{j-1}+\vy_j\vy_j^{\dagger}, \quad (\mr_{11})_{jj}\equiv \|\vy_j\|, \qquad 1\leq j\leq k. \end{align}\] The projected columns are orthogonal, \(\vy_j^T\vy_i = 0\) for \(i\neq j\). Hence \(\begin{bmatrix}\vy_1 &\ldots &\vy_j\end{bmatrix}\) represents an orthogonal basis for \(\range\begin{bmatrix}\vx_1 & \ldots & \vx_j\end{bmatrix}\), \(1\leq j\leq k\).
3.3.0.2 Functions on matrix columns
Let \(\mathcal{S}=\{s_1,\ldots, s_k\}\subset\{1,\ldots, n\}\) be the set of column indices of \(\mx_1\), with nested subsets \[\begin{align} \mathcal{S}_0\equiv\emptyset,\qquad \mathcal{S}_j=\{s_1,\ldots, s_j\}, \qquad 1\leq j\leq k-1. \end{align}\] Then we can write \(\mx_{\mathcal{S}}= \mx_1\), with leading column submatrices \(\mx_{\mathcal{S}_j}\in\real^{m\times j}\), \(1\leq j\leq k-1\). Define the orthogonal projectors onto \(\range(\mx_{\mathcal{S}_{j}})\), \[\begin{align} \mP_0\equiv\vzero_{m\times m}, \qquad \mP_{j} \equiv \mx_{\mathcal{S}_{j}}\mx_{\mathcal{S}_{j}}^\dagger, \qquad 1\leq j\leq k-1. \end{align}\] Lemma 3 implies \[\begin{align} \vol(\mx_{\mathcal{S}_j}) = \vol(\mx_{\mathcal{S}_{j-1}})\| (\mi-\mP_{j-1})\mx_{\{s_j\}}\|, \qquad 1\leq j\leq k-1. \end{align}\]
This motivates the following definition.
Definition 5. Let \(g\) be a non-negative, non-decreasing function with \(g(0) = 0\), and define \[\trace_g(\mx_{\mathcal{S}}) \equiv \sum_{j=1}^k g\left(\| (\mi - \mP_{j-1}) \mx_{\{s_j\}}\|\right)= \sum_{j=1}^k{g((\mr_{11})_{jj})}.\] An alternative expression is \[\trace_g(\mx_{\mathcal{S}}) = \sum_{j=1}^kg\left(\frac{\vol(\mx_{\mathcal{S}_j})}{\vol(\mx_{\mathcal{S}_{j-1}})}\right) = \sum_{j=1}^kg\left(\sqrt{\frac{\det(\mx_{\mathcal{S}_j}^\top\mx_{\mathcal{S}_j})}{\det(\mx_{\mathcal{S}_{j-1}}^\top\mx_{\mathcal{S}_{j-1}})}}\right),\] where we interpret \(0/0 = 0\).
The function \(\trace_g\) equals the log-volume if \(g(x) = \log(x)\).
Theorem 5. The function \(\trace_g\) in Definition 5 is non-negative, non-decreasing, and submodular on column subsets of \(\mx\in\rmn\) with \(\trace_g(\mx_{\emptyset}) = 0\).
Proof. Since \(g\) is non-negative and non-decreasing, so is \(\trace_g\). To show that \(\trace_g\) is submodular, let \(\mathcal{S}\subset\mathcal{T}\subset \{1,\ldots, n\}\), and \(j \in \{1,\ldots,n\} - \mathcal{T}\). Since \(g\) is nonnegative with \(g(0) = 0\), it follows that \(\trace_g(\mx_{\emptyset}) = 0\). Define the orthogonal projectors \(\mP_{\mathcal{S}} \equiv\mx_{\mathcal{S}}\mx_{\mathcal{S}}^\dagger\) and \(\mP_{\mathcal{T}} \equiv \mx_{\mathcal{T}}\mx_{\mathcal{T}}^\dagger\). The Loewner inequality \(\mP_{\mathcal{S}} \preceq \mP_{\mathcal{T}}\) implies \[\begin{align} \trace_g(\mx_{\mathcal{T}\cup \{j\}}) - \trace_g(\mx_{\mathcal{T}}) = & g\left(\|(\mi - \mP_{\mathcal{T}})\mx_{\{j\}}\|\right) \le g\left( \|(\mi - \mP_{\mathcal{S}})\mx_{\{j\}}\|\right)\\ =& \trace_g(\mx_{\mathcal{S}\cup \{j\}}) - \trace_g(\mx_{\mathcal{S}}).\end{align}\] ◻
Corollary 1. Let \(\mx\in\rmn\) have a pivoted QR decomposition \(\mx\mP=\mq\mr\) where \(\mP\in\rnn\) is permutation matrix, \(\mq\in\rmm\) is an orthogonal matrix, and \(\mr\in\rmn\) is upper trapezoidal with non-negative diagonal elements.
The function \(\trace(\mr)\) is submodular, non-decreasing, and non-negative.
Proof. This follows from Theorem 5 with \(g(x)=x\). ◻
4 A greedy algorithm for column subset selection↩︎
We show that the Businger-Golub QR factorization (Algorithm 3) is a greedy algorithm in the sense of Algorithm 1 for maximizing two functions: the log-volume of a full column-rank matrix (Section 4.1), and the trace of the upper triangular matrix in a QR decomposition of a general matrix (Section 4.2).
![Figure 3: Conceptual version of Businger-Golub QR [11], GoVL13?](/article/2607.13823/_pdflatex/qhngkdsa.png "fig:")
Algorithm 3 views the Businger-Golub QR factorization [11], GoVL13? from a modified Gram-Schmidt perspective, with a focus on selecting \(k\) columns rather than computing a QR factorization.
4.1 Log-volume for full column-rank matrices↩︎
We show that Algorithm 3 is a greedy algorithm for maximizing the log-volume of a full column-rank matrix with sufficiently large singular values (Theorem 6). For full-column matrices with general singular singular values, we give bounds for scaled matrices (Corollary 2), and discuss the ramifications of the scaling (Remark [r952]).
Theorem 6 below shows that Algorithm 3 is a greedy algorithm and it approximates the maximal log-volume with an error of at most 37 percent when applied to full column-rank matrices with sufficiently large singular values.
Theorem 6. Let \(\mx\in\rmn\) have \(\rank(\mx)=n\) and singular values \(\sigma_1(\mx)\geq \cdots \geq\sigma_n(\mx)\geq 1\). Let \(k\) be a fixed integer with \(0<k<n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
Algorithm 3 is a greedy algorithm in the sense of Algorithm 1 for maximizing the log-volume, and its output matrix \(\mx_k\equiv\mx_{\mathcal{S}_k}\) satisfies \[\label{e95t2a} 0\leq\log{\vol(\mx_k^*)}-\log{\vol(\mx_k)}\leq (1-\tfrac{1}{k})^k\;\log{\vol(\mx_k^*)}\qquad{(1)}\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e}< .37\).
The upper bound in (?? ) is tight for matrices \(\mx\) with orthonormal columns.
Proof. Theorem 4 and Lemma 4 imply that log-volume is a non-negative, non-decreasing submodular function on column submatrices of \(\mx\in\rmn\). This justifies that we can initialize the volume of the empty matrix as \(\vol(\emptyset)=1\), so that \(\log(\vol(\emptyset))=0\).
Algorithm 3 selects \(k\) columns \(\mx_{\mathcal{S}_k}\) of \(\mx\). Let \(\mx_{\mathcal{S}_k}=\mq\mr\) be a QR decomposition where \(\mq\in\real^{m\times k}\) has orthonormal columns and \(\mr\in\real^{k\times k}\) is nonsingular, with non-increasing diagonal elements and \(r_{11}\geq \ldots \geq r_{kk}> 0\).
First we show that Algorithm 3 is a greedy algorithm in the sense of Algorithm 1.
4.1.0.1 Induction basis
For \(j=1\), let \(\mx_{\mathcal{S}_1}\in\real^m\) be a column of \(\mx\) with QR decomposition is \(\mx_{\mathcal{S}_1}=\mq_1r_{11}\) where \(\mq_1\in\real^m\) has \(\|\mq_1\|=1\) and \(r_{11}\geq 0\). Algorithm 3 selects the column with largest two-norm, so that among all columns of \(\mx\), \[r_{11}^2=\|\mx_{\mathcal{S}_1}\|^2= \det(\mx_{\mathcal{S}_1}^T\mx_{\mathcal{S}_1})= \vol(\mx_{\mathcal{S}_1})^2\] is maximal. Let \(\mP_1\equiv \mx_{\mathcal{S}_1}\mx_{\mathcal{S}_1}^{\dagger}=\mq_1\mq_1^T\) be the orthogonal projector onto \(\range(\mx_{\mathcal{S}_{1}})\).
4.1.0.2 Induction hypothesis
Assume that for \(1\leq j-1<n\), Algorithm 3 is a greedy algorithm and has selected the \(j-1\) linearly independent columns \(\mx_{\mathcal{S}_{j-1}}=\mq_{j-1} \mr_{j-1}\) where \(\mq_{j-1}\in\real^{m\times (j-1)}\) has orthonormal columns, \(\mr_{j-1}\in\real^{(j-1)\times (j-1)}\) is nonsingular upper triangular, and \[r_{11}^2 \cdots r_{j-1,j-1}^2=\det(\mr_{j-1}^T\mr_{j-1})= \det(\mx_{\mathcal{S}_{j-1}}^T\mx_{\mathcal{S}_{j-1}})= \vol(\mx_{\mathcal{S}_{j-1}})^2.\] Let \(\mP_{j-1}\equiv\mx_{\mathcal{S}_{j-1}}\mx_{\mathcal{S}_{j-1}}^{\dagger}=\mq_{j-1}\mq_{j-1}^T\) be the orthogonal projector onto \(\range(\mx_{\mathcal{S}_{j-1}})\).
4.1.0.3 Induction step
Let \(\mx_{\mathcal{S}_j}\equiv\begin{bmatrix} \mx_{\mathcal{S}_{j-1}} & \vx_t\end{bmatrix}=\mq_j\mr_j\) where \[\mq_j=\begin{bmatrix} \mq_{j-1} & \vq_j\end{bmatrix}\in\real^{m\times j}, \qquad \mr_j\equiv \begin{bmatrix} \mr_{j-1} & *\\ \vzero & r_{jj}\end{bmatrix}\in\real^{j\times j},\] \(\mq_j\) has orthonormal columns, and \(\vq_jr_{jj}=(\mi-\mP_{j-1})\vx_j\). Among all columns \(i\in\mathcal{T}_{j-1}\), Algorithm 3 selects the column \(\vx_t\) with maximal \(r_{jj}=\|(\mi-\mP_{j-1})\vx_t\|\). With \(\vy_j\equiv (\mi-\mP_{j-1})\vx_t\), set \[\mP_j\equiv \mP_{j-1}+\vy_j\vy_j^{\dagger}=\mq_{j-1}\mq_{j-1}^T +\vq_j\vq_j^T=\mq_j\mq_j^T=\mx_{\mathcal{S}_j}\mx_{\mathcal{S}_j}^{\dagger}.\] Because \(\mr_j\) is square, \[\begin{align} \det(\mx_{\mathcal{S}_j}^T\mx_{\mathcal{S}_j})&=\det(\mr_j^T\mr_j)=\det(\mr_j)^2 =\det(\mr_{j-1})^2\>r_{jj}^2 =\det(\mx_{\mathcal{S}_{j-1}}^T\mx_{\mathcal{S}_{j-1}})\> r_{jj}^2. \end{align}\] Since Algorithm 3 chooses a column \(\vx_t\) with largest \(r_{jj}\), it is a greedy algorithm in the sense of Algorithm 1. Now apply Theorem 1 to the log-volume.
The upper bound in (?? ) is tight for matrices \(\mx\) with orthonormal columns, since \(\log{\vol(\mb)}=0\) for any submatrix \(\mb\in\real^{m\times k}\) of \(\mx\). ◻
For \(k=1\), \(\mx_k=\mx_k^*\) and the bound (?? ) is zero.
If \(\mx\) has singular values less than one, we multiply by a scalar to make all singular values are greater or equal to 1. This kind of manipulation is not always recommended [4] as it makes the relative error not unique. However, since Algorithm 3 is invariant under scalar multiplication, we feel that this is justified.
Corollary 2. Let \(\mx\in\rmn\) have \(\rank(\mx)=n\), and let \(k\) be a fixed integer with \(0<k<n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\). The output matrix \(\mx_k\equiv\mx_{\mathcal{S}_k}\) from Algorithm 3 satisfies \[0\leq \log{\vol(\mx_k^*)}-\log{\vol(\mx_k)}\leq (1-\tfrac{1}{k})^k\;\log{\vol(\mx_k^*/\sigma_n(\mx))},\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e}< .37\).
Proof. To scale matrices so that the log-volume is a non-negative non-decreasing function, let the singular values \(\sigma_1(\mx)\geq \cdots \geq \sigma_n(\mx)\geq\tfrac{1}{\alpha}>0\) for some \(\alpha>0\). Since \(\sigma_1(\alpha\mx)\geq\cdots\geq \sigma_n(\alpha\mx)\geq 1\), Lemma 4 implies that \(\log\vol(\alpha\mx)\geq 0\) and non-decreasing, that is: If \(\mathcal{S}\subset\mathcal{T}\subset\{1,\ldots, n\}\), then \(0\leq\log{\vol(\alpha\mx_{\mathcal{S}})}\leq \log{\vol(\alpha\mx_{\mathcal{T}})}\).
Furthermore, Algorithm 3 is invariant under multiplication by a scalar \(\alpha>0\) and so is the absolute error, \[\begin{align} \log{\vol(\alpha\mx_k^*)}-\log{\vol(\alpha\mx_k)}&=(k\log{\alpha}+\log{\vol(\mx_k^*)})-(k\log{\alpha}+\log{\vol(\mx_k)})\\ &= \log{\vol(\mx_k^*)}-\log{\vol(\mx_k)}. \end{align}\] Choose \(\alpha=1/\sigma_n(\mx)\). This ensures that \(\log{\vol(\mb)}\geq 0\) for any submatrix \(\mb\in\real^{m\times k}\) of \(\mx\). At last apply Theorem 6 to the scaled matrix \(\mx/\sigma_n(\mx)\). ◻
If all singular values of \(\mx\) are the same, i.e. \(\sigma_1(\mx)=\sigma_n(\mx)\), then any submatrix \(\mb\in\real^{m\times k}\) of \(\mx\) maximizes the volume, and satisfies \(\log{\vol(\mx_k^*)}=\log{\vol(\mb)}\). Thus, Algorithm 3 terminates with an optimal solution in \(k\) steps. The scaling in Corollary 2 ensures a zero error bound, with \(\log{\vol(\mx_k^*/\sigma_n(\mx))}=0\), and thus the tightness of the bound.
If \(\sigma_n(\mx)<1\), the scaling in Corollary 2 guarantees that the log-volume of the scaled matrix is non-negative. If \(\sigma_n(\mx)>1\), then the scaling lowers the bound on the absolute error.
4.2 Trace for general matrices↩︎
For general matrices that do not necessarily have full column rank, we show that Algorithm 3 maximizes the trace of the upper triangular matrix in the QR decomposition. Theorem 7 below uses the notation for upper triangular matrices from Definition 1.
Theorem 7. Let \(\mx\in\rmn\) and let \(k\) be a fixed integer with \(0<k\leq\rank(\mx)\). Let \(\mx_k^*\in\real^{m\times k}\) be a matrix that maximizes \(\trace(\mathcal{R}(\mb))\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
Algorithm 3 is a greedy algorithm in the sense of Algorithm 1 for maximzing the trace, and its output matrix \(\mx_k\equiv\mx_{\mathcal{S}_k}\) satisfies \[0\leq \trace(\mathcal{R}(\mx_k^*))-\trace(\mathcal{R}(\mx_k)) \leq (1-\tfrac{1}{k})^k\,\trace(\mathcal{R}(\mx_k^*)),\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e} <.37\).
Thus, for general matrices, Algorithm 3 selects columns \(\mx_k\) whose upper triangular matrix \(\mathcal{R}(\mx_k)\) maximizes the trace with a relative error of at most 37 percent.
5 A 1-interchange algorithm for column subset selection↩︎
We show that the Gu-Eisenstat QR algorithm (Algorithm 4) is a 1-interchange algorithm in the sense of Algorithm 2 for maximizing two functions: the log-volume of a full column-rank matrix (Section 5.1), and the trace of the upper triangular matrix in a QR factorization of a general matrix (Section 5.2).
![Figure 4: Conceptual version of Gu-Eisenstat QR [12]](/article/2607.13823/_pdflatex/honjzgau.png "fig:")
Algorithm 4 represents a conceptual version of the Gu-Eisenstat QR factorization [12], with a focus on selecting \(k\) columns rather than computing a QR factorization. Algorithm 4 starts with an initial set of \(k\) columns that are linearly independent. Without this initialization, the first swap could still produce a rank deficient submatrix with zero volume, thus the algorithm would make no progress and terminate. The relaxation factor \(\eta=1\) represents an ideal version, where Algorithm 4 is continued until the volume of the leading \(k\) columns stops increasing.
5.1 Log-volume for full column-rank matrices↩︎
We show that Algorithm 4 is a 1-interchange algorithm for maximizing the log-volume of a full column-rank matrix in two ways: in the sense of Algorithm 2 if the relaxation factor \(\eta=1\) (Theorem 8 and Corollary 3); and in the sense of Theorem 3 if the relaxation factor \(\eta>1\) (Theorem 9) in which case we can bound the number of swaps (Theorem 10).
Theorem 8 below shows that Algorithm 4 with relaxation factor \(\eta=1\) applied to a full column-rank matrix with sufficiently large singular values is an instance of the \(R\)-interchange algorithm with \(R=1\) [1], that approximates the maximal log-volume with an error of at most 50 percent.
Theorem 8. Let \(\mx\in\rmn\) have \(\rank(\mx)=n\) with singular values \(\sigma_1(\mx)\geq \cdots \geq \sigma_n(\mx)\geq 1\), and let \(k\) be a fixed integer with \(0<k< n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
Algorithm 4 with relaxation factor \(\eta=1\) is a 1-interchange algorithm in the sense of Algorithm 2 for maximizing the log-volume, and its output matrix \(\mx_k\equiv \mx_{\mathcal{S}}\) satisfies \[\label{e95t3a} 0\leq \log{\vol(\mx_k^*)}-\log{\vol(\mx_k)}\leq \frac{k-1}{2k-1}\;\log{\vol(\mx_k^*)},\qquad{(2)}\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\).
The upper bound in (?? ) is tight for matrices \(\mx\) with orthonormal columns.
Proof. This follows from Theorem 2 applied to the function \(f(\mathcal{S}) = \log\vol(\mx_{\mathcal{S}})\) and Lemma 4. ◻
For \(k=1\), \(\mx_k=\mx_k^*\) and the bound (?? ) is equal to zero.
Corollary 3 presents a bound, based on scaling, for full column-rank matrices with general singular values, and the observations from Remark [r952] apply here as well.
Corollary 3. Let \(\mx\in\rmn\) with \(\rank(\mx)=n\), and let \(k\) be a fixed integer with \(0<k<n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\). The output matrix \(\mx_k\equiv\mx_{\mathcal{S}}\) from Algorithm 4 with relaxation factor \(\eta=1\) satisfies \[0\leq \log{\vol(\mx_k^*)}-\log{\vol(\mx_k)}\leq \frac{k-1}{2k-1}\;\log{\vol(\mx_k^*/\sigma_n(\mx))},\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\).
Proof. This follows from Theorem 8 applied to the scaled matrix \(\mx/\sigma_n(\mx)\). ◻
Theorem 9 shows that Algorithm 4 with relaxation factor \(\eta > 1\) is a 1-interchange algorithm for maximizing the log-volume, and implies that the larger \(\eta\), the larger the error in the log-volume.
Theorem 9. Let \(\mx\in\rmn\) have \(\rank(\mx)=n\), and let \(k\) be a fixed integer with \(0<k<n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
Algorithm 4 with relaxation factor \(\eta>1\) is a 1-interchange algorithm in the sense of Theorem 3, and its output matrix \(\mx_k\equiv\mx_{\mathcal{S}}\) satisfies \[0\leq \log{\vol(\mx_k^*)}-\log{\vol(\mx_k)} \leq \frac{k-1}{2k-1}\log\vol(\mx_k^*/\sigma_n(\mx)) + \frac{k^2}{2k-1}\,\log(\eta),\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\) and \(\frac{k^2}{2k-1}\leq k\).
Proof. This follows from Theorem 3 with additive relaxation factor \(\theta_a = \log \eta\). ◻
An initialization of Algorithm 4 with Algorithm 3 is recommended in [37], to reduce the number of swaps. However [1] cautions that a greedy initialization may not improve worst case behavior. Nevertheless, Theorem 10 below bounds the number of swaps.
Theorem 10. Let \(\mx\in\rmn\) have \(\rank(\mx)=n\), and let \(k\) be a fixed integer with \(0<k<n\). Let \(\mx_k^*\in\real^{m\times k}\) be a column submatrix of \(\mx\) that maximizes \(\log{\vol(\mb)}\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
If line 1 of Algorithm 4 with relaxation factor \(\eta>1\) is initialized with the output matrix \(\mx_{\mathcal{S}_k}\) from Algorithm 3, the maximal number of swaps is \[\label{eqn:gueisbound} q \le \frac{\log(\vol(\mx_k^*/\sigma_n(\mx))}{e\,\log(\eta)} \le \frac{k}{e}\log_\eta(\kappa_2(\mx)).\qquad{(3)}\]
Proof. Let \(\mx_k\) be the output from Algorithm 3. Corollary 2 implies \[\begin{align} \log\left(\vol(\mx_k^*)/\vol(\mx_k)\right)\leq \log\left(\vol(\mx_k^*)^{1/e}/\sigma_n(\mx)^{k/e}\right) \end{align}\] Since the logarithm is a monotone function, \[{\vol(\mx_k^*)^\gamma}{\sigma_n(\mx)^{k/e}} \le \vol(\mx_k) \le \vol(\mx_k^*), \qquad \gamma \equiv 1-1/e.\] Each swap of Algorithm 4 improves the volume by a factor \(\eta > 1\). After \(q\) swaps, the improvement is \(\eta^q\), hence \[\eta^q \, {\vol(\mx_k^*)^\gamma}{\sigma_n(\mx)^{k/e}} \le \eta^q\,\vol(\mx_k)\leq \vol(\mx_k^*).\] Rearranging the lower and upper bounds gives \[\eta^q \le \vol(\mx_k^*)^{1/e} /\sigma_n(\mx)^{k/e} = \vol(\mx_k^*/\sigma_n(\mx))^{1/e}.\] Thus, the maximal number of swaps is \[q \le \frac{1}{e}\frac{\log(\vol(\mx_k^*/\sigma_n(\mx))}{\log(\eta)} \le \frac{k}{e}\log_\eta(\kappa_2(\mx)).\]\(\qquad\) ◻
Theorem 10 implies that the number of swaps decreases as the relaxation factor increases. In contrast to the bound \(q \le k \log_\eta(\sqrt{n})\) from [12] and [37], the bound in Theorem 10 is independent of \(n\) and tighter if \(\kappa_2(\mx)^{1/e} < \sqrt{n}\); that is, if the matrix is well-conditioned. Combining the bounds shows that the number of swaps cannot exceed \[q \le k \min\{ \log_\eta(\sqrt{n}), \log_\eta(\kappa_2(\mx))/e\}.\]
5.2 Trace for general matrices↩︎
For general matrices that do not necessarily have full column-rank, we show that Algorithm 4 maximizes the trace of the upper triangular matrix in a QR decomposition. Theorem 11 below uses the notation for upper triangular matrices from Definition 1.
Theorem 11. Let \(\mx\in\rmn\) and let \(k\) be a fixed integer with \(0<k\leq\rank(\mx)\). Let \(\mx_k^*\in\real^{m\times k}\) be a matrix that maximizes \(\trace(\mathcal{R}(\mb))\) among all submatrices \(\mb\in\real^{m\times k}\) of \(\mx\).
Algorithm 4 is a 1-interchange algorithm in the sense of Algorithm 2 for maximizing the trace, and its output matrix \(\mx_k\equiv\mx_{\mathcal{S}}\) satisfies \[0 \le \trace(\mathcal{R}(\mx_k^*)) - \trace(\mathcal{R}(\mx_k)) \le \frac{k-1}{2k-1}\,\trace(\mathcal{R}(\mx_k^*)),\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\).
Thus, for general matrices, Algorithm 4 selects columns \(\mx_k\) whose upper triangular matrix \(\mathcal{R}(\mx_k)\) maximizes the trace with a relative error of at most 50 percent. A comparison with Theorem 7 shows that, in this sense, Algorithm 4 is less accurate than Algorithm 3.
6 Greedy and 1-interchange algorithms for submatrix selection↩︎
We analyze algorithms for selecting principal submatrices of maximal volume from a symmetric positive-definite matrix.
After reviewing existing work (Section 6.1), we show that pivoted Cholesky (Algorithm 5) is a greedy algorithm in the sense of Algorithm 1 for maximizing two functions: the log-volume of a positive definite matrix (Section 6.2), and the trace of the squared Cholesky factor of a positive semi-definite matrix (Section 6.3). Then we show that strong rank-revealing Cholesky (Algorithm 6) is a 1-interchange algorithm in the sense of Algorithm 2 for maximizing the log-volume of a positive definite matrix (Section 6.4).
6.1 Existing work↩︎
Finding a \(k\times k\) submatrix with maximal volume \(\ma\in\rnn\) is NP-hard [27], [38].
For matrices \(\ma\in\rnn\) with \(\rank(\ma)\geq k\), a \(k\times k\) submatrix with largest volume is attained by a principal submatrix if \(\ma\) is diagonally dominant [39], or if \(\ma\) is symmetric positive semi-definite [39]. Low-rank approximations have been implemented by pivoted Cholesky factorizations [40], and randomized pivoted Cholesky factorizations [41], [42].
The pivoted Cholesky algorithm, in Algorithm 5, coincides with adaptive cross approximation [43] when applied to symmetric positive semi-definite matrices [40]. ACA computes a low-rank approximation, but can be viewed as a greedy method for volume maximization, and amounts to Gaussian elimination process with rook pivoting [43], [38]. It seems to be the same as the fast Greedy MAP inference [5]. Other algorithms for volume maximization have been proposed in [37], [44]. Log-volume maximization via pivoted Cholesky is implemented as fast Greedy MAP inference [5].
![Figure 5: Conceptual version of pivoted Cholesky [14], GoVL13?](/article/2607.13823/_pdflatex/aedvgbzh.png "fig:")
6.2 Greedy algorithm for positive definite matrices↩︎
After relating Algorithm 5 to Algorithm 3 (Theorem 12), we show that Algorithm 5 is a greedy algorithm for maximizing the log-volume of symmetric positive definite matrices (Theorem 13).
Algorithm 5 presents the conceptual version of the pivoted Cholesky algorithm GoVL13?, [14], with a focus on selecting a \(k\times k\) principal submatrix rather than computing a Cholesky factorization.
Theorem 12 below shows that Algorithm 3 applied to a Cholesky factor \(\mx\) of \(\ma\) implements Algorithm 5. The Cholesky factor \(\mx\) does not have to be a triangular matrix.
Theorem 12. Let \(\ma=\mx^T\mx\) where \(\mx\in\rnn\), and let \(k\) be a fixed integer with \(0<k\leq \rank(\mx)\). Let \(\mx_{\mathcal{S}_k}\in\real^{n\times k}\) be the output matrix of Algorithm 3. Then the output matrix of Algorithm 5 satisfies \[\ma[\mathcal{S}_k]=\mx_{\mathcal{S}_k}^T\mx_{\mathcal{S}_k}.\]
Proof. With an induction proof over the dimension of the principal submatrices, we show that Algorithm 3 selects the same indices, apart from ties, as Algorithm 5.
6.2.0.1 Induction basis
For \(j=1\), Algorithm 3 selects a column \(\mx_{\mathcal{S}_1}\) of \(\mx\) with largest two-norm, so that among all columns of \(\mx\), \[\|\mx_{\mathcal{S}_1}\|^2= \det(\mx_{\mathcal{S}_1}^T\mx_{\mathcal{S}_1})= \vol(\ma[\mathcal{S}_1])\] is maximal. The corresponding column and diagonal element of \(\ma=\mx^T\mx\) are \(\va_t\equiv \mx^T\mx_{\mathcal{S}_1}\), and \(a_{tt}=\mx_{\mathcal{S}_1}^T\mx_{\mathcal{S}_1}\), respectively. Then \[\begin{align} \mx^T(\mi-\mP_1)\mx= \ma-\mx^T\mx_{\mathcal{S}_1}(\mx_{\mathcal{S}_1}^T\mx_{\mathcal{S}_1})^{-1}\mx_{\mathcal{S}_1}^T\mx =\ma-\va_t a_{tt}^{-1} \va_t^T=\ma_1, \end{align}\] where \(\mP_1\) is the orthogonal projector onto \(\range(\mx_{\mathcal{S}_1})\). With \(\mx=\begin{bmatrix}\mx_{\mathcal{S}_1}&\mx_{\mathcal{T}_1}\end{bmatrix}\) this gives \[\ma_1=\begin{bmatrix} 0 & \vzero \\ \vzero & \mx_{\mathcal{T}_1}^T(\mi-\mP_1)\mx_{\mathcal{T}_1}. \end{bmatrix}\]
6.2.0.2 Induction hypothesis
Assume that for \(1\leq j-1<n\), Algorithm 3 has selected the \(j-1\) linearly independent columns \(\mx_{\mathcal{S}_{j-1}}\), and \[\det(\mx_{\mathcal{S}_{j-1}}^T\mx_{\mathcal{S}_{j-1}})= \vol(\ma[\mathcal{S}_{j-1}]).\] Let \(\mP_{j-1}\equiv\mx_{\mathcal{S}_{j-1}}\mx_{\mathcal{S}_{j-1}}^{\dagger}=\mq_{j-1}\mq_{j-1}^T\) be the orthogonal projector onto \(\range(\mx_{\mathcal{S}_{j-1}})\). Then \[\begin{align} \mx^T(\mi-\mP_{j-1})\mx&= \ma-\mx^T\mx_{\mathcal{S}_{j-1}}(\mx_{\mathcal{S}_{j-1}}^T\mx_{\mathcal{S}_{j-1}})^{-1}\mx_{\mathcal{S}_{j-1}}^T\mx\\ &=\begin{bmatrix} \vzero_{(j-1)\times (j-1)} & \vzero \\ \vzero & \mx_{\mathcal{T}_{j-1}}^T(\mi-\mP_{j-1})\mx_{\mathcal{T}_{j-1}}\end{bmatrix}=\ma_{j-1}. \end{align}\]
6.2.0.3 Induction step
Let \(\mx_{\mathcal{S}_j}\equiv\begin{bmatrix} \mx_{\mathcal{S}_{j-1}} & \vx_t\end{bmatrix}\), where Algorithm 3 selects a column \(\vx_t\in\mx_{\mathcal{T}_{j-1}}\) with maximal \(\|(\mi-\mP_{j-1})\vx_t\|\). With \(\vy_j\equiv (\mi-\mP_{j-1})\vx_t\), the corresponding column \(\va_t\) and diagonal element \(a_{tt}\) of \(\ma_{j-1}\) are respectively \[\va_t=\begin{bmatrix}\vzero_{(j-1)\times 1}\\ \mx_{\mathcal{T}_{j-1}}^T(\mi-\mP_{j-1})\vy_j\end{bmatrix}\quad \text{and} \quad a_{tt}=\vy_j^T\vy_j.\] Update \(\mx_{\mathcal{T}_{j-1}}\equiv\begin{bmatrix} \vx_t& \mx_{\mathcal{T}_{j}} \end{bmatrix}\) and abbreviate \[(\mi-\mP_{j-1})\mx_{\mathcal{T}_{j-1}}= \begin{bmatrix} \vy_j & \my_j\end{bmatrix}\qquad\text{where}\quad \my_j\equiv (\mi-\mP_{j-1})\mx_{\mathcal{T}_j}.\] Then the downdate is \[\ma_j=\ma_{j-1}-\va_t a_{tt}^{-1}\va^T= \begin{bmatrix} \vzero_{j\times j} & \vzero \\ \vzero & \mz_j\end{bmatrix}\] where \[\begin{align} \mz_j&\equiv \my_j^T\my_j-\my_j^T\vy_j(\vy_j^T\vy_j)^{-1}\vy_j^T\my_j= \my_j^T\left(\mi-\vy_j\vy_j^{\dagger}\right)\my_j =\mx_{\mathcal{T}_j}^T(\mi-\mP_j)\mx_{\mathcal{T}_j} \end{align}\] and \(\mP_j\equiv \mP_{j-1}+\vy_j\vy_j^{\dagger} =\mx_{\mathcal{S}_j}\mx_{\mathcal{S}_j}^{\dagger}\) is the orthogonal projector onto \(\range(\mx_{\mathcal{S}_j})\). Thus, \[\ma_j=\begin{bmatrix} \vzero_{j\times j} & \vzero \\ \vzero & \mx_{\mathcal{T}_{j}}^T(\mi-\mP_{j})\mx_{\mathcal{T}_{j}}\end{bmatrix}.\] The proof of Theorem 6 implies \[\begin{align} \vol(\ma[\mathcal{S}_j])=\det(\mx_{\mathcal{S}_j}^T\mx_{\mathcal{S}_j})& =\det(\mx_{\mathcal{S}_{j-1}}^T\mx_{\mathcal{S}_{j-1}})\> \|\vy_j\|^2=\vol[(\ma[\mathcal{S}_{j-1}])\|\vy_j\|^2. \end{align}\] ◻
Theorem 13 below shows that Algorithm 5 is a greedy algorithm for maximizing the log-volume of symmetric positive definite matrices.
Theorem 13. Let \(\ma\in\rnn\) be symmetric positive definite, and let \(k\) be a fixed integer with \(0<k<n\). Let \(\ma_k^*\in\real^{k\times k}\) be a principal submatrix of \(\ma\) that maximizes \(\log\vol(\mb)\) among all principal submatrices \(\mb\in\real^{k\times k}\) of \(\ma\).
Then Algorithm 5 is a greedy algorithm in the sense of Algorithm 1 for maximizing the log-volume and its output matrix \(\ma_k\equiv\ma[\mathcal{S}_k]\) satisfies \[\label{e95piva} 0 \le {\log\vol(\ma_k^*)-\log\vol(\ma_k)} \le \left(1-\frac{1}{k}\right)^k \log\vol(\ma_k^*/\lambda_n(\ma)),\qquad{(4)}\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e} <.37\).
The upper bound in (?? ) is tight for matrices \(\ma\) that are non-zero multiples of the identity matrix.
Symmetric positive definite matrices \(\ma\) arise in applications, such as sensor placement for Gaussian processes, where \(\ma = \mi + \sigma_n^{-2}\mk\), with \(\mk\) being a positive semidefinite kernel matrix and \(\sigma_n\) the noise standard deviation [3], [45], [46]. In such applications, Theorem 13 applies with \(\lambda_n(\ma) \ge 1\).
6.3 Trace for positive semi-definite matrices↩︎
For matrices that are semi-definite, log-volume can be decreasing. Thus, similar to Definition 5, we introduce a submodular function \(f\). Let \(\mm\) be a symmetric positive semidefinite matrix, partitioned as [36] \[\mm = \bmat{\mb & \mc\\ \mc^T & \md} = \bmat{ \mi \\ \mc^T \mb^\dagger & \mi } \bmat{\mb \\ & \md - \mc^T \mb^\dagger \mc} \bmat{\mi & \mb^\dagger\mc \\ & \mi} .\] We call \(\mm/\mb = \md - \mc^T \mb^\dagger \mc\) the generalized Schur complement of \(\mb\) in \(\mm\), and note that it is positive semidefinite since \(\mm\) is positive semidefinite.
With the notation in Definition 2, we define the function \(f(\emptyset)=0\) and for \(|\mathcal{S}| = k\), \[f(\ma[\mathcal{S}]) \equiv \sum_{j=1}^k \ma[\mathcal{S}_j]/\ma[\mathcal{S}_{j-1}] = \sum_{j=1}^k\frac{\vol(\ma[\mathcal{S}_j])}{\vol(\ma[\mathcal{S}_{j-1}])},\] where \(\ma[\mathcal{S}_j]/\ma[\mathcal{S}_{j-1}] \in [0,\infty)\) are generalized Schur complements, \(1 \le j \le k\), and we interpret \(0/0 = 0\). This definition is analogous to Definition 5 for general matrices. If \(\ma[\mathcal{S}] = \ml\ml^T\) is a Cholesky factorization, then \(f(\ma[\mathcal{S}]) = \trace(\ml^2) = \sum_{j=1}^k(\ml)_{jj}^2\).
The function \(f\) is non-decreasing, because \(\ma\) is positive semi-definite.
Theorem 14. Let \(\ma\in\rnn\) be symmetric positive semi-definite, and let \(k\) be a fixed integer with \(0<k\leq \rank(\ma)\). Let \(\ma[\mathcal{S}_k^*]\in\real^{k\times k}\) be a principal submatrix of \(\ma\) that maximizes \(f\) over all principal submatrices \(\mb\in\real^{k\times k}\) of \(\ma\).
Then Algorithm 5 is a greedy algorithm in the sense of Algorithm 1 for maximizing \(f\) and its output matrix \(\ma[\mathcal{S}_k]\in\real^{k\times k}\) satisfies \[0 \le f(\ma[\mathcal{S}_k^*]) - f(\ma[\mathcal{S}_k]) \le( 1- \tfrac{1}{k})^k \> f(\ma[\mathcal{S}_k^*]),\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e}< .37\).
The bound in Theorem 14 can also be expressed in terms of Cholesky factors.
Corollary 4. Let \(\ma\in\rnn\) be symmetric positive semi-definite, and let \(0<k\leq \rank(\ma)\) be a fixed integer. Let \(\ma_k^*\in\real^{k\times k}\) be a principal submatrix of \(\ma\) that maximizes \(\log\vol(\mb)\) among all principal submatrices \(\mb\in\real^{k\times k}\) of \(\ma\), and let \(\ma_k\equiv \ma[\mathcal{S}_k]\) be the output matrix of Algorithm 5.
The Cholesky factors from \(\ma_k^*=\ml_k^*(\ml_k^*)^T\) and \(\ma_k=\ml_k\ml_k^T\) satisfy \[\begin{align} 0\leq \trace((\ml_k^*)^2)-\trace((\ml_k)^2)\leq (1-\tfrac{1}{k})^k\,\trace((\ml_k^*)^2), \end{align}\] where \((1-\tfrac{1}{k})^k\leq \tfrac{1}{e} <.37\).
6.4 Interchange algorithm for positive definite matrices↩︎
After relating Algorithm 6 to Algorithm 4 (Theorem 15), we show that Algorithm 6 is a 1-interchange algorithm for maximizing the log-volume of a symmetric positive definite matrix in two ways: in the sense of Algorithm 2 if the relaxation factor \(\eta=1\) (Corollary 5); and in the sense of Theorem 3 if the relaxation factor \(\eta>1\) (Theorem 16).
A conceptual version of [13] is presented in Algorithm 6, with a focus on selecting a \(k\times k\) principal submatrix rather than computing a Cholesky factorization.
Theorem 15 below shows that Algorithm 4 applied to a Cholesky factor \(\mx\) of \(\ma\) implements Algorithm 6. The Cholesky factor \(\mx\) does not have to be a triangular matrix.
Theorem 15. Let \(\ma=\mx^T\mx\) where \(\mx\in\rnn\), and let \(k\) be a fixed integer with \(0<k\leq \rank(\mx)\). Let \(\mx_{\mathcal{S}_k}\in\real^{n\times k}\) be the output matrix of Algorithm 4. Then the output matrix of Algorithm 6 satisfies \[\ma[\mathcal{S}_k]=\mx_{\mathcal{S}_k}^T\mx_{\mathcal{S}_k}.\]
Proof. This is analogous to the proof of Theorem 12. ◻
![Figure 6: Conceptual version of strong rank-revealing Cholesky [13]](/article/2607.13823/_pdflatex/vgcnlomj.png "fig:")
Corollary 5 below shows that Algorithm 6 with relaxation factor \(\eta=1\) is a 1-interchange algorithm for maximizing the log-volume.
Corollary 5. Let \(\ma\in\rnn\) be symmetric positive definite, and let \(k\) be a fixed integer with \(0<k<n\). Let \(\ma_k^*\in\real^{m\times k}\) be a principal submatrix of \(\ma\) that maximizes \(\log{\vol(\mb)}\) among all principal submatrices \(\mb\in\real^{k\times k}\) of \(\ma\).
Then Algorithm 6 with relaxation factor \(\eta=1\) is a 1-interchange algorithm in the sense of Algorithm 2 for maximizing the log-volume, and its output matrix \(\ma_k\equiv\mx_{\mathcal{S}}\) satisfies \[0\leq \log{\vol(\ma_k^*)}-\log{\vol(\ma_k)}\leq \frac{k-1}{2k-1}\;\log{\vol(\ma_k^*/\lambda_n(\mx))},\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\).
Theorem 16 belowshows that Algorithm 6 with relaxation factor \(\eta > 1\) is a 1-interchange algorithm for maximizing the log-volume, and implies that the larger \(\eta\), the larger the error in the log-volume.
Theorem 16. Let \(\ma\in\rnn\) be symmetric positive definite, and let \(0<k<n\) be a fixed integer. Let \(\ma_k^*\in\real^{k\times k}\) be a principal submatrix of \(\ma\) that maximizes \(\log\vol(\mb)\) among all principal submatrices \(\mb\in\real^{k\times k}\) of \(\ma\).
Then Algorithm 6 with relaxation factor \(\eta>1\) is a 1-interchange algorithm in the sense of Theorem 3 for maximizing the log-volume, and its output matrix \(\ma_k\equiv\ma[\mathcal{S}]\) satisfies \[0 \le \log\vol(\ma_k^*)-\log\vol(\ma_k) \le \frac{k}{2k-1}\log\vol(\ma_k^*/\lambda_n(\ma)) + \frac{k^2}{2k-1}\log \sqrt{\eta},\] where \(\frac{k-1}{2k-1}<\frac{1}{2}\) and \(\frac{k^2}{2k-1}\leq k\).
7 Possible directions for future work↩︎
We have established connections between column subset selection and submodular functions, and shown that log-volume is a submodular function on the columns of a matrix.
Avenues for future work include connections to submodular optimization for algorithm development. Since log-volume is in general not monotone, one could look towards greedy algorithms for non-monotone functions in [20], [21]. There are other types of algorithms for submodular functions such as lazy greedy, and stochastic greedy [46]. The distributed greedy method in [3] appears to be related to the notion of tournament pivoting in rank-revealing QR factorizations [47]. There may also be opportunities for randomized QR factorizations such as [48] and [49]. Another option is to relax submodularity and and replace it by approximate submodularity [2] instead. Extensions of the submodular bounds to other algorithms such as local maximum volume submatrix [37] and criteria such as S-opt [50], [51] are also of interest.
References↩︎
Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA. ,↩︎
The work was funded, in part, by the National Science Foundation through awards DMS-2411198 and CCF-2209510, and by the Department of Energy through awards DE-SC0023188 an DE-SC0026310.↩︎
For QR factorizations that are expressed in a way to guarantee non-negative diagonal elements.↩︎