A DEIM Induced CUR Factorization
We derive a CUR matrix factorization based on the Discrete Empirical Interpolation Method (DEIM). For a given matrix $A$, such a factorization provides a low rank approximate decomposition of the form $A \approx C U R$, where $C$ and $R$ are subsets of the columns and rows of $A$, and $U$ is constructed to make $CUR$ a good approximation. Given a low-rank singular value decomposition $A \approx V S W^T$, the DEIM procedure uses $V$ and $W$ to select the columns and rows of $A$ that form $C$ and $R$. Through an error analysis applicable to a general class of CUR factorizations, we show that the accuracy tracks the optimal approximation error within a factor that depends on the conditioning of submatrices of $V$ and $W$. For large-scale problems, $V$ and $W$ can be approximated using an incremental QR algorithm that makes one pass through $A$. Numerical examples illustrate the favorable performance of the DEIM-CUR method, compared to CUR approximations based on leverage scores.

\section{Conclusions}  

\vspace*{-5pt}

\noindent
The Discrete Empirical Interpolation Method (DEIM) is an index selection procedure that gives simple, deterministic CUR factorizations of the matrix ${\bf A}$.\ \ 
Since DEIM utilizes (approximate) singular vectors, we propose an effective one-pass incremental approximate QR factorization that can efficiently compute dominant singular vectors for data sets with rapidly decaying singular values; this method could prove useful in a variety of other settings.    
The accuracy of the resulting rank-$k$ CUR factorization can be bounded in terms of $\sigma_{k+1}$, the error in the best rank-$k$ approximation to ${\bf A}$.\ \  Our analysis of the 2-norm error $\|{\bf A}-{\bf C}{\bf U}{\bf R}\|$ applies to all CUR approximations that use the optimal central factor ${\bf U}={\bf C}^I{\bf A}{\bf R}^I$, and hence can give insight into the performance of other index selection algorithms, such as leverage scores, uniform random sampling, or column-pivoted QR factorization.  Numerical examples illustrate that the DEIM-CUR approach can deliver very good low rank approximations, compared to row selection based on dominant leverage scores.

\section*{Acknowledgements}

\vspace*{-7pt}
We thank Inderjit Dhillon, Petros Drineas, Ilse Ipsen, Michael Mahoney, 
and Nick Trefethen for a number of helpful discussions.  
We are also grateful to Gunnar Martinsson for recommending experiments with the column-pivoted QR selection algorithm, and to an anonymous referee for encouraging us to seek the improved analysis and growth example for the DEIM error constants at the end of Section~\ref{sec:theory}.





   



\bibliographystyle{plain}

\section{DEIM} \label{sec:deim}
The DEIM point selection algorithm was first presented in~\cite{CS10} in the
context of model order reduction for nonlinear dynamical systems, and is a discrete variant of the Empirical Interpolation Method originally proposed in~\cite{BMNP}.   
The DEIM procedure operates on the singular vector matrices ${\bf V}$ and ${\bf W}$ independently to select the row indices ${\bf p}$ and column indices ${\bf q}$.  We explain the process for selecting ${\bf p}$; applying the same steps to ${\bf W}$ yields~${\bf q}$.
To derive the method, we elaborate upon the interpolatory projectors introduced in the last section.

\begin{definition}
Given a full rank matrix ${\bf V} \in \mbox {\Bb R}^{m\times k}$ and a set of distinct indices ${\bf p}\in\mbox {\Bb N}^k$, the \emph{interpolatory projector} for ${\bf p}$ onto ${\rm Ran}({\bf V})$ is
\begin{equation} \label{eq:intproj}
\mbox{\boldmath$\cal P$} \equiv {\bf V}({\bf P}^T{\bf V})^{-1}{\bf P}^T,
\end{equation}
where ${\bf P} = {\bf I}(@:@,{\bf p})\in\mbox {\Bb R}^{m\times k}$, 
provided ${\bf P}^T{\bf V}$ is invertible.
\end{definition}

In general $\mbox{\boldmath$\cal P$}$ is an oblique projector, and it has an important property not generally enjoyed by orthogonal projectors:  for any  ${\bf x}\in\mbox {\Bb R}^m$, 
\[ (\mbox{\boldmath$\cal P$}{\bf x})({\bf p}) \,=\, {\bf P}^T\mbox{\boldmath$\cal P$}{\bf x} \,=\, {\bf P}^T{\bf V}({\bf P}^T{\bf V})^{-1}{\bf P}^T{\bf x} \,=\, {\bf P}^T{\bf x} \,=\, {\bf x}({\bf p}),\]
so the projected vector $\mbox{\boldmath$\cal P$}{\bf x}$ matches ${\bf x}$ in the ${\bf p}$ entries, justifying the name ``interpolatory projector.''


The DEIM algorithm processes the columns of 
\[ {\bf V} = \left[\begin{array}{cccc}{\bf v}_1 & {\bf v}_2 & \cdots & {\bf v}_k\end{array}\right] \]
one at a time, starting from the leading singular vector ${\bf v}_1$.  Each step processes the next singular vector to produce the next index.  The first index $p_1$ corresponds to the largest magnitude entry in ${\bf v}_1$:
\[ |{\bf v}_1(p_1)| = \|{\bf v}_1\|_\infty.\]
Now define ${\bf p}_1 \equiv [p_1]$, and let
\[ \mbox{\boldmath$\cal P$}_1 \equiv {\bf v}_1 ({\bf P}_1^T{\bf v}_1)^{-1}{\bf P}_1^T\]
denote the interpolatory projector for ${\bf p}_1$ onto ${\rm Ran}({\bf v}_1)$.
The second index $p_2$ corresponds to the largest entry in ${\bf v}_2$, after the interpolatory projection in the ${\bf v}_1$ direction has been removed:
\begin{eqnarray*}
 {\bf r}_2\!\! &\equiv&\!\! {\bf v}_2 - \mbox{\boldmath$\cal P$}_1 {\bf v}_2\\[.25em]
 |{\bf r}_2(p_2)|\!\! &=&\!\! \|{\bf r}_2\|_\infty.
\end{eqnarray*}
Notice that ${\bf r}_2(p_1) = 0$, since $\mbox{\boldmath$\cal P$}_1{\bf v}_2$ matches ${\bf v}_2$ in the $p_1$ position, a consequence of interpolatory projection.  This property ensures the process will never produce duplicate indices.

Now suppose we have $j-1$ indices, with
\[ {\bf p}_{j-1} \equiv \left[\begin{array}{c} p_1 \\ \vdots \\ p_{j-1}\end{array}\right], \quad
    {\bf P}_{j-1} \equiv {\bf I}(@:@,{\bf p}_{j-1}), \quad
     {\bf V}_{j-1} \equiv [\begin{array}{ccc} {\bf v}_1 & \cdots & {\bf v}_{j-1}\end{array}], \quad
    \mbox{\boldmath$\cal P$}_{j-1} \equiv {\bf V}_{j-1}^{}({\bf P}_{j-1}^T{\bf V}_{j-1}^{})^{-1}{\bf P}_{j-1}^T.
    \]
To select $p_{j}$, remove from ${\bf v}_j$ its interpolatory projection onto indices ${\bf p}_{j-1}$ and take the largest remaining entry:
\begin{eqnarray*}
 {\bf r}_{j} \!\! &\equiv&\!\! {\bf v}_j - \mbox{\boldmath$\cal P$}_{j-1} {\bf v}_j\\[.25em]
 |{\bf r}_{j}(p_{j})|\!\! &=&\!\! \|{\bf r}_{j}\|_\infty.
\end{eqnarray*}
Implementations should not explicitly construct these projectors; see the pseudocode in Algorithm~\ref{fig:Deim_Arnoldi} for details.

Those familiar with partially pivoted LU decomposition will notice,
on a moment's reflection, 
that this index selection scheme 
is exactly equivalent to the index selection of partial pivoting.   
This arrangement is equivalent to the ``left looking" variant of
LU factorization~\cite[sect.~5.4]{DDSV98}, but with two important differences.  
First, there are no explicit row interchanges in DEIM, as there are
in LU factorization.   Second, the original basis vectors 
(columns of ${\bf V}$) are not replaced
with the residual vectors, as happens in traditional LU decomposition.  
(In the context of model reduction,
it is preferable to keep the nice orthogonal basis intact for use 
as a reduced basis.)
We will exploit this connection with partially pivoted LU factorization 
to analyze the approximation properties of DEIM.

Since the DEIM algorithm processes the singular
vectors sequentially, from most to least significant, it 
 introduces new singular vector information in a 
coherent manner as it successively selects the $k$~indices.  Contrast 
this to index selection strategies based on 
leverage scores, where all singular vectors are 
incorporated at once via row norms of ${\bf V}$ and ${\bf W}$; to account for the fact that higher singular vectors are less significant, such approaches often instead compute leverage scores using only a few of the leading singular vectors.%
\footnote{A potential limitation of the DEIM approach is that ${\bf r}_j$ could have multiple entries that have nearly the same magnitude, but only one index is selected at the $j$th step; if the other large-magnitude entries in ${\bf r}_j$ are not significant in subsequent ${\bf r}_\ell$ vectors, the corresponding indices will not be selected.  One can imagine modifications of the selection algorithm to account for such situations, e.g., by processing multiple singular vectors at a time.}


\begin{algorithm}[t!]
\begin{center}
\fbox{
\begin{tabular}{p{10cm}}
\ \\[-25pt]
\begin{tabbing}
defghasijkl\=bbb\=ccc\= \kill
$\ \ \ $ {\bf Input:} \> $ {\bf V}$,  an $m \times k$ matrix ($m\ge k$)\\[5pt]
$\  $  {\bf Output:}  \>  ${\bf p}$, an integer vector with $k$ distinct entries in $ \{1,\ldots,m \}$ \\[7pt]
$\ \ \ $          \> $ {\bf v} = {\bf V}(:,1)$  \\
$\ \ \ $          \> $[\sim , p_1] = {\rm max}( | {\bf v} |)$\\
$\ \ \ $           \>$ {\bf p} = [p_1] $ \\
$\ \ \ $           \>{\bf for} $j=2,3,\ldots, k$  \\
$\ \ \ $               \>\> $ {\bf v} = {\bf V}(:,j)$  \\
$\ \ \ $               \>\> $ {\bf c} =  {\bf V}({\bf p},1:j-1)^{-1} {\bf v}({\bf p})$  \\
$\ \ \ $               \>\> $ {\bf r} = {\bf v} - {\bf V}(:,1:j-1) {\bf c}$  \\
$\ \ \ $               \>\>$[\sim , p_j] = {\rm max}(| {\bf r} |)$\\ 
$\ \ \ $               \>\>$ {\bf p} = [{\bf p};\ p_j] $ \\
$\ \ \ $             \>{\bf end}\\
\end{tabbing}\\[-20pt]\end{tabular}}
\end{center}

\vspace*{-7pt}
\caption {DEIM point selection algorithm.}
\label{fig:Deim_Arnoldi}
\end{algorithm}


For the interpolatory projector $\mbox{\boldmath$\cal P$}_j$ to exist at the $j$th step, ${\bf P}_{j-1}^T{\bf V}_{j-1}$ must be nonsingular.  The linear independence of the columns of ${\bf V}$ assures this.  In the following, ${\bf e}_j$ denotes the $j$th column of the identity matrix.

\begin{lemma} \label{DEIM_inv}
Let ${\bf P}_j = [ {\bf e}_{p_1}, {\bf e}_{p_2}, \ldots, {\bf e}_{p_j}]$ and
let ${\bf V}_j = [ {\bf v}_1, {\bf v}_2, \ldots, {\bf v}_j]$ for $1 \le j \le k$.
If $ { {{\rm rank }\,} } ({\bf V}) = k$, then ${\bf P}_j^T {\bf V}_j$ is nonsingular for $1 \le j \le k$.
\end{lemma}
\noindent{\it Proof:} 
Suppose ${\bf P}_{j-1}^T {\bf V}_{j-1} $ is nonsingular and
let ${\bf r}_j = {\bf v}_j - {\bf V}_{j-1}({\bf P}_{j-1}^T {\bf V}_{j-1})^{-1} {\bf P}_{j-1}^T {\bf v}_j.$
Then $\|{\bf r}_j\|_{\infty}  > 0$, for otherwise
$ {\bf 0} = {\bf v}_j - {\bf V}_{j-1} {\bf c}_{j-1}$, in violation of the assumption that $ { {{\rm rank }\,} }  ({\bf V}) = k$.
Thus
\begin{equation}
\label{eta_formula}
   0 < |{\bf e}_{p_j}^T {\bf r}_j| =| {\bf e}_{p_j}^T {\bf v}_j - {\bf e}_{p_j}^T {\bf V}_{j-1}({\bf P}_{j-1}^T {\bf V}_{j-1})^{-1} {\bf P}_{j-1}^T {\bf v}_j |,
\end{equation}
where $p_j$ is the $j$th DEIM interpolation point.
Now factor
\begin{equation}
\label{PV_formula}
   {\bf P}_j^T {\bf V}_j =  
                   \left[ \begin{array}{cc} {\bf P}_{j-1}^T {\bf V}_{j-1} & {\bf P}_{j-1}^T{\bf v}_j \\
                                             {\bf e}_{p_j}^T {\bf V}_{j-1} & {\bf e}_{p_j}^T {\bf v}_j  
                            \end{array}
                     \right]=
                   \left[ \begin{array}{cc} {\bf I}_{j-1} & {\bf 0} \\
                                             {\bf e}_{p_j}^T {\bf V}_{j-1}({\bf P}_{j-1}^T {\bf V}_{j-1})^{-1} & 1  
                            \end{array}
                     \right]
                   \left[ \begin{array}{cc} {\bf P}_{j-1}^T {\bf V}_{j-1} & {\bf P}_{j-1}^T{\bf v}_j \\
                                             {\bf 0} & \nu_j  
                            \end{array}
                     \right],
\end{equation}
where 
\[
   \nu_j = {\bf e}_{p_j}^T {\bf v}_j - {\bf e}_{p_j}^T {\bf V}_{j-1}({\bf P}_{j-1}^T {\bf V}_{j-1})^{-1} {\bf P}_{j-1}^T {\bf v}_j .
\]
The inequality~(\ref{eta_formula}) implies $\nu_j \ne 0$ and hence equation~(\ref{PV_formula})
implies ${\bf P}_j^T {\bf V}_j$ is nonsingular.  Since ${\bf e}_{p_1}^T{\bf v}_1 \ne 0$, this argument
provides an inductive proof that  ${\bf P}_j^T {\bf V}_j $  is nonsingular for  $1 \le j \le k$.   { { \hfill\rule{3mm}{3mm} } } 


\section{Computational Examples} \label{sec:examples}

This section presents some computational evidence illustrating the excellent approximation
properties of the DEIM-CUR factorization, consistent with the error analysis in Section~\ref{sec:theory}.  
For each of our three examples, we compare the accuracy of the DEIM-CUR factorization with several schemes based on leverage scores.  To remove random variations from our experiments, in most cases we  select columns and rows having the highest leverage scores; for the first example, we include results for random leverage score sampling.
For Example~1 we also study the effect of inaccurate singular vectors on the DEIM selection, and compare the accuracy of DEIM-CUR to CUR approximations based on the column-pivoted QR algorithm.

\newpage

\noindent \textbf{Example 1}.  \emph{Low-rank approximation of a sparse, nonnegative matrix}

\medskip
\noindent
The first example builds a matrix ${\bf A}\in\mbox {\Bb R}^{300,000\times 300}$ of the form
\begin{equation} \label{eq:sparsex2}
{\bf A} = \sum_{j=1}^{10} {2\over j}\,{\bf x}_j^{}{\bf y}_j^T 
           + \sum_{j=11}^{300} {1\over j}\,{\bf x}_j^{}{\bf y}_j^T,
\end{equation}
where ${\bf x}_j \in \mbox {\Bb R}^{300,000}$ and ${\bf y}_j \in \mbox {\Bb R}^{300}$ are sparse vectors with random nonnegative entries (in MATLAB, ${\bf x}_j = {\tt sprand(300000,1,0.025)}$ and ${\bf y}_j = {\tt sprand(300,1,0.025)}$).  In this instantiation, ${\bf A}$ has 15,971,584 nonzeros, i.e., about 18\% of all entries are nonzero.  The form~(\ref{eq:sparsex2}) is not a singular value decomposition, since $\{{\bf x}_j\}$ and $\{{\bf y}_j\}$ are not orthonormal sets; however, this decomposition suggests the structure of the SVD:  the singular values decay like $1/j$, and with the first ten singular values weighted more heavily to give a notable drop between $\sigma_{10}$ and $\sigma_{11}$.
We begin these experiments by computing ${\bf V}$ and ${\bf W}$ using 
MATLAB's economy-sized SVD routine (\verb|[V,S,W] = svd(A,'0')|).


\begin{figure}[b!]
\begin{center}
\includegraphics[scale=0.65]{dirty_sparse4_err_q0}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:sparse_ex4_cur}
Accuracy of CUR approximations for the sparse, nonnegative matrix~(\ref{eq:sparsex2}) using $k$ columns and rows, constructed by DEIM-CUR and two leverage score strategies:  ``LS (all)'' selects rows and columns with highest leverage scores computed using all 300~singular vectors; ``LS (10)'' only uses the leading ten singular vectors.  The ``DEIM($\widehat{{\bf V}},\widehat{{\bf W}}$)'' curve (nearly atop the ``DEIM'' curve) uses approximate singular vectors, described later.}
\end{figure}

Figure~\ref{fig:sparse_ex4_cur} compares the error $\|{\bf A}-{\bf C}{\bf U}{\bf R}\|$ for DEIM-CUR and methods that take ${\bf C}$ and ${\bf R}$ as the columns and rows of ${\bf A}$ with the highest leverage scores.  These scores are computed using either all right and left singular vectors (300 of each), or using only the leading ten right and left singular vectors.  Both  approaches perform rather worse than DEIM-CUR, which closely tracks the optimal value $\sigma_{k+1}$.  

To gain insight into these results, we examine the interpolation constants $\eta_p$ and $\eta_q$ for all three approaches.  
Figure~\ref{fig:sparse_ex4_eta} shows that these  constants are largest for leverage scores based on all the singular vectors; using only ten singular vectors improves both the interpolation constants and the accuracy of the approximation (as seen in Figure~\ref{fig:sparse_ex4_cur}).  
The DEIM-CUR method gives better interpolation constants and more accurate approximations.  

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{sparse_ex4_eta_ls13}\quad
\includegraphics[scale=0.45]{sparse_ex4_eta_d1}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:sparse_ex4_eta}
Error constants $\eta_p = \|({\bf P}_k^T{\bf V}_k)^{-1}\|$ and $\eta_q = \|({\bf W}_k^T{\bf Q})^{-1}\|$ for rows and columns selected using
 two leverage score strategies (left plot) and the DEIM algorithm (right plot), for the matrix~(\ref{eq:sparsex2}). }
\end{figure}
 
A CUR factorization can also be obtained by randomly sampling columns and rows of ${\bf A}$, with the probability of selection weighted by leverage scores~\cite{MD09}.  We apply this approach on the current example, selecting $k=30$ rows and columns of ${\bf A}$ with a probability given by the leverage scores computed from the leading ten singular vectors (normalized to give a probability distribution).  Figure~\ref{fig:sparse_ex5_cur} gives the results of ten independent experiments, showing that while sampling can sometimes yield better results than the deterministic leverage score approach, overall the approximations are still inferior to those from DEIM-CUR.

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.65]{sparse_ex5_cur}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:sparse_ex5_cur}
Accuracy of CUR approximations for~(\ref{eq:sparsex2}) generated by randomly sampling rows and columns with probability weighted by leverage scores computed from the leading ten singular vectors.  All ten~trials (gray lines) perform similarly to the deterministic ``LS (10)'' approach, and worse than the DEIM-CUR approximation.}
\end{figure}

How robust is the DEIM-CUR approximation to errors in the singular vectors?
To investigate, we compute 
$\widehat{{\bf V}}\approx {\bf V}$ and $\widehat{{\bf W}}\approx {\bf W}$ using the
Incremental QR algorithm detailed in Section~\ref{sec:qr} (with ${\it tol}=10^{-4}$)
and the Randomized SVD algorithm described by
Halko, Martinsson, and Tropp~\cite[p.~227]{HMT11}.  
To give extreme examples of the latter, we compute $\widehat{{\bf V}}$ and $\widehat{{\bf W}}$ through 
\emph{only one or two} applications each of ${\bf A}$ and ${\bf A}^T$.%
\footnote{
This corresponds to $q=0$ and $q=1$ in the notation of~\cite[p.~227]{HMT11}.
Let the columns of ${\bf Q}\in \mbox {\Bb R}^{m\times 2 k_{\rm max}}$ form an orthonormal basis for 
$({\bf A}\bA^T)^q{\bf A}\mbox{\boldmath$\Omega$}$, where $\mbox{\boldmath$\Omega$}\in \mbox {\Bb R}^{n\times 2 k_{\rm max}}$ is a 
random matrix with i.i.d.\ Gaussian entries and we take $k_{\rm max} = 30$.
Then the leading $k_{\rm max}$~columns of ${\bf V}$ and ${\bf W}$ are approximated 
by taking the SVD of ${\bf Q}^*{\bf A} \in \mbox {\Bb R}^{2k_{\rm max} \times n}$. 
}
As Figure~\ref{fig:dirty_sparse4_angle} illustrates, in both cases the angle between the exact and approximate leading singular subspaces is significant, particularly as $k$ grows.
This drift in the subspaces has little effect on the accuracy of the DEIM approximations.

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{dirty_sparse4_angle_op}\quad
\includegraphics[scale=0.45]{dirty_angles}

\begin{picture}(0,0)
\put(-202,40){\small \emph{SVD via Incremental QR, ${\it tol}=10^{-4}$}}
\put(38,145){\rotatebox{10}{\small \emph{one application of ${\bf A}$ and ${\bf A}^{\kern-1pt T}$}}}
\put(47,60){\rotatebox{33}{\small \emph{two applications of ${\bf A}$ and ${\bf A}^{\kern-1pt T}$}}}
\put(45,40){\small \emph{Randomized SVD}}
\end{picture}
\end{center}

\vspace*{-32pt}
\caption{\label{fig:dirty_sparse4_angle}
The angle between the leading $k$-dimensional exact singular subspaces 
${\rm Ran}({\bf V}_k)$ and ${\rm Ran}({\bf W}_k)$
(generated by MATLAB's {\tt svd} command) and approximate
singular subspaces ${\rm Ran}(\widehat{{\bf V}}_k)$ and ${\rm Ran}(\widehat{{\bf W}}_k)$
for the matrix~(\ref{eq:sparsex2}).
On the left, $\widehat{{\bf V}}_k$ and $\widehat{{\bf W}}_k$ are generated using the Incremental QR algorithm described in Section~\ref{sec:qr}, with ${\it tol} = 10^{-4}$;
on the right, $\widehat{{\bf V}}_k$ and $\widehat{{\bf W}}_k$ are
generated using randomized SVD algorithm~\cite{HMT11}  using one and two applications of ${\bf A}$ and ${\bf A}^T$.}
\end{figure}

\begin{itemize}
\item The DEIM approximation using the Incremental QR algorithm is quite robust, choosing at most 3~different row indices and 2~different column indices for $k=1,\ldots,30$, with a relative discrepancy in $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$ of at most 9.27\% (and this realized only at step $k=30$).
\item When ${\bf A}$ and ${\bf A}^T$ are applied once in the Randomized SVD algorithm, the DEIM indices differ considerably from those drawn from exact singular vectors (e.g., for $k=30$, 20 of  30 row indices and 3 of 30 column indices differ), yet the quality of the approximation $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$ remains almost the same (relative difference of at most 10.45\%); see the dashed line in Figure~\ref{fig:sparse_ex4_cur}.
\item When ${\bf A}$ and ${\bf A}^T$ are applied twice, the DEIM indices are nearly identical (e.g., for $k=30$, 0 of 30 row indices and 2 of 30 column indices differ). On the scale of the plot in Figure~\ref{fig:sparse_ex4_cur}, $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$ could not be distinguished from the DEIM-CUR errors using exact singular vectors; the maximum relative discrepancy is 2.21\%.
\end{itemize}

\begin{figure}[b!]
\begin{center}
\includegraphics[scale=0.52]{test10_deimqrcomp}
\begin{picture}(0,0)
\put(-200,-13){\small ratio (DEIM-CUR error)/(QR-CUR error)}
\put(-215,175){\small DEIM-CUR better} 
\put(-103,175){\small QR-CUR better} 
\put(0,82){\small \rotatebox{90}{rank, $k$}}
\end{picture}

\vspace*{2em}
\includegraphics[scale=0.45]{test10_deim_p}\qquad
\includegraphics[scale=0.45]{test10_qr_p}

\begin{picture}(0,0)
\put(-163,0){\small $\log_{10}(\eta_p)$ for DEIM-CUR}
\put( 82,0){\small $\log_{10}(\eta_p)$ for QR-CUR}
\put(-243,82){\small \rotatebox{90}{rank, $k$}}
\put(0,82){\small \rotatebox{90}{rank, $k$}}
\end{picture}
\end{center}

\vspace*{-10pt}
\caption{\label{fig:test10}
Comparison of DEIM-CUR and QR-CUR performance for
100~sparse random $300,000\times 300$ matrices of the form~(\ref{eq:sparsex2}).  
The top plot shows a histogram of the ratio of $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$ for DEIM-CUR and QR-CUR.  
The bottom plots compare the error constant $\eta_p = \|({\bf P}_k^T{\bf V}_k)^{-1}\|$ for DEIM-CUR (left) and QR-CUR (right); 
note the logarithmic scale of the horizontal axes in the lower plots.} 
\end{figure}

Thus far we have only compared the DEIM-CUR approximations to CUR factorizations obtained from leverage scores, which also use singular vector information, thus illustrating how DEIM can use the same raw materials to better effect.
Next we compare DEIM-CUR to approximations computed using a different approach based on QR factorization of ${\bf A}$; see, e.g., \cite{CGMR05,Ste99}.
Begin by computing a column-pivoted QR factorization of ${\bf A}$; the first $k$ selected columns  give the indices ${\bf q}$, from which we extract ${\bf C}_k = {\bf A}(@@:@@,{\bf q})$.
Next, a column-pivoted QR factorization of ${\bf C}_k^T$ is performed; the first $k$ selected columns of ${\bf C}_k^T$ 
give the indices ${\bf p}$, from which we build ${\bf R}_k = {\bf A}({\bf p},@@:@@)$.
We refer to this technique as ``QR-CUR.''

Figure~\ref{fig:test10} compares the results for 100~trials 
involving sparse random matrices of dimension $300,000\times 300$ 
having the form of our first experiment~(\ref{eq:sparsex2}).  DEIM-CUR and QR-CUR 
produce factorizations with similar accuracy, which we illustrate with a histogram of the ratio of $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$ for DEIM-CUR to QR-CUR, for $k=1,\ldots, 100$.  
(DEIM-CUR produces a smaller error when the ratio is less than one.)
While these errors are similar, the error constants $\eta_p$ and $\eta_q$ 
for the two methods are quite different.  
The bottom plots in Figure~\ref{fig:test10} compare histograms 
of $\log_{10} \eta_p$.  
For DEIM-CUR, the $\eta_p$ values are quite consistent
across the 100~random ${\bf A}$, while for QR-CUR the $\eta_p$ values
are both larger and rather less consistent.  (The figures for $\eta_q$ are qualitatively identical, but about an order of magnitude smaller for both methods.)  

\medskip
The advantage of DEIM-CUR over approximations based on leverage scores remains when the singular values decrease more sharply.
Modify~(\ref{eq:sparsex2}) to give a more significant drop between $\sigma_{10}$ and $\sigma_{11}$:
\begin{equation} \label{eq:sparsex1000}
{\bf A} = \sum_{j=1}^{10} {1000\over j}\,{\bf x}_j^{}{\bf y}_j^T 
           + \sum_{j=11}^{300} {1\over j}\,{\bf x}_j^{}{\bf y}_j^T.
\end{equation}
As seen in Figures~\ref{fig:sparse_ex3_cur} and~\ref{fig:sparse_ex3_eta}, the DEIM-CUR approach again delivers excellent approximations, while selecting the rows and columns with highest leverage scores does not perform nearly as well.  (In Figure~\ref{fig:sparse_ex3_eta}, note the significant jump in the ``LS~(10)'' error constant $\eta_q$ corresponding to those $k$ values where $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|/\sigma_{k+1}$ is large.)
 
\begin{figure}[ht!]
\begin{center}
\includegraphics[scale=0.65]{sparse_ex3_cur}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:sparse_ex3_cur}
Accuracy of CUR approximations using $k$ rows and columns, for DEIM-CUR and two leverage score strategies for the sparse, nonnegative matrix~(\ref{eq:sparsex1000}).  ``LS (all)'' selects rows and columns having the highest leverage scores computed using all 300~singular vectors; ``LS (10)'' uses the leading 10~singular vectors.}
\end{figure}

\begin{figure}[ht!]
\begin{center}
\includegraphics[scale=0.45]{sparse_ex3_eta_ls13}\quad
\includegraphics[scale=0.45]{sparse_ex3_eta_d1}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:sparse_ex3_eta}
Error constants $\eta_p = \|({\bf P}_k^T{\bf V}_k)^{-1}\|$ and $\eta_q = \|({\bf W}_k^T{\bf Q}_k)^{-1}\|$ for rows and columns selected using two 
leverage score strategies (left plot) and the DEIM algorithm (right plot), for the sparse matrix ${\bf A}$ given in~(\ref{eq:sparsex1000}). }
\end{figure}

\newpage
\noindent \textbf{Example 2}.  \emph{TechTC term document data}

\medskip
\noindent
The second example, adapted from Mahoney and Drineas~\cite{MD09}, computes the CUR factorization of a term document matrix with data drawn from the Technion Repository of Text Categorization Datasets (TechTC)~\cite{GM04}.  The rows of the data matrix correspond to websites (consolidated from multiple webpages), while the columns correspond to ``features'' (words from the text of the webpages).  The $(j,k)$ entry of ${\bf A}$ reflects the importance of the feature text on the given website; most entries are zero. For this experiment we use TechTC-100 test set~26, which concatenates a data set relating to Evansville, Indiana (id~10567) with another for Miami, Florida (id~11346).  Following Mahoney and Drineas~\cite{MD09}, we omit all features with four or fewer characters from the data set, leaving a matrix with 139~rows and 15,170~columns.  Each row of ${\bf A}$ is then scaled to have unit 2-norm.  Ideally a CUR factorization not only gives an accurate low-rank approximation to ${\bf A}$, but also selects rows corresponding to representative webpages from each geographic area, and columns corresponding to meaningful features.  

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.65]{techtc_rs_cur}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:techtc_rs_cur} Accuracy of CUR factorizations for the TechTC example, selecting rows and columns using top leverage scores for all singular vectors and the leading two singular vectors, and DEIM.}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.65]{techtc_rs_comp_ls2}
\end{center}

\vspace*{-20pt}
\caption{\label{fig:techtc_rs_comp_ls2} The columns selected by DEIM for the TechTC example, compared to leverage scores from the leading two singular vectors.}
\vspace*{-10pt}
\end{figure}

Figure~\ref{fig:techtc_rs_cur} compares DEIM-CUR approximations to row and column selection based on highest leverage scores (from all singular vectors, or the two leading singular vectors).  The DEIM-CUR approximations are typically more accurate than those based on leverage scores, but all approaches give errors roughly two times larger than the slowly-decaying optimal value of $\sigma_{k+1}$. 
How do the DEIM columns (features) compare to those with the highest leverage scores?  Figure~\ref{fig:techtc_rs_comp_ls2} shows the leverage scores associated with each column of ${\bf A}$ (based on the two leading singular vectors), along with the first 30~columns selected by DEIM.  While the columns with highest leverage scores were found by DEIM, there are DEIM columns with marginal leverage scores, and vice versa.  This data is more easily parsed in Table~\ref{tbl:techtc}, which lists the features corresponding to the first 20~DEIM columns.  (To ease comparison, we normalize leverage scores so that the maximum value is one.)  The leading features identified by DEIM, including ``evansville'' (first DEIM point), ``florida'' (second), ``miami'' (sixth), and ``indiana'' (nineteenth), indeed reveal key geographic terms.  These terms scored at least as high when ranked by leverage scores based on two leading singular vectors; when all singular vectors are used, the scores of these terms generally drop, relative to other features.  Overall, one notes that DEIM selects a significantly different set of indices than those valued by leverage scores, and, as seen in Figure~\ref{fig:techtc_rs_cur}, tends to provide a somewhat better low-rank approximation.



\begin{table}[t!]
\caption{\label{tbl:techtc}
The features selected by DEIM-CUR for the TechTC data set, compared to the (scaled) leverage scores using the leading two singular vectors, and all singular vectors.}

\vspace*{-5pt}
\begin{center}\small
\begin{tabular}{cr|cc|cc|l} \hline
\multicolumn{1}{c}{DEIM } &
\multicolumn{1}{c}{index} &
\multicolumn{2}{|c}{LS (2)} &
\multicolumn{2}{|c}{LS (all) } &
\multicolumn{1}{|c}{} \\
\multicolumn{1}{c}{rank } &
\multicolumn{1}{c}{$q_j$} &
\multicolumn{1}{|c}{rank } &
\multicolumn{1}{c}{score} &
\multicolumn{1}{|c}{rank } &
\multicolumn{1}{c}{score} &
\multicolumn{1}{|c}{feature} \\ \hline
     1 &      10973 &      1 &  1.000 &      4 &  0.875 & evansville\\
     2 &          1 &      2 &  0.741 &      8 &  0.726 & florida\\
     3 &       1547 &     13 &  0.031 &      2 &  0.948 & spacer\\
     4 &        109 &      8 &  0.055 &     66 &  0.347 & contact\\
     5 &        209 &     12 &  0.040 &     32 &  0.458 & service\\
     6 &         50 &      4 &  0.116 &      6 &  0.739 & miami\\
     7 &        824 &     46 &  0.007 &      5 &  0.809 & chapter\\
     8 &       1841 &     33 &  0.010 &     20 &  0.537 & health\\
     9 &        171 &      5 &  0.113 &     13 &  0.617 & information\\
    10 &        234 &     16 &  0.026 &     37 &  0.436 & events\\
    11 &        595 &     84 &  0.004 &     15 &  0.576 & church\\
    12 &         60 &     15 &  0.026 &     67 &  0.347 & email\\
    13 &        945 &     10 &  0.047 &     30 &  0.474 & services\\
    14 &       1670 &    129 &  0.002 &      1 &  1.000 & bullet\\
    15 &        216 &     35 &  0.009 &     38 &  0.430 & music\\
    16 &         78 &      3 &  0.246 &     24 &  0.492 & south\\
    17 &        213 &     19 &  0.018 &    110 &  0.259 & their\\
    18 &        138 &     14 &  0.030 &     43 &  0.408 & please\\
    19 &       6110 &      7 &  0.060 &     95 &  0.280 & indiana\\
    20 &       1152 &     70 &  0.005 &    152 &  0.221 & member
\end{tabular}\end{center}
\vspace*{-2em}
\end{table}


\medskip
\noindent \textbf{Example 3}.  \emph{Tumor detection in genetics data}

\medskip
\noindent
Our final example uses the GSE10072 cancer genetics data set from the National Institutes of Health, previously investigated by Kundu, Nambirijan, and Drineas~\cite{KND}.  The matrix ${\bf A}\in\mbox {\Bb R}^{22,283\times 107}$ contains data for 22,283~probes applied to 107~patients.  The $(j,k)$ entry of ${\bf A}$ reflects how strongly patient $k$ responded to probe $j$.  This experiment seeks probes that segment the population into two clusters: the 58~patients with tumors, and the 49~without.%
\footnote{The data is available from \url{http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE10072}.}
To center the data, we subtract the mean of each row from all entries in that row.  
As shown in~\cite{KND}, the leading two principal vectors of this matrix segment the population very well. 

Like the TechTC data, the singular values of ${\bf A}$ decay slowly, as seen in Figure~\ref{fig:gene_cur}.  Once again the DEIM-CUR procedure produces a more accurate low-rank approximation than obtained by selecting the rows and columns with highest leverage scores, whether those are computed using all the singular vectors, or just the leading two or ten.


\begin{figure}[ht]
\begin{center}
 \includegraphics[scale=0.65]{gene_cur}
 \end{center}
 
 \vspace*{-17pt}
 \caption{\label{fig:gene_cur}
 Accuracy of CUR factorizations for a genetics data set.
DEIM-CUR consistently outperforms factorizations
 derived by taking the rows and columns with largest leverage scores,
 regardless of whether these scores are drawn from all singular vectors,
 the leading ten singular vectors, or the leading two singular vectors.
 }
 \end{figure}


Table~\ref{tbl:gene_deim} reports the first 15~rows selected by the DEIM-CUR process, along with the corresponding leverage scores based on two, ten, and all singular vectors.  Do the probes selected by DEIM  discriminate the patients with tumors (``sick'') from those without (``well'')? 
To investigate, for each selected probe we count the number of large positive entries corresponding to sick and well patients.%
\footnote{In particular, we call an entry of the mean-centered matrix ${\bf A}$ large if its value exceeds one.  Of the 22,283 probes, for only 23~probes do at least~30 of the~58 sick patients have such large entries; for only~95 probes do at least~30 of the~49 well patients have large entries.  There is no overlap between the probes that are strongly expressed by the sick and well patients.}  Some but not all of the DEIM-CUR probes effectively select only sick or well patients.
Contrast these results with Table~\ref{tbl:gene_ls2}, which shows the probes with highest leverage scores (drawn from the leading two singular vectors).  Only four of these probes were also selected by the DEIM procedure (even if we continue the DEIM procedure to select the maximum number, $n=107$, of indices).  This discrepancy is quite different from the good agreement between DEIM and leverage score indices for the TechTC data in Table~\ref{tbl:techtc},  despite the similar dimensions and the comparably slow decay of the singular values.  


\begin{table}[t!]
\caption{\label{tbl:gene_deim}
Genetics example: the probes selected by DEIM-CUR, compared to the (scaled) leverage scores using the leading two singular vectors, ten singular vectors, and all singular vectors.}
\begin{center}\small
\begin{tabular}{cr|ll|cc|cc|cc|cc} \hline
\multicolumn{1}{c}{DEIM } &
\multicolumn{1}{c}{index} &
\multicolumn{1}{|c}{probe} &
\multicolumn{1}{c}{gene} &
\multicolumn{1}{|c}{number} &
\multicolumn{1}{c}{number} &
\multicolumn{2}{|c}{LS (2) } &
\multicolumn{2}{|c}{LS (10)} &
\multicolumn{2}{|c}{LS (all) } \\
\multicolumn{1}{c}{rank } &
\multicolumn{1}{c}{$q_j$} &
\multicolumn{1}{|c}{set } &
\multicolumn{1}{c}{name} &
\multicolumn{1}{|c}{sick } &
\multicolumn{1}{c}{well} &
\multicolumn{1}{|c}{rank } &
\multicolumn{1}{c}{score} &
\multicolumn{1}{|c}{rank } &
\multicolumn{1}{c}{score} &
\multicolumn{1}{|c}{rank } &
\multicolumn{1}{c}{score}  \\ \hline
     1 &       9565 & 210081\_at &  AGER &   2 &  45 &     1 &  1.000 &    45 &  0.504 &    386 &  0.123 \\
     2 &      14270 & 214895\_s\_at & ADAM10  &   8 &   3 &   211 &  0.173 &  1171 &  0.108 &   3344 &  0.036 \\
     3 &       8650 & 209156\_s\_at & COL6A2  &   5 &   6 & 15156 &  0.005 &   252 &  0.245 &    708 &  0.091 \\
     4 &      11057 & 211653\_x\_at & AKR1C2  &  18 &   1 &  6440 &  0.017 &    11 &  0.656 &    146 &  0.185 \\
     5 &      14153 & 214777\_at & IGKV4-1  &  27 &   3 &   281 &  0.148 &    19 &  0.607 &    106 &  0.209 \\
     6 &      18976 & 219612\_s\_at & FGG  &  17 &  17 &  2591 &  0.039 &     2 &  0.956 &      4 &  0.825 \\
     7 &       3831 & 204304\_s\_at & PROM1  &  16 &   4 &   992 &  0.073 &    70 &  0.417 &     32 &  0.345 \\
     8 &       3351 & 203824\_at & TSPAN8  &  17 &   4 &  9687 &  0.011 &    21 &  0.582 &     31 &  0.355 \\
     9 &       4275 & 204748\_at & PTGS2  &  18 &  14 &   424 &  0.118 &    13 &  0.624 &     42 &  0.313 \\
    10 &       1437 & 201909\_at & RPS4Y1  &  21 &  34 &  8232 &  0.013 &     3 &  0.913 &      5 &  0.736 \\
    11 &      14150 & 214774\_x\_at & TOX3  &  34 &   0 &    95 &  0.262 &    49 &  0.492 &    102 &  0.210 \\
    12 &      10518 & 211074\_at & FOLR1  &   7 &   4 &  9482 &  0.011 &   926 &  0.124 &    213 &  0.159 \\
    13 &       9580 & 210096\_at & CYP4B1  &   6 &  44 &     8 &  0.797 &    65 &  0.431 &     54 &  0.284 \\
    14 &       4002 & 204475\_at & MMP1 &  27 &   0 &    34 &  0.406 &    24 &  0.564 &     21 &  0.465 \\
    15 &      13990 & 214612\_x\_at & MAGEA  &  16 &   0 &   489 &  0.110 &   134 &  0.323 &     35 &  0.339 
    \end{tabular}\end{center}
    \vspace*{-18pt}
\end{table}



\begin{table}[t!]
\caption{\label{tbl:gene_ls2}
Genetics example: the probes with top (scaled) leverage scores, derived from the first two singular vectors.}

\begin{center}\small
\begin{tabular}{crr|ll|cc|c} \hline
\multicolumn{1}{c}{LS (2) } &
\multicolumn{1}{c}{index} &
\multicolumn{1}{c}{LS (2)} &
\multicolumn{1}{|c}{probe} &
\multicolumn{1}{c}{gene} &
\multicolumn{1}{|c}{number} &
\multicolumn{1}{c}{number} &
\multicolumn{1}{|c}{DEIM } \\
\multicolumn{1}{c}{rank } &
\multicolumn{1}{c}{$q_j$} &
\multicolumn{1}{c}{score} &
\multicolumn{1}{|c}{set} &
\multicolumn{1}{c}{name} &
\multicolumn{1}{|c}{sick } &
\multicolumn{1}{c}{well} &
\multicolumn{1}{|c}{rank } \\ \hline
     1 &  9565 &  1.000 & 210081\_at & AGER &   2 &  45 &    1  \\ 
     2 & 13766 &  0.922 & 214387\_x\_at & SFTPC&   6 &  48 &  ---  \\ 
     3 & 11135 &  0.907 & 211735\_x\_at & SFTPC&   5 &  48 &  73  \\ 
     4 &  9361 &  0.899 & 209875\_s\_at & SPP1&  50 &   2 &  ---  \\ 
     5 &  5509 &  0.896 & 205982\_x\_at & SFTPC&   5 &  48 &  ---  \\ 
     6 &  9103 &  0.835 & 209613\_s\_at & ADH1B&   2 &  47 &  ---  \\ 
     7 & 14827 &  0.834 & 215454\_x\_at & SFTPC&   0 &  46 &  ---  \\ 
     8 &  9580 &  0.797 & 210096\_at & CYP4B1&   6 &  44 &   13  \\ 
     9 &  4239 &  0.754 & 204712\_at & WIF1&   5 &  43 &  70  \\ 
    10 &  3507 &  0.724 & 203980\_at & FABP4&   2 &  44 &  ---  \\ 
    11 & 18594 &  0.717 & 219230\_at & TMEM100&   2 &  38 &  ---  \\ 
    12 &  9102 &  0.684 & 209612\_s\_at & ADH1B&   2 &  46 &  ---  \\ 
    13 & 13514 &  0.626 & 214135\_at & CLDN18&   3 &  47 &  ---  \\ 
    14 &  5393 &  0.626 & 205866\_at & FCN3&   0 &  39 &  ---  \\ 
    15 &  4727 &  0.614 & 205200\_at & CLEC3B&   0 &  39 &  ---  \\ 
\end{tabular}\end{center}
    \vspace*{-18pt}

    \end{table}
    


While the rows selected from leverage scores did not produce as accurate an approximation, $\|{\bf A}-{\bf C}_k{\bf U}_k{\bf R}_k\|$, as DEIM, these probes do a much more effective job of discriminating patients with tumors from those without.  Indeed, for~14 of the top~15 probes, the tumor-free patients express strongly, while the patients with tumors do not;
in the remaining case, the opposite occurs.


\section{CUR Factorization} \label{sec:factor}
We are concerned with large matrices ${\bf A}\in\mbox {\Bb R}^{m\times n}$ that represent nearly low-rank data, which can therefore be expressed as
\begin{equation} \label{eq:CURF}
  {\bf A} = {\bf C} {\bf U} {\bf R}  + {\bf F},
\end{equation}
with $\|{\bf F}\|$ small relative to $\|{\bf A}\|$.\ \ 
The matrix ${\bf C}\in\mbox {\Bb R}^{m\times k}$ is formed by extracting $k$ columns from ${\bf A}$, and ${\bf R}\in\mbox {\Bb R}^{k\times n}$ from $k$ rows of ${\bf A}$.  The selected row and column indices are
stored in the vectors ${\bf p}, {\bf q}\in\mbox {\Bb N}^k$, so that ${\bf C} = {\bf A}(:,{\bf q})$ and ${\bf R} = {\bf A}({\bf p},:)$.  
Our choice for ${\bf p}$ and ${\bf q}$ is guided by knowledge of the  rank-$k$ SVD (or an approximation to it).  Before detailing the method for selecting these indices, we discuss how, given ${\bf p}$ and ${\bf q}$, one should construct ${\bf U}$ so that ${\bf C}{\bf U}{\bf R}$ satisfies desirable approximation properties.

As motivation, suppose for the moment that ${\bf A}$ has exact rank~$k$, and ${\bf C}$ and ${\bf R}$ are full-rank subsets of the columns and rows of ${\bf A}$.\ \ 
Now let ${\bf Y}\in\mbox {\Bb R}^{m\times k}$ and ${\bf Z}\in\mbox {\Bb R}^{n\times k}$ be any matrices that satisfy ${\bf Y}^T {\bf C} = {\bf R} {\bf Z} = {\bf I}\in\mbox {\Bb R}^{k\times k}$.  Then ${\bf C}{\bf Y}^T$ is a projector onto ${\rm Ran}({\bf C})={\rm Ran}({\bf A})$ and $({\bf Z} {\bf R})^T$ is a projector onto ${\rm Ran}({\bf R}^T)={\rm Ran}({\bf A}^T)$, where ${\rm Ran}(\cdot)$ denotes the range (column space).  It follows that ${\bf C}{\bf Y}^T {\bf A} = {\bf A}$ and $({\bf Z}{\bf R})^T{\bf A}^T = {\bf A}^T$.\ \   Putting ${\bf U} \equiv {\bf Y}^T{\bf A}{\bf Z}$ gives 
\[ {\bf C}{\bf U}{\bf R} = {\bf C}{\bf Y}^T{\bf A}{\bf Z}{\bf R} = {\bf A}{\bf Z}{\bf R} = {\bf A}.\]
Thus, any choice of ${\bf Y}$ and ${\bf Z}$ that satisfies ${\bf Y}^T {\bf C} = {\bf R} {\bf Z} = {\bf I}$ gives a ${\bf U}$ such that  ${\bf C}{\bf U}{\bf R}$ exactly recovers ${\bf A}$.\ \ 
In general different choices for ${\bf Y}$ and ${\bf Z}$ give different ${\bf U} = {\bf Y}^T {\bf A} {\bf Z}$.

Now consider the general case~(\ref{eq:CURF}).  Once ${\bf p}$, ${\bf q}$, ${\bf Y}$, and ${\bf Z}$ have been specified, then
\[
   {\bf U} = {\bf Y}^T {\bf A} {\bf Z}  \quad {\rm and } \quad  {\bf F} \equiv {\bf A} - {\bf C} {\bf U} {\bf R}.
\]   
One might design ${\bf Y}$ and ${\bf Z}$ so that ${\bf C}{\bf U}{\bf R}$ matches the selected columns ${\bf C} = {\bf A}(:,{\bf q})$ 
and rows ${\bf R} = {\bf A}({\bf p},:)$ of ${\bf A}$ exactly.   This can be accomplished with \emph{interpolatory projectors}, which we discuss in detail in the next section.  For now, let ${\bf P} = {\bf I}(@:@,{\bf p})\in\mbox {\Bb R}^{m\times k}$ and ${\bf Q} = {\bf I}(@:@,{\bf q})\in\mbox {\Bb R}^{n\times k}$ be submatrices of the identity, so that ${\bf P}^T {\bf a} = {\bf a}({\bf p})$ and ${\bf b}^T {\bf Q} = {\bf b}({\bf q})^T$ for arbitrary vectors
${\bf a}$ and ${\bf b}$ of appropriate dimensions. 
Now define ${\bf Y}^T = ({\bf P}^T {\bf C})^{-1} {\bf P}^T$ and ${\bf Z} = {\bf Q} ({\bf R}{\bf Q})^{-1}$ (presuming ${\bf P}^T{\bf C}$ and ${\bf R}{\bf Q}$ are invertible).  Then since ${\bf C} = {\bf A}(:,{\bf q})$ and ${\bf R} = {\bf A}({\bf p},:)$,
\[
   {\bf P}^T {\bf C} = {\bf C}({\bf p},:) = {\bf A}({\bf p},{\bf q}) 
   \quad \mbox{and}\quad
   {\bf R} {\bf Q} = {\bf R}(:,{\bf q}) = {\bf A}({\bf p},{\bf q}),\]
so
\[ {\bf U} = {\bf Y}^T {\bf A} {\bf Z} 
       = ({\bf P}^T{\bf C})^{-1}{\bf P}^T {\bf A}{\bf Q} ({\bf R}{\bf Q})^{-1} = {\bf A}({\bf p},{\bf q})^{-1}{\bf A}({\bf p},{\bf q}) {\bf A}({\bf p},{\bf q})^{-1} = {\bf A}({\bf p},{\bf q})^{-1}.\]
This CUR approximation matches the ${\bf q}$ columns and ${\bf p}$ rows of ${\bf A}$,
\[
   {\bf A}(:,{\bf q}) = {\bf C} {\bf U} {\bf R}(:,{\bf q}) \ \ {\rm and} \ \
   {\bf A}({\bf p},:) = {\bf C}({\bf p},:) {\bf U} {\bf R},
\]
and, in our experiments, usually delivers a very good approximation.
However, a CUR factorization with better theoretical approximation properties
results from orthogonal projection, as originally suggested by Stewart~\cite[p.~320]{Ste99}; see also, e.g., Mahoney and Drineas \cite{MD09}.
Given a selection of indices ${\bf p}$ and ${\bf q}$, again put
\[
   {\bf C} = {\bf A}(:,{\bf q}) \quad\mbox{and}\quad {\bf R} = {\bf A}({\bf p},:).
\]
Assume that ${\bf C}$ and ${\bf R}$ 
 both have full rank $k$, and now let 
${\bf Y}^T = {\bf C}^I \equiv ({\bf C}^T {\bf C})^{-1} {\bf C}^T $ 
and ${\bf Z} = {\bf R}^I \equiv {\bf R}^T ({\bf R} {\bf R}^T)^{-1} $ denote left and right inverses of ${\bf C}$ and ${\bf R}$.  
These choices also satisfy ${\bf Y}^T{\bf C} = {\bf I}$ and ${\bf R}{\bf Z}={\bf I}$, but now ${\bf C}{\bf Y}^T = {\bf C}\bC^I$ and 
${\bf Z}{\bf R} = {\bf R}^I{\bf R}$ are \emph{orthogonal projectors}. 
We compute
\[
   {\bf U} = {\bf Y}^T {\bf A}{\bf Z} = {\bf C}^I {\bf A} {\bf R}^I,
\]
yielding a CUR factorization that can be viewed as a two step process: 
first the columns of ${\bf A}$ are projected onto Ran(${\bf C}$), 
then the result is projected onto the row space of ${\bf R}$:
\[
   1) \  \  {\bf M} = {\bf C}\bC^I {\bf A}, \quad 2)\ \  {\bf C} {\bf U} {\bf R}  =  {\bf M} {\bf R}^I {\bf R}.
\] 
Both steps are optimal with respect to the 2-norm error, which is the primary
source of the excellent approximation properties of this approach.

Several strategies for selecting ${\bf p}$ and ${\bf q}$ have been
proposed.%
\footnote{In the theoretical computer science literature, one often takes ${\bf C}$ and/or ${\bf R}$ to have rank larger than $k$, but then builds ${\bf U}$ with rank $k$.  By selecting these extra columns and/or rows, one seeks to get within some factor $1+\varepsilon$ of the optimal approximation; see, e.g.,~\cite{BW}.}
The approach presented in the next section is simple to implement
and has complexity $ m@k $ and $n@k $ to select the
indices ${\bf p}$ and ${\bf q}$, provided the leading $k$ right and left singular vectors of ${\bf A}$ are available.
Thus the overall complexity is dominated by the construction of the
rank-$k$ SVD  ${\bf A} \approx {\bf V} {\bf S} {\bf W}^T$, where ${\bf V}^T {\bf V} = {\bf W}^T {\bf W} = {\bf I}\in\mbox {\Bb R}^{k\times k}$
and ${\bf S} = {\rm diag}(\sigma_1, \sigma_2, \ldots , \sigma_{k})$ is the
$k \times k$ matrix of dominant singular values $\sigma_1 \ge \sigma_2  \ge \cdots \ge \sigma_{k}$.

\section{Introduction}
This work presents a new CUR matrix factorization
based upon the Discrete Empirical Interpolation Method (DEIM).\ \ 
A CUR factorization is a low rank approximation
of a matrix ${\bf A} \in \mbox {\Bb R}^{m \times n}$ of the form ${\bf A} \approx {\bf C} {\bf U} {\bf R}$, where
${\bf C} = {\bf A}(:,{\bf q}) \in \mbox {\Bb R}^{m \times k}$ is a subset of the columns of ${\bf A}$ and
${\bf R} = {\bf A}({\bf p},:) \in \mbox {\Bb R}^{k \times n} $ is a subset of the rows of ${\bf A}$.\ \  
(We generally assume $m\ge n$ throughout.)
The $k \times k$ matrix ${\bf U}$ 
is constructed to assure that ${\bf C} {\bf U} {\bf R}$ is a good approximation to  ${\bf A} $.
Assuming the best rank-$k$ singular value decomposition (SVD) ${\bf A} \approx {\bf V} {\bf S} {\bf W}^T$ is available, the algorithm uses the DEIM index selection procedure, ${\bf q} = {\rm DEIM}({\bf V})$ and ${\bf p} = {\rm DEIM}({\bf W})$, to determine ${\bf C}$ and ${\bf R}$.
The resulting approximate factorization is nearly as accurate as the best rank-$k$ SVD, with
\[
  \| {\bf A} - {\bf C} {\bf U} {\bf R} \| \le (\eta_p + \eta_q)\,\sigma_{k+1},
\]
where $\sigma_{k+1}$ is the first neglected singular value of ${\bf A}$,  $\eta_p \equiv \| {\bf V}({\bf p},:\,)^{-1} \|$, and $\eta_q \equiv \| {\bf W}({\bf q},:\,)^{-1} \|$.  

Here and throughout, $\|\cdot\|$ denotes the vector 2-norm and the matrix norm it induces, and $\|\cdot\|_F$ is the Frobenius norm.   We use MATLAB notation to index vectors and matrices, so that, e.g., ${\bf A}({\bf p},:)$ denotes the $k$ rows of ${\bf A}$ whose indices are specified by the entries of the vector ${\bf p}\in\mbox {\Bb N}^k$, while ${\bf A}(:,{\bf q})$ denotes the $k$ columns of ${\bf A}$ indexed by ${\bf q}\in\mbox {\Bb N}^k$.

The CUR factorization is an important tool for handling large-scale data sets, offering two advantages over the SVD:  when ${\bf A}$ is sparse, so too are ${\bf C}$ and ${\bf R}$, unlike the matrices ${\bf V}$ and ${\bf W}$ of singular vectors; and the columns and rows that comprise ${\bf C}$ and ${\bf R}$ are representative of the data (e.g., sparse, nonnegative, integer valued, etc.).\ \  The following simple example, adapted from Mahoney and Drineas~\cite[Fig.~1b]{MD09}, illustrates the latter advantage.  Construct ${\bf A}\in\mbox {\Bb R}^{2\times n}$ so that its first $n/2$ columns have the form
\[ \left[\begin{array}{c} x_1 \\ x_2 \end{array}\right]\]
and the remaining $n/2$ columns have the form
\[    {\sqrt{2}\over2} \left[\begin{array}{cc} -1 & 1 \\ 1 & 1\end{array}\right] 
    \left[\begin{array}{c} x_1 \\ x_2 \end{array}\right],
\]
where in both cases $x_1 \sim N(0,1)$ and $x_2 \sim N(0,4^2)$ are independent samples of normal random variables, i.e.,  
the columns of ${\bf A}$ are drawn from two different multivariate normal distributions.  Figure~\ref{fig:mv_example} shows that the two left singular vectors, though orthogonal by construction, fail to represent the true nature of the data; in contrast, the first two columns selected by the DEIM-CUR procedure give a much better overall representation.  While trivial in this two-dimensional case, one can imagine the utility of such approximations for high-dimensional data.  We shall illustrate the advantages of CUR approximations with further computational examples in Section~\ref{sec:examples}.

\begin{figure}[b!]
\begin{center}
\includegraphics[scale=0.475]{mv_example1}\qquad\quad
\includegraphics[scale=0.475]{mv_example2}
\end{center}

\vspace*{-10pt}
\caption{\label{fig:mv_example} Comparison of singular vectors (left, scaled, in red) 
and DEIM-CUR columns (right, in blue) for a data set drawn from two multivariate normal distributions having different principal axes.}
\end{figure}


CUR-type factorizations originated with ``pseudoskeleton'' approximations~\cite{GTZ97} and pivoted, truncated QR decompositions~\cite{Ste99}; in recent years many new algorithms have been proposed in the numerical linear algebra and theoretical computer science literatures. 
Some approaches seek to maximize the \emph{volume} of the decomposition~\cite{GTZ97,TKB12}.  
Numerous other algorithms instead use \emph{leverage scores}~\cite{BW,DMM08,MD09,WZ13}.  These methods typically first compute a singular value decomposition%
\footnote{We use the nonstandard notation ${\bf V} {\bf S} {\bf W}^T$ for the SVD to avoid conflicts with ${\bf U}$ in the standard CUR notation.} 
${\bf A} = {\bf V}{\bf S}{\bf W}^T$ (or an approximation to it), with ${\bf V}\in\mbox {\Bb R}^{m\times n}$, ${\bf W}\in\mbox {\Bb R}^{n\times n}$.\ \ 
The leverage score for the $j$th row ($k$th column) of ${\bf A}$ is the squared two-norm of the $j$th row of ${\bf V}$ ($k$th row of ${\bf W}$).  When scaled by the number of singular vectors, these leverage scores give probability distributions for randomly sampling the columns and rows to form ${\bf C}$ and ${\bf R}$.  This approach leads to probabilistic bounds on $\|{\bf A} - {\bf C} {\bf U} {\bf R} \|_F$~\cite{DMM08,MD09}.  In cases where ${\bf A}$ has small singular values (precisely the case where one would seek a low-rank factorization), the singular vectors can be sensitive to perturbations to ${\bf A}$, making the leverages scores unstable~\cite{IW}.  Thus leverage scores are often computed using only the leading few singular vectors, but the choice of how many vectors to keep can be somewhat ad hoc.

The algorithm described in Sections~\ref{sec:factor} and~\ref{sec:deim} is entirely deterministic and involves few (if any) parameters.  The method is supported by an error analysis in Section~\ref{sec:theory} that also applies to a broad class of CUR factorizations.  
This section includes an improved bound on the error constants 
$\eta_p$ and $\eta_q$ for DEIM row and column selection, which also applies to  the analysis of DEIM-based model order reduction~\cite{CS10}.
In Section~\ref{sec:qr} we propose a novel incremental QR algorithm for approximating the SVD
(and potentially also approximating leverage scores).  Section~\ref{sec:examples} illustrates the performance of this new CUR factorization on several examples. 


In many applications one cares primarily about key columns or rows of ${\bf A}$, rather than an explicit ${\bf A}={\bf C}{\bf U}{\bf R}$ factorization.
The DEIM technique, which identifies rows and columns of ${\bf A}$ independently, can easily be used to 
select only columns or rows, leading to an ``interpolatory decomposition'' of the form 
${\bf A} = {\bf C}\widehat{{\bf U}}$ or ${\bf A} = \widehat{{\bf U}}{\bf R}$;
such factorizations have the advantage that $\widehat{{\bf U}}$ can be much better conditioned than the ${\bf U}$ matrix in the CUR factorization.  For further details about general interpolatory decompositions, see~\cite[\S1]{CGMR05}.

\section{Incremental QR Factorization} \label{sec:qr}

The DEIM point selection process presumes access to the first $k$ left and right singular vectors of ${\bf A}\in\mbox {\Bb R}^{m\times n}$.  If either $m$ or $n$ is of modest size (say ${}\le 1000$) and ${\bf A}$ can be stored as a dense matrix, \ library software for computing the ``economy sized'' SVD, e.g., \verb|[V,S,W] = svd(A,'econ')| in MATLAB, usually performs very well. 
For larger scale problems, the leading $k$ singular vectors can be computed using iterative methods, such as the Krylov subspace-based ARPACK software~\cite{LSY98} (used by MATLAB's {\tt svds} command), PROPACK~\cite{Lar05}, IRLBA~\cite{BR13}, or the Jacobi--Davidson algorithm~\cite{Hoc01}.  Randomized SVD algorithms provide an appealing alternative with probabilistic error bounds~\cite{HMT11}. 
Here we describe another approach that satisfies a deterministic error bound (Lemma~\ref{lem:incqr}) and \emph{only requires one pass through the matrix ${\bf A}$}, a key property for massive data sets that cannot easily be stored in memory.



\begin{algorithm}[h!]
\begin{center}
\fbox{\begin{tabular}{p{10cm}}
\ \\[-25pt]
\begin{tabbing}
defghasijkl\=bbbb\=cccc\= \kill
$ \ \ $ {\bf Input:} \> $ {\bf A}$,  an $m \times n$ matrix \\
$ \ \ $ \> {\it tol}, a positive scalar controlling the accuracy of the factorization\\[5pt]
$\  $  {\bf Output:}  \> $\bQ$, an $m \times k$ matrix with orthonormal columns\\
$\ \ \ \ \ \ \ $ \> $\bR$, a $k \times n$ rectangular matrix  \\
$\ \ \ \ \ \ \ $ \>\>\>  with ${\bf A} \approx  \bQ \bR$ \\[5pt]
$\ \ \ $ Choose $k \ll \min(m,n)$\\
$\ \ \ $ Compute the QR factorization ${\bf A}(:,1:k) = \bQ \bR$, with $\bQ \in \mbox {\Bb R}^{n\times k}$ and $\bR\in\mbox {\Bb R}^{k\times m}$\\  
$ \ \ \ $  rownorms($i$) = $\|\bR(i,:)\|^2$ for $i = 1, \ldots, k$ \\
$\ \ \ $   $j = k+1$\\
$\ \ \ $ {\bf while} $j \le n$\\
$ \ \ \ $  \> ${\bf a} = {\bf A}(@:@,j) ; \ \br = \bQ^T{\bf a}; \ \  {\bf f} = {\bf a} - \bQ \br; \ \   \rho = \|{\bf f} \|; \ \ \bq = {\bf f}/\rho$\\
$\ \ \ $  \>   $\bQ = [\bQ, \ \bq];  \ \   \bR = \left[ \begin{array}{cc} \bR  &  \br\\
                                                                {\bf 0} & \rho \\
                                               \end{array} \right]$ \\[5pt]
$\ \ \ $  \> rownorms($i$) = rownorms($i$) + $\br(i)^2$ for $i = 1,\ldots, k$ \\                                               
$\ \ \ $  \> rownorms($k+1$) ${}= \rho^2; \ \  $ \\
$ \ \ \ $ \>  FnormR = {\rm sum}(rownorms); \\
$ \ \ \ $ \> $[\sigma, i_{\rm min}] = {\rm min}({\rm rownorms}(1:k+1));$ \\[5pt]
$ \ \ \ $ \> {\bf if} $\sigma > ({\it tol}^2)*({\rm FnormR} - {\rm rownorms}(i_{\rm min})$ \\
$ \ \ \ $ \>\> \% \emph{no deflation} \\
$ \ \ \ $ \>\> $ k = k + 1;$ \\
$ \ \ \ $ \> {\bf else} \% \emph{deflation required}\\
$ \ \ \ $ \>\> {\bf if} $i_{\rm min}  < k+1$\\
$ \ \ \ $ \>\>\> $\bR(i_{\rm min},@:@) = \bR(k+1,@:@); \ \bQ(@:@,i_{\rm min}) = \bQ(@:@, k+1)$ \\
$ \ \ \ $ \>\>\> $\ {\rm rownorms}(i_{min}) = {\rm rownorms}(k+1) $ \\
$ \ \ \ $ \>\> {\bf end} \\
$ \ \ \ $ \>\> \%  \emph{delete the minimum norm row of $\bR$}\\ 
$ \ \ \ $ \>\> $\bQ = \bQ(@:@,1:k); \bR = \bR(1:k,@:@) $ \\
$ \ \ \ $ \> {\bf end} \\
$ \ \ \ $ \> $j = j+1$\\
$ \ \ \ $  {\bf end} \\
\end{tabbing}\\[-23pt]\end{tabular}}
\end{center}

\vspace*{-10pt}
\caption{Incremental QR low rank approximate factorization}
\label{fig:QRincremental}
\end{algorithm}

This approach is based on an {\it incremental} low
rank ${\bf A}\approx{\bf Q} \bR $ approximation, where $\bQ\in\mbox {\Bb R}^{n\times k}$ has orthonormal columns and $\bR\in\mbox {\Bb R}^{k\times m}$ is upper triangular.  
(In this section only, $\bQ$ and $\bR$ denote different quantities from elsewhere in the paper.)
Take the dense (economy sized) SVD $\bR = \widehat{{\bf V}} {\bf S} {\bf W}^T$, and put ${\bf V} = \bQ \widehat{{\bf V}} $ to get
\begin{equation} \label{eq:AQRSVD}
   {\bf A} \approx  \bQ \bR = {\bf V} {\bf S} {\bf W}^T.
\end{equation}
Incremental algorithms for building the QR factorization and SVD
have been proposed by Stewart~\cite{Ste99}, Baker, Gallivan, and Van Dooren~\cite{BGV12} and many others, as
surveyed in~\cite{Bak04}; these ideas are also closely related to rank-revealing QR factorizations~\cite{GE96}.  Algorithm~\ref{fig:QRupdate} differs from those of Stewart in its use of internal pivoting and threshold truncation in place of Stewart's column pivoting.   This distinction enables a
one-pass algorithm that is closely related to~\cite[Algorithm~1]{BGV12}.

The proposed method is presented in Algorithm~\ref{fig:QRincremental}, which proceeds at each step by orthogonalizing a column of ${\bf A}$ against the previously orthogonalized columns. 
The rank of the resulting factors is controlled through an update-and-delete procedure that is illustrated in Figure~\ref{fig:QRupdate}.  
After orthogonalizing a column of ${\bf A}$, 
the algorithm checks if any row of $\bR$ has small relative norm; 
if such a row exists, the corresponding column of ${\bf Q}$ 
makes little contribution to the factorization, so that column of ${\bf Q}$ and row of $\bR$ can be deleted at only a small loss of accuracy in the factorization.  
(Future columns of ${\bf A}$ will not be orthogonalized
against the vector deleted from ${\bf Q}$, so this direction can re-emerge if a later column in ${\bf A}$ warrants it.)

Robust implementations of Algorithm~\ref{fig:QRupdate} should replace the classical Gram--Schmidt operations 
\[
   \br = \bQ^T {\bf a}, \ \  {\bf f} = {\bf a} - \bQ \br
\]
with a re-orthogonalization step, as 
suggested by Daniel, Gragg, Kaufman, and Stewart~\cite{DGKS}:
\begin{eqnarray}
  \br &=& \bQ^T {\bf a} \nonumber \\
  {\bf f} &=& {\bf a} - \bQ \br \nonumber \\
       &\ & {\bf c} = \bQ^T {\bf f} \label{eq:reorth1} \\
       &\ &  {\bf f} = {\bf f} - \bQ{\bf c} \label{eq:reorth2}\\
       &\ &  \br = \br + {\bf c} \label{eq:reorth3}\\
  \rho &=& \| {\bf f} \| \nonumber \\
  \bq  &=&  {\bf f}/\rho. \nonumber 
\end{eqnarray}
The extra steps~(\ref{eq:reorth1})--(\ref{eq:reorth3}) generally provide a $\bQ$ that is numerically orthogonal to working precision.
Pathological cases are easily overcome with some additional slight modifications;
see~\cite{GIR05} for a complete analysis.
Because this algorithm uses the classical Gram--Schmidt method, one can easily block it for parallel efficiency.


\begin{figure}[t!]
\begin{center}
\begin{picture}(300,50)
\put(-40,-20){${\bf A}(@:@,1:j) = {}$}
\put(-25,-80){
\begin{minipage}[t]{2.5in}\begin{center} \hspace*{30pt}\includegraphics[scale=0.35]{QRupdateA}\\  Partial ${\bf Q}\bR$ factorization\end{center} \end{minipage}}

\put(160,-20){\small ${\bf A}(:,1\!:\!j\!+\!1) = {}$}
\put(177,-80){\begin{minipage}[t]{2.5in}\begin{center}\hspace*{30pt} \includegraphics[scale=0.35]{QRupdateB}\\  Extend with Gram--Schmidt\end{center} \end{minipage}}

\put(-95,-230){\begin{minipage}[t]{3in}\begin{center} \includegraphics[scale=0.35]{QRupdateC}\\  Find $\bq_i$ with $\|\bR(i,:)\|^2 < \tol^{@2}\,(\|\bR\|_F^2 - \|\bR(i,:)\|^2)$\end{center} \end{minipage}}

\put(90,-230){\begin{minipage}[t]{2in}\begin{center} \includegraphics[scale=0.35]{QRupdateD}\\ Replace $\bq_i$, $\bR(i,:)$ \end{center} \end{minipage}}

\put(255,-230){\begin{minipage}[t]{2in}\begin{center} \includegraphics[scale=0.35]{QRupdateE}\\ Truncate last column\\ of $\bQ$ and row of $\bR$\end{center} \end{minipage}}
\end{picture}
\vspace*{3.35in}
\end{center}
\caption{\label{fig:QRupdate}
Diagram illustrating the QR update procedure.}
\end{figure}




\subsection{Incremental QR Error Bounds}

At step $j$ the truncation criterion in Algorithm~\ref{fig:QRincremental} will delete row $\br_i^T =  {\bf e}_i^T \bR_j \ $ if
\[
 \| \br_i \| \le \tol\, \| \widehat{\bR}_j \|_F,
\]
where $\br_i^T$ is the row of minimum norm and $\widehat{\bR}_j$ denotes $\bR_j$ with the $i$th row deleted.  This strategy has a straightforward error analysis, 
which, in light of the approximation~(\ref{eq:AQRSVD}), also 
implies an error bound on the resulting SVD.

\begin{lemma}\label{lem:incqr}
Let $\bR_j$ be the triangular factor at step $j$ of 
Algorithm~\ref{fig:QRincremental}, and $\bQ_j $ the corresponding orthonormal columns in the approximate QR factorization ${\bf A}_j \approx \bQ_j \bR_j$, where ${\bf A}_j$ consists of
the first $j$ columns of ${\bf A}$. Then
\[
 \| {\bf A}_j - \bQ_j \bR_j \|_F \le \tol \cdot d_j \cdot \| \bR_j \|_F, 
 \]
 where $d_j$ is the number of coloumn/row deletions that have been made up to and including step $j$.
 (Note that $\bQ_j \in \mbox {\Bb R}^{m\times (j-d_j)}$, $\bR\in\mbox {\Bb R}^{(j-d_j)\times j}$, and $d_n = \min\{m,n\} - k$, where $k =  { {{\rm rank }\,} } (\bQ_n\bR_n)$.)
\end{lemma}

Before proving this lemma, we note that it gives a bound on the error in the resulting
approximate SVD of ${\bf A}$.  Suppose $d_n$ deletions are made when this algorithm computes
the approximate factorization ${\bf A} \approx {\bf Q}{\bf R}$ with tolerance $\tol$. 
Given the SVD ${\bf R} = \widehat{{\bf V}}{\bf S}{\bf W}^*$, set ${\bf V} \equiv {\bf Q}\widehat{{\bf V}}$.
Then
\[ \|{\bf A} - {\bf V}{\bf S}{\bf W}^*\|_F \le \tol \cdot d_n \cdot  \|\bR\|_F.\]
{\it Proof of Lemma~\ref{lem:incqr}:}
The proof shall be by induction.  Let ${\bf E}_j = {\bf A}_j - \bQ_j  \bR_j$ and assume
\begin{equation} \label{eq:inductive}
\| {\bf E}_j \|_F \le \tol \cdot d_j \cdot \| \bR_j \|_F.
\end{equation}
Orthogonalize column $j+1$ of ${\bf A}$ using Gram--Schmidt to obtain
\[
 {\bf A}_{j+1} = \bQ_{j+1}\bR_{j+1} + [{\bf E}_j, {\bf 0}].
\]
If no deflation occurs at this step, the bound holds trivially
since 
\[
   \| {\bf E}_{j+1} \|_F = \|[{\bf E}_j, {\bf 0}]\|_F \le tol\cdot d_j \cdot \|\bR_j\|_F
   \le \tol \cdot d_{j+1} \cdot \| \bR_{j+1} \|_F,
\]
because $d_{j+1}= d_j$ and $\|\bR_j \|_F \le \| \bR_{j+1}\|_F$.

Suppose $\bR_j$ has dimension $k \times j$ (i.e., $k=j-d_j$).
Let $i$ be the index of the row of minimum norm and let $\widehat{\bR}_{j+1}$
be obtained by deleting the $i$th row of $\bR_{j+1}$.\ \ 
If $\br_i^T  = {\bf e}_i^T \bR_{j+1} $ satisfies $\|\br_i^T \| \le tol\cdot \| \widehat{\bR}_{j+1} \|_F$
then deflation occurs.
Deleting column $i$ of $\bQ_{j+1}$ and row $i$ of $\bR_{j+1}$
replaces $\bQ_{j+1}$ and $\bR_{j+1}$ with $\widehat{\bQ}_{j+1}$ and $\widehat{\bR}_{j+1}$.  Then
\[
\widehat{\bQ}_{j+1} \widehat{\bR}_{j+1} = \bQ_{j+1} ( \bR_{j+1} - {\bf e}_i^{} \br_i^T),
\]
and 
\[
  {\bf A}_{j+1} = \bQ_{j+1} (\bR_{j+1} - {\bf e}_i^{} \br_i^T) + [{\bf E}_j, {\bf 0}] + \bQ_{j+1}{\bf e}_i^{} \br_i^T.
\]
Hence the deletion gives the overall error
\[
  {\bf E}_{j+1} = {\bf A}_{j+1} - \widehat{\bQ}_{j+1} \widehat{\bR}_{j+1} = [{\bf E}_j, {\bf 0}] + 
   \bQ_{j+1}{\bf e}_i^{}\br_i^T.
\]
Therefore, when $i < k+1$, the inductive assumption~(\ref{eq:inductive}) implies
\[
   \| {\bf E}_{j+1}\|_F \le \| {\bf E}_j \|_F + \|\br_i^T \| \le \tol \cdot ( d_j \cdot \| \bR_j \|_F + \|\widehat{\bR}_{j+1} \|_F) \le \tol \cdot (d_j+1) \cdot \| \widehat{\bR}_{j+1}\|_F,
\]
since $\widehat{\bR}_{j+1}$ contains row $k+1$ of $\bR_{j+1}$, 
which must have a norm larger than the row marked for deletion.   
Since row $k+1$ of $\widehat{\bR}_{j+1}$ consists of just one
nonzero element, 
\[
 \| \widehat{\bR}_{j+1} \|^2_F \ge \| \widehat{\bR}_j \|_F^2 + \rho_{k+1,j+1}^2 \ge \| \bR_j \|_F^2,
 \]
 where $\rho_{k+1,j+1}$ is the element  $\bR_{j+1}(k+1,j+1) $  and $\widehat{\bR}_j$
 is the matrix $\bR_j$ with $i$th row deleted.
 If $i=k+1$, then the last row of $\bR_{j+1} $ is deleted and the desired inequality must hold, since $\bR_j$ is a submatrix of $\widehat{\bR}_{j+1}$.  At the end of this process, replace $\bR_{j+1}$ and $\bQ_{j+1}$ with $\widehat{\bR}_{j+1}$ and $\widehat{\bQ}_{j+1}$ to obtain the approximation
 \[ \|{\bf A}_{j+1}\| \le \tol \cdot d_{j+1}\cdot\|\bR_{j+1}\|_F,\]
 since $d_{j+1} = d_j + 1$.

The error bound for the base case $j=1$ clearly holds, completing the induction. { { \hfill\rule{3mm}{3mm} } } 

\medskip
The approximate QR factorization that results from this algorithm could be used directly for the approximation of leverage scores.  The perturbation theory of Ipsen and Wentworth~\cite{IW} describes how the tolerance in our algorithm will affect the accuracy of the resulting leverage scores.
We also note that for extra expediency this one-pass QR algorithm could be stopped when $\|\widehat{\bR}_j\|_F \approx \|{\bf A}\|_F$ (at the cost of an extra pass through ${\bf A}$ to compute $\|{\bf A}\|_F$), or applied to only a random sampling of $k$ columns of ${\bf A}$.  (Drma\v{c} and Gugercin propose a different random approach to DEIM index selection, based on sampling rows of ${\bf V}$ to compute DEIM indices~\cite{DG}.)

\section{CUR Approximation Properties} \label{sec:theory}

While the theory presented in this section was designed to
bound $\|{\bf A}-{\bf C}{\bf U}{\bf R}\|$ for the DEIM-CUR method, 
the analysis applies to \emph{any} CUR factorization with 
full rank ${\bf C} \in \mbox {\Bb R}^{m\times k}$ and 
${\bf R}\in \mbox {\Bb R}^{k\times n}$, and ${\bf U} = {\bf C}^I{\bf A}{\bf R}^I$, 
regardless of the procedure used for selecting the columns 
and rows.%
\footnote{We are grateful to Ilse Ipsen for noting
the applicability of this analysis to all such ${\bf C}{\bf U}{\bf R}$ factorizations,
and for also pointing out that, given knowledge of \emph{all} the singular values and vectors of ${\bf A}$, our Lemma~\ref{ProjErr} can be sharpened via application of~\cite[Thm.~9.1]{HMT11}.  Indeed, Ipsen observes that the interpolatory projector proof of Lemma~\ref{ProjErr} can be adapted to simplify the multipage proof of \cite[Thm.~9.1]{HMT11}.}

Consider a CUR factorization that uses row indices ${\bf p} \in \mbox {\Bb N}^k$ and column indices ${\bf q}\in\mbox {\Bb N}^k$, and set
\[ {\bf P} = {\bf I}(@@:@@,{\bf p}) = [{\bf e}_{p_1}, \ldots, {\bf e}_{p_k}]\in\mbox {\Bb R}^{m\times k}, \qquad
    {\bf Q} = {\bf I}(@@:@@,{\bf q}) = [{\bf e}_{q_1}, \ldots, {\bf e}_{q_k}]\in\mbox {\Bb R}^{n\times k}.
 \]
The first step in this analysis bounds the mismatch
between ${\bf A}$ and its interpolatory projection $\mbox{\boldmath$\cal P$}{\bf A}$.



\begin{lemma}
\label{DEIM_bound}
Assume ${\bf P}^T{\bf V}$ is invertible and let $\mbox{\boldmath$\cal P$} = {\bf V} ({\bf P}^T{\bf V})^{-1}{\bf P}^T$ be the interpolatory
projector~$(\ref{eq:intproj})$. If ${\bf V}^T{\bf V}={\bf I}$,  then  any ${\bf A}\in\mbox {\Bb R}^{m\times n}$ satisfies
\[
  \| {\bf A} - \mbox{\boldmath$\cal P$}{\bf A} \| \le \|({\bf P}^T {\bf V})^{-1} \| \|({\bf I} - {\bf V} {\bf V}^T){\bf A} \|.
\]
Additionally, if ${\bf V}$ consists of the leading $k$ left singular vectors of ${\bf A}$, then
\[
  \| {\bf A} - \mbox{\boldmath$\cal P$}{\bf A} \|  =  \|({\bf I} - \mbox{\boldmath$\cal P$}) {\bf A}\| \le \|({\bf P}^T {\bf V})^{-1} \|\, \sigma_{k+1}.
\]
\end{lemma}


\noindent{\it Proof:} First note 
that $\mbox{\boldmath$\cal P$} {\bf V} = {\bf V} ({\bf P}^T {\bf V})^{-1} {\bf P}^T {\bf V} = {\bf V}$, so that $({\bf I} - \mbox{\boldmath$\cal P$}){\bf V} = {\bf 0}$.
Therefore
\[
  \| {\bf A} - \mbox{\boldmath$\cal P$}{\bf A} \|  = \|({\bf I} - \mbox{\boldmath$\cal P$}) {\bf A}\| 
  = \|({\bf I} - \mbox{\boldmath$\cal P$})({\bf I} - {\bf V} {\bf V}^T) {\bf A}\| \le 
   \|({\bf I} - \mbox{\boldmath$\cal P$})\| \|({\bf I} - {\bf V} {\bf V}^T) {\bf A}\|.
\]
It is well known that 
\[
   \|{\bf I} - \mbox{\boldmath$\cal P$}\| = \| \mbox{\boldmath$\cal P$} \|  = \|({\bf P}^T {\bf V})^{-1} \| 
\]
so long as $\mbox{\boldmath$\cal P$} \ne {\bf 0} \ {\rm or} \ {\bf I}$; see, e.g.,~\cite{Szy06}.
This establishes the first result.  The second follows from the fact that
\[
    \|({\bf I} - {\bf V} {\bf V}^T) {\bf A}\| = \| {\bf A} - {\bf V} {\bf S} {\bf W}^T\| = \sigma_{k+1} 
\]
when ${\bf V}$ consists of the leading $k$ left singular vectors of ${\bf A}$.
 { { \hfill\rule{3mm}{3mm} } } 

\medskip
Now let ${\bf V} {\bf S} {\bf W}^T \approx {\bf A}$ be a rank-$k$ SVD of ${\bf A}$.   
(The singular vectors play a crucial role in this analysis, even if ${\bf p}$ and ${\bf q}$ were selected using some scheme that did not reference them.)
In addition to the interpolatory projector 
$\mbox{\boldmath$\cal P$} = {\bf V}({\bf P}^T{\bf V})^{-1}{\bf P}^T$ that operates 
on the left of ${\bf A}$, we shall also use 
$\mbox{\boldmath$\cal Q$} = {\bf Q}({\bf W}^T {\bf Q})^{-1} {\bf W}^T$, which operates
on the right of ${\bf A}$.\ \ 
Assuming that ${\bf P}^T{\bf V}$ and ${\bf W}^T{\bf Q}$ are invertible, define the error constants
\[ \eta_p \equiv \|({\bf P}^T{\bf V})^{-1}\|, \qquad
   \eta_q \equiv \|({\bf W}^T{\bf Q})^{-1}\|. \]
Lemma~\ref{DEIM_bound} implies
\begin{equation} \label{eq:iebnd}
   \| {\bf A} ({\bf I} - \mbox{\boldmath$\cal Q$}) \| \le \eta_q @@\sigma_{k+1}
   \quad {\rm and } \quad
      \| ({\bf I} - \mbox{\boldmath$\cal P$}) {\bf A} \| \le \eta_p @@\sigma_{k+1}.
\end{equation}
The next lemma shows that these bounds on the
error of the interpolatory projection of ${\bf A}$ onto the select columns and rows also apply to the \emph{orthogonal} projections of ${\bf A}$ onto the same column and row spaces.

\begin{lemma}
\label{ProjErr}
Suppose the row and column indices ${\bf p}$ and ${\bf q}$ give
full rank matrices ${\bf C} = {\bf A}(@@:@@,{\bf q}) = {\bf A}{\bf Q}\in\mbox {\Bb R}^{m\times k}$ and ${\bf R} = {\bf A}({\bf p},@@:@@)={\bf P}{\bf A}\in\mbox {\Bb R}^{k\times n}$, 
with finite error constants $\eta_p$ and $\eta_q$, and suppose that $k<\min\{m,n\}$.  Then
\[
   \| ({\bf I} - {\bf C} {\bf C}^I) {\bf A} \| \le \eta_q @@\sigma_{k+1}
   \quad {\rm and}\quad
   \| {\bf A}({\bf I} - {\bf R}^I {\bf R}) \| \le \eta_p @@\sigma_{k+1}.
\]
\end{lemma}
{\it Proof:}
Using the formula ${\bf C}= {\bf A}{\bf Q}$, we have 
${\bf C}^I = ({\bf C}^T{\bf C})^{-1}{\bf C}^T = ({\bf Q}^T{\bf A}^T{\bf A}{\bf Q})^{-1}({\bf A}{\bf Q})^T$, 
so the orthogonal projection of ${\bf A}$ onto ${\rm Ran}({\bf C})$ is
\[
    {\bf C} {\bf C}^I {\bf A}  = ({\bf A} {\bf Q}({\bf Q}^T {\bf A}^T {\bf A} {\bf Q})^{-1} {\bf Q}^T {\bf A}^T){\bf A}
      = {\bf A} ({\bf Q} ({\bf Q}^T {\bf A}^T {\bf A} {\bf Q})^{-1} {\bf Q}^T {\bf A}^T{\bf A}).
\]
Hence the error in the orthogonal projection of ${\bf A}$ is
\begin{equation} \label{eq:lemid1}
    ({\bf I} - {\bf C} {\bf C}^I) {\bf A} = {\bf A}({\bf I} - \mbox{\boldmath$\Phi$}), \ \ {\rm where} \ \ \mbox{\boldmath$\Phi$} =  {\bf Q}({\bf Q}^T {\bf A}^T {\bf A} {\bf Q})^{-1} {\bf Q}^T {\bf A}^T{\bf A}.
\end{equation}
Note that $\mbox{\boldmath$\Phi$}$ is an oblique projector onto ${\rm Ran}({\bf Q})$, so  $\mbox{\boldmath$\Phi$} {\bf Q} = {\bf Q}$.  
Therefore, $\mbox{\boldmath$\Phi$} \mbox{\boldmath$\cal Q$} = \mbox{\boldmath$\cal Q$}$, since
\[
\mbox{\boldmath$\Phi$} \mbox{\boldmath$\cal Q$} =  \mbox{\boldmath$\Phi$} {\bf Q} ({\bf W}^T{\bf Q})^{-1} {\bf W}^T =  {\bf Q} ({\bf W}^T{\bf Q})^{-1} {\bf W}^T = \mbox{\boldmath$\cal Q$}.
\]
This implies that
\[
  {\bf A}({\bf I} - \mbox{\boldmath$\Phi$}) = {\bf A}({\bf I} - \mbox{\boldmath$\Phi$})({\bf I} - \mbox{\boldmath$\cal Q$}) =  ({\bf I} - {\bf C} {\bf C}^I) {\bf A} ({\bf I} - \mbox{\boldmath$\cal Q$}),
\]
and so from~(\ref{eq:lemid1}) we have
\begin{eqnarray*}
\|({\bf I} - {\bf C} {\bf C}^I) {\bf A} \| &=& \| {\bf A}({\bf I} - \mbox{\boldmath$\Phi$})\| \\
     &=& \|  ({\bf I} - {\bf C} {\bf C}^I) {\bf A} ({\bf I} - \mbox{\boldmath$\cal Q$})\| \\
     &\le& \|  {\bf I} - {\bf C} {\bf C}^I \| \|{\bf A} ({\bf I} - \mbox{\boldmath$\cal Q$})\| \\
     &\le& \eta_q@@\sigma_{k+1}.
\end{eqnarray*}
The last line follows from the bound~(\ref{eq:iebnd}) and the fact that $ \|  {\bf I} - {\bf C} {\bf C}^I \| = 1$, since ${\bf C}\bC^I$ is an orthogonal projector and $k<\min\{m,n\}$.

A similar argument shows that 
\[
  {\bf A}({\bf I} - {\bf R}^I {\bf R}) = ({\bf I} - \mbox{\boldmath$\Psi$}) {\bf A}  
\]
where $\mbox{\boldmath$\Psi$} = {\bf A} {\bf A}^T {\bf P} ({\bf P}^T{\bf A}\bA^T{\bf P})^{-1}{\bf P}^T$, and also that
\[
  ({\bf I} - \mbox{\boldmath$\Psi$}){\bf A} = ({\bf I} - \mbox{\boldmath$\cal P$})({\bf I} - \mbox{\boldmath$\Psi$}){\bf A} = 
({\bf I} - \mbox{\boldmath$\cal P$}) {\bf A} ({\bf I} - {\bf R}^I {\bf R}),
\]
from which follows the error bound
\[
  \|{\bf A}({\bf I} - {\bf R}^I {\bf R})\|   \le \|({\bf I} - \mbox{\boldmath$\cal P$}) {\bf A} \| \|{\bf I} - {\bf R}^I {\bf R}\| \le \eta_p @@\sigma_{k+1}.\qquad  { { \hfill\rule{3mm}{3mm} } } 
\]

\medskip
The main result on approximation of ${\bf A}$ by ${\bf C}{\bf U}{\bf R}$
readily follows from combining this last lemma with a basic CUR analysis technique used by Mahoney and Drineas~\cite[eq.~(6)]{MD09}.
 

\begin{theorem}
\label{basic_bound}
Given ${\bf A}\in\mbox {\Bb R}^{m\times n}$ and $1\le k<\min\{m,n\}$, let ${\bf C} = {\bf A}(@@:@@,{\bf q})\in\mbox {\Bb R}^{m\times k}$ and ${\bf R} = {\bf A}({\bf p},@@:@@)\in\mbox {\Bb R}^{k\times n}$ 
with finite error constants $\eta_p$ and $\eta_q$, and set 
${\bf U} = {\bf C}^I{\bf A}{\bf R}^I$.  Then
\[
   \| {\bf A} - {\bf C} {\bf U} {\bf R} \| \le (\eta_p + \eta_q)@@ \sigma_{k+1}.
\]
\end{theorem}
\noindent{\it Proof:} 
From the definitions,
\[
  {\bf A} - {\bf C} {\bf U} {\bf R} = {\bf A} - {\bf C} {\bf C}^I {\bf A} {\bf R}^I {\bf R} 
                    = ({\bf I} - {\bf C} {\bf C}^I ){\bf A} + {\bf C}\bC^I {\bf A} ({\bf I} - {\bf R}^I {\bf R}).
\]
Applying Lemma~\ref{ProjErr},
\begin{eqnarray*}
  \|{\bf A} - {\bf C} {\bf U} {\bf R} \| 
                    &\le& \|({\bf I} - {\bf C} {\bf C}^I ){\bf A} \| +  \|{\bf C}\bC^I \| \| {\bf A} ({\bf I} - {\bf R}^I {\bf R}) \| \\
                    &\le& \eta_q@@\sigma_{k+1} +  \eta_p@@ \sigma_{k+1}\\
                    &=& ( \eta_p  +  \eta_q )@@\sigma_{k+1},
\end{eqnarray*}
since $\|{\bf C}\bC^I\| = 1$.   { { \hfill\rule{3mm}{3mm} } } 

\medskip
Theorem~\ref{basic_bound} shows that ${\bf C}{\bf U}{\bf R}$ is within 
a factor of $\eta_p +\eta_q$ of the optimal rank-$k$ 
approximation, hence these error constants suggest a 
way to assess a wide variety of column/row selection 
schemes. 
The quality of the approximation is controlled by the conditioning of the selected $k$ rows of the dominant $k$ (exact) singular vectors.
If those singular vectors are available as part of the
column/row selection process, then Theorem~\ref{basic_bound}
provides an \emph{a posteriori} bound requiring only the fast ($ {\cal O} (k^3)$) computation of $\eta_p$ and $\eta_q$, and thus
could suggest methods for adjusting either $k$ or the point selection process to reduce the error constants.  In this context, notice that if ${\bf V}{\bf S}{\bf W}^T$ is only an \emph{approximation} to the optimal rank-$k$ SVD with ${\bf V}$ and ${\bf W}$ having orthonormal columns (as computed, for example, using the incremental QR algorithm described in the next section), the preceding analysis gives
\begin{eqnarray}
\|{\bf A} - {\bf C}{\bf U}{\bf R}\| 
   \!\!&\le&\!\! \|({\bf I}-{\bf C}\bC^I){\bf A}\| + \|{\bf A}({\bf I}-{\bf R}^I{\bf R})\| \nonumber \\
   \!\!&=&\!\! \|{\bf A}({\bf I}-\mbox{\boldmath$\cal Q$})\| + \|({\bf I}-\mbox{\boldmath$\cal P$}){\bf A}\| \nonumber \\
   \!\!&\le&\!\! \|({\bf W}^T{\bf Q})^{-1}\| \|{\bf A}({\bf I}-{\bf W}\bW^T)\| + \|({\bf P}^T{\bf V})^{-1}\| \|({\bf I}-{\bf V}\bV^T){\bf A}\|, \label{eq:approxbnd}
\end{eqnarray}
showing how $\sigma_{k+1}$ in Theorem~\ref{basic_bound} is replaced by the error in the approximate SVD through $\|{\bf A}({\bf I}-{\bf W}\bW^T)\|$ and $\|({\bf I}-{\bf V}\bV^T){\bf A}\|$. 
In this case $\|({\bf W}^T{\bf Q})^{-1}\|$ and $\|({\bf P}^T{\bf V})^{-1}\|$ are computed using the approximate singular vectors in ${\bf V}$ and ${\bf W}$, rather than the exact singular vectors in the theorem.
Alternatively, if one has probabilistic bounds for $\eta_p$ and $\eta_q$, 
then Theorem~\ref{basic_bound}  immediately gives a probabilistic bound for $\| {\bf A} - {\bf C} {\bf U} {\bf R} \|$.

Numerical examples in Section~\ref{sec:examples} compare how the error constants evolve as $k$ increases for the DEIM-CUR 
factorization and several other factorizations based on leverage scores.  

\subsection{Interpretation of the bound for DEIM-CUR}

For DEIM-CUR, we can ensure  the hypotheses of Theorem~\ref{basic_bound} are satisfied and bound the error constants.  Suppose the DEIM points are selected using the exact rank-$k$ SVD ${\bf A} \approx {\bf V}{\bf S}{\bf W}^T$.\ \ Lemma~\ref{DEIM_inv} ensures that the matrices
${\bf P}^T{\bf V}$ and ${\bf W}^T{\bf Q}$ are invertible, so $\eta_p$ and $\eta_q$ are finite.  The DEIM strategy also gives full rank ${\bf C}$ and ${\bf R}$ matrices, presuming $k\le {\rm rank}({\bf A})$.  To see this, note that for any unit vector ${\bf y}\in\mbox {\Bb R}^k$,
\[
   {\bf C}{\bf y} =  {\bf A} {\bf Q}@{\bf y} = {\bf V} {\bf S} {\bf W}^T {\bf Q}@{\bf y} + {\bf E}{\bf Q}@ {\bf y},
\]
where ${\bf E} = {\bf A} - {\bf V}{\bf S}{\bf W}^T$.\ \ 
Since ${\bf V}^T {\bf E} = {\bf 0},$
\[
   \| {\bf C}{\bf y}\|^2  = \|  {\bf A} {\bf Q}@{\bf y}\|^2  =\|  {\bf V} {\bf S} {\bf W}^T {\bf Q}{\bf y}\|^2  + \| {\bf E}{\bf Q}@{\bf y} \|^2.
\]
Since  $\| {\bf W}^T {\bf Q}{\bf y}\| \ge \|{\bf y}\|/\|({\bf W}^T{\bf Q})^{-1}\| = 1/\eta_q$,
\[
   \| {\bf C}{\bf y}\|  \ge \|  {\bf V} {\bf S} {\bf W}^T {\bf Q}{\bf y}\| \ge \sigma_k / \eta_q > 0.
\]
Thus ${\bf C}$ must be full rank.   A similar argument shows ${\bf R}$ to be full rank as well.
\medskip

The examples in Section~\ref{sec:examples} illustrate that $\eta_p$ and $\eta_q$ are often quite modest for the DEIM-CUR approach, e.g., $ {\cal O} (100)$.\ \  However, worst-case  bounds permit significant growth in $k$ that is generally not observed in practice.  
We begin by stating a bound on this growth developed by
Chaturantabut and Sorensen~\cite[Lemma~3.2]{CS10}.

\begin{lemma} 
For the DEIM selection scheme derived above,
\[\eta_p \le {(1+\sqrt{2m})^{k-1}\over \|{\bf v}_1\|_\infty}, \qquad
   \eta_q \le {(1+\sqrt{2n})^{k-1} \over \|{\bf w}_1\|_\infty},\]
where ${\bf v}_1$ and ${\bf w}_1$ denote the first columns of ${\bf V}$ and ${\bf W}$.
\end{lemma}

Motivated by recent work by Drma\v{c} and Gugercin~\cite{DG} on a modified DEIM-like algorithm for model reduction, we can improve this bound considerably.
\begin{lemma} \label{lemma:new_eta_bound}
For the DEIM selection scheme derived above,
\[
   \eta_p < \sqrt{mk\over 3} \ {2^k}, \qquad
   \eta_q < \sqrt{nk\over 3} \ {2^k}.
\]
\end{lemma}
\noindent{\it Proof:} 
We shall prove the result for $\eta_p$; the result for $\eta_q$ follows similarly.
As usual, let ${\bf V} \in \mbox {\Bb R}^{m \times k}$ have orthonormal columns, 
and let ${\bf p} = {\rm DEIM}({\bf V})$ denote the row index vector derived 
from the DEIM selection scheme described above.  
Let ${\bf P} = {\bf I}(\,:\,,{\bf p})$ so that ${\bf P}^T {\bf V} = {\bf V}({\bf p},\,:\,)$.  

Without loss of generality, assume the DEIM index selection gives 
${\bf p} = [1, 2, \ldots, k]^T$.  (Otherwise, introduce a permutation matrix to the argument that follows.)
As described in section~\ref{sec:deim}, the DEIM index selection 
is precisely the index selection of LU decomposition with partial pivoting,
so one can write
\[
    {\bf V} = {\bf L}\kern1pt {\bf T},
\]
where the nonsingular matrix ${\bf T}\in\mbox {\Bb R}^{k\times k}$ is upper triangular 
and ${\bf L}\in\mbox {\Bb R}^{m\times k}$ is unit lower triangular with $ | {\bf L}(i,j) | \le 1$, 
${\bf L}(j,j) = 1, \ 1 \le j \le k$ and ${\bf L}(i,j) = 0, \ j > i $.

Let ${\bf L}_1 \equiv {\bf L}(1:k,1:k)$.  Then ${\bf V}({\bf p},\,:\,) = {\bf L}_1 {\bf T}$ and thus
\[
 \eta_p \equiv \| ({\bf P}^T {\bf V})^{-1} \| = \|({\bf L}_1 {\bf T})^{-1} \| \le \| {\bf T}^{-1}\| \| {\bf L}_1^{-1}\|.
\]
(The linear independence of the columns of ${\bf V}$ ensure that ${\bf L}_1$ and ${\bf T}$ are invertible.)
Upper bounds for $\|{\bf T}^{-1}\|$ and $\| {\bf L}_1^{-1}\|$ will give an upper bound for $\eta_p$.

To bound $\| {\bf T}^{-1}\|$, let ${\bf y}\in\mbox {\Bb R}^k$ be a unit vector such that 
$\| {\bf T}^{-1} {\bf y} \| = \| {\bf T}^{-1}\|$.
Then
\[
   \| {\bf T}^{-1} \| = \| {\bf T}^{-1} {\bf y} \| = \| {\bf V} {\bf T}^{-1} {\bf y} \| = \| {\bf L} {\bf y} \|.
\]
Now
\[
   \| {\bf L} {\bf y} \| \le \sqrt{m}\, \|{\bf L}{\bf y}\|_\infty = \sqrt{m}\, |{\bf e}_j^T {\bf L} {\bf y} |,
\]
for some index $j\in\{1,\ldots, m\}$.\ \ 
By the Cauchy--Schwarz inequality and the bound $|{\bf L}(i,j)|\le 1$, 
\[
  |{\bf e}_j^T {\bf L} {\bf y} | \le \|{\bf L}^T {\bf e}_j\| \| {\bf y} \| \le \sqrt{k}\cdot 1,
\]
and so it follows that $\| {\bf T}^{-1} \| \le \sqrt{mk}$.

The inverse of ${\bf L}_1$ can be bounded using forward substitution.
Let ${\bf L}_1 {\bf z} = {\bf y}$, where $\|{\bf y}\| = 1$ and $\| {\bf z} \| =  \| {\bf L}_1^{-1}\|$.
Forward substitution provides
\begin{eqnarray*}
   \zeta_1 &=& \gamma_1 \\
   \zeta_i &=& \gamma_i - \sum_{j=1}^{i-1} \lambda_{ij} \zeta_j, \rlap{\qquad $i = 2,\ldots, k$,}\\
\end{eqnarray*}
where $\zeta_i = {\bf z}(i)$, $\gamma_i = {\bf y}(i)$ and $\lambda_{ij} = {\bf L}(i,j)$.
We now use induction to prove
\[
    | \zeta_i | \le 2^{i-1}, \rlap{\qquad $1 \le i \le k$.}
\]
First note that $\zeta_1 = \gamma_1$, so $|\zeta_1| \le |\gamma_1| \le 1 = 2^0$
to establish the base case.  Assume for some $i \ge 1$ that
\[
    | \zeta_j | \le 2^{j-1}, \rlap{\qquad $1 \le j \le i$.}
\]
Then
\begin{eqnarray*}
   | \zeta_{i+1} | = \bigg|  \gamma_{i+1} - \sum_{j=1}^{i} \lambda_{ij} \zeta_j \bigg|
                   &\le& |\gamma_{i+1}| + \sum_{j=1}^{i} |\lambda_{ij}| |\zeta_j | \\
                   &\le&  1 + \sum_{j=1}^{i} 1 \cdot 2^{j-1}
                   =  1 + \sum_{j=0}^{i-1}  2^j 
                   = 1 + (2^i - 1) = 2^i
\end{eqnarray*}
to complete the induction.
Now since $\|{\bf z}\| = \|{\bf L}_1^{-1}\|$,
\[
   \| {\bf L}_1^{-1} \|^2 = {\bf z}^T {\bf z} = \sum_{i=1}^k |\zeta_i|^2 \le \sum_{i=0}^{k-1} 4^i = (4^k - 1)/3.
\]
Thus $ \| {\bf L}_1^{-1} \| < 2^k/\sqrt{3}$, which, together with the bound on the inverse of ${\bf T}$, provides
the final result for $m>k$:
\[
     \eta_p \equiv \| ({\bf P}^T {\bf V})^{-1} \| < \sqrt{mk\over 3}\, 2^k.
\]
If $m=k$, then $\eta_p=1$, and the result holds trivially.
\quad  { { \hfill\rule{3mm}{3mm} } } 

\medskip
Note that this proof only relies on the orthonormality of the columns of ${\bf V}$ and ${\bf W}$,
and hence it applies when the DEIM selection scheme is applied to approximate 
singular vectors, as in CUR error bound in~(\ref{eq:approxbnd}).

Lemma~\ref{lemma:new_eta_bound}  
was inspired by the proof technique developed by 
Drma\v{c} and Gugercin~\cite{DG} to bound $\|({\bf P}^T {\bf V})^{-1} \|$,
when ${\bf P}$ is selected by applying a pivoted rank-revealing 
QR factorization scheme to ${\bf V}$.\ \   
Note that this new bound is on the same order of magnitude as the 
Drma\v{c}--Gugercin scheme.  In practice, their scheme
seems to give slightly smaller growth that is more consistent over a wide range of examples.
Neither scheme experienced exponential growth over very extensive testing.
For the DEIM approach, this absence of exponential growth is closely related
to decades of experience with Gaussian elimination with partial pivoting.   
Element growth in ${\bf T}$ is bounded by a factor of $2^{k-1}$ 
(for a $k \times k$ matrix), and there is an example that achieves this 
growth.   Nevertheless, this algorithm is almost exclusively used to solve
linear systems because such growth is never experienced.%
\footnote{See, for example, the extensive numerical tests 
involving random matrices described 
in~\cite[lecture~22]{TB97} and~\cite{TS90}.
Interestingly, in the experiments of Trefethen and Schreiber~\cite{TS90},
random matrices with orthonormal columns tend to have slightly larger 
growth factors than Gaussian matrices, though both cases are very 
far indeed from the exponential upper bound.}
Indeed, a similar near worst case example can be constructed for DEIM,
although this growth has not been observed in practice.

\medskip
\noindent{\bf A Growth Example:}
We now construct an orthonormal matrix ${\bf V}$ with the property
\begin{equation}\label{eq:growthex}
   {1\over \sqrt{8}}\ 2^{k}  <  \eta_p \equiv \|({\bf P}^T {\bf V})^{-1} \| < \sqrt{m k\over 3}\ 2^k
\end{equation}
where ${\bf P}^T {\bf V} = {\bf V}({\bf p},\,:\,)$ with ${\bf p} = {\rm DEIM}({\bf V})$.  
To construct ${\bf V}$, begin by defining
\[ {\bf L} := \left[\begin{array}{cccc}
1 & & & \\
-1 & 1 & & \\
\vdots & \ddots & \ddots & \\
-1 & \cdots & -1 & 1 \\
-1 & \cdots & -1 & -1 \\
\vdots &     &   &  \vdots \\
-1 & \cdots & -1 & -1
           \end{array}\right]\in \mbox {\Bb R}^{m\times k}.\]
Now construct ${\bf V}{\bf T}_1 \equiv {\bf L}$ as an economy-sized QR factorization of
${\bf L}$ (with no column pivoting).  Since the columns of ${\bf L}$ are
linearly independent by construction, ${\bf T}_1\in\mbox {\Bb R}^{k\times k}$ is invertible;
define ${\bf T} \equiv {\bf T}_1^{-1}$, so that ${\bf V} = {\bf L}\kern1pt{\bf T}$.
(Note that ${\bf T}$ plays the same role it does in the proof of Lemma~\ref{lemma:new_eta_bound}.)
If the DEIM procedure is applied to ${\bf V}$, then by construction
${\bf p} = [1, 2, \ldots, k]$  (in exact arithmetic):
during the DEIM procedure,  the relations
\[
     \mbox{\boldmath$\ell$}_j\kern1pt\tau_{jj} = {\bf v}_j -  {\bf V}_{j-1}^{} ( {\bf P}_{j-1}^T {\bf V}_{j-1}^{})^{-1} {\bf P}_{j-1}^T {\bf v}_j, \rlap{\qquad $j > 1$}
\]
hold, with $\mbox{\boldmath$\ell$}_j  = {\bf L}(\,:\,,j)$, $\tau_{jj} = {\bf T}(j,j)$, 
${\bf v}_j = {\bf V}(\,:\,,j)$ , ${\bf P}_{j-1} = {\bf I}(\,:\,,1:j-1)$
and ${\bf V}_{j-1} = {\bf V}(\,:\,,1:j-1)$.  Thus ${\bf p}(j) = j , \ j>1 $ and it is easily seen that ${\bf p}(1) = 1$.

Note that $ {\bf V} {\bf T}^{-1} = {\bf L} $ implies ${\bf T}^{-T} {\bf T}^{-1} = {\bf L}^T {\bf L}$ hence $\| {\bf T} \| = 1/\sigma_k$,
where $\sigma_k $ is the smallest singular value of ${\bf L}$.\ \ 
Let ${\bf y}$ be the corresponding right singular vector, so that
\[
   \sigma_k^2 = {\bf y}^T {\bf L}^T {\bf L} {\bf y}.
\]
We claim that $\sigma_k \ge \sqrt{2}$.\ \ 
To see this, write ${\bf L}$ in the form
\[
   {\bf L} = \left[ \begin{array}{c}
                     {\bf I}_k \\
                     {\bf 0}
                \end{array}
           \right]
         - \left[ \begin{array}{c}
                     {\bf L}_0 \\
                      {\bf E}
                \end{array}
           \right],
\]
where 
\[ {\bf I}_k = \left[\begin{array}{cccc}
                   1 \\ & 1 \\ & & \ddots \\ & & & 1
            \end{array}\right], 
    \quad
    {\bf L}_0 = \left[\begin{array}{cccc}
               0 \\ 1 & \ddots \\ \vdots & \ddots & \ddots \\
               1 & \cdots & 1 & 0\end{array}\right],
     \quad 
     {\bf E} = \left[\begin{array}{cccc}
               1 & 1 & \cdots & 1 \\ 
               1 & 1 & \cdots & 1 \\ 
               \vdots & \vdots & \ddots & \vdots \\
               1 & 1 & \cdots  & 1\end{array}\right]
            = {\bf f}\kern1pt {\bf e}^T,
\]
for ${\bf f} = [1,\ldots, 1]^T \in \mbox {\Bb R}^{m-k}$.\ \ 
\begin{eqnarray*} 
 {\bf L}^T {\bf L} &=& {\bf I}_k - {\bf L}_0 - {\bf L}_0^T + {\bf L}_0^T {\bf L}_0 + {\bf E}^T {\bf E} \\[.25em]
           &=& {\bf I}_k - ({\bf e} {\bf e}^T - {\bf I}_k) + {\bf L}_0^T{\bf L}_0 + (m-k){\bf e}\bfe^T \\[.25em]
           &=& 2\kern1pt {\bf I}_k + {\bf L}_0^T {\bf L}_0 + (m - k -1) {\bf e} {\bf e}^T \\[.25em]
           &=& 2\kern1pt {\bf I}_k +  {\bf M},
\end{eqnarray*}
where ${\bf M}:= {\bf L}_0^T {\bf L}_0 + (m - k -1) {\bf e} {\bf e}^T$ 
is symmetric and positive semidefinite whenever $m>k$.
Thus
\[
   \sigma_k^2 = {\bf y}^T {\bf L}^T {\bf L} {\bf y} = {\bf y}^T (2{\bf I}_k +  {\bf M}) {\bf y} \ge 2
\]
and hence it follows that
\[
   \| {\bf T} \| \le 1/\sqrt{2}.
\]
This implies
\[
   \|({\bf P}^T {\bf V})^{-1} \| = \| {\bf T}^{-1} {\bf L}_1^{-1} \| \ge \frac{\| {\bf L}_1^{-1} \|}{\| {\bf T} \|} \ge
     \sqrt{2}\ \|{\bf L}_1^{-1} \|.
\]
To complete the lower bound, we must analyze $\|{\bf L}_1^{-1}\|$.  Forward substitution gives ${\bf L}_1^{-1} {\bf e}_1 = [1, 1, 2, 4, \ldots , 2^{k-2}]^T$ and thus
\[
    \|{\bf L}_1^{-1} \| > \|{\bf L}_1^{-1} {\bf e}_1\| = \sqrt{1 + (4^{k-1} - 1)/3} \ > 2^{k-2}.
\]
We arrive at the lower bound
\[
   \eta_p \equiv \|({\bf P}^T {\bf V})^{-1} \| 
     \ge  \sqrt{2}\ \|{\bf L}_1^{-1} \| \  > \sqrt{2}\cdot 2^{k-2},
\]
and thus for this choice of ${\bf V}$, the DEIM error constant satisfies
\[
  {1\over\sqrt{8}}\ 2^k  < \eta_p < \sqrt{{m k\over 3}} \ 2^k.
\]

This example is interesting because it relies on the behavior of the
classic example for growth in LU decomposition~\cite[lecture~22]{TB97}.
However, in this case the pathological growth is caused by ${\bf L}$ 
and not by ${\bf T}$.

2015-09-22T02:00:52.000
1407
1407.5516
en
https://arxiv.org/abs/1407.5516
We derive a CUR matrix factorization based on the Discrete Empirical Interpolation Method (DEIM). For a given matrix $A$, such a factorization provides a low rank approximate decomposition of the form $A \approx C U R$, where $C$ and $R$ are subsets of the columns and rows of $A$, and $U$ is constructed to make $CUR$ a good approximation. Given a low-rank singular value decomposition $A \approx V S W^T$, the DEIM procedure uses $V$ and $W$ to select the columns and rows of $A$ that form $C$ and $R$. Through an error analysis applicable to a general class of CUR factorizations, we show that the accuracy tracks the optimal approximation error within a factor that depends on the conditioning of submatrices of $V$ and $W$. For large-scale problems, $V$ and $W$ can be approximated using an incremental QR algorithm that makes one pass through $A$. Numerical examples illustrate the favorable performance of the DEIM-CUR method, compared to CUR approximations based on leverage scores.
Numerical Analysis (math.NA)
A DEIM Induced CUR Factorization
Qwen/Qwen-72B
1. YES
2. YES