# A DEIM-CUR factorization with iterative SVDs11footnotemark: 1

Perfect Y. Gidisu<EMAIL_ADDRESS>Michiel E. Hochstenbach
<EMAIL_ADDRESS>
###### Abstract

A CUR factorization is often utilized as a substitute for the singular value
decomposition (SVD), especially when a concrete interpretation of the singular
vectors is challenging. Moreover, if the original data matrix possesses
properties like nonnegativity and sparsity, a CUR decomposition can better
preserve them compared to the SVD. An essential aspect of this approach is the
methodology used for selecting a subset of columns and rows from the original
matrix. This study investigates the effectiveness of _one-round sampling_ and
iterative subselection techniques and introduces new iterative subselection
strategies based on iterative SVDs. One provably appropriate technique for
index selection in constructing a CUR factorization is the discrete empirical
interpolation method (DEIM). Our contribution aims to improve the
approximation quality of the DEIM scheme by iteratively invoking it in several
rounds, in the sense that we select subsequent columns and rows based on the
previously selected ones. That is, we modify $A$ after each iteration by
removing the information that has been captured by the previously selected
columns and rows. We also discuss how iterative procedures for computing a few
singular vectors of large data matrices can be integrated with the new
iterative subselection strategies. We present the results of numerical
experiments, providing a comparison of one-round sampling and iterative
subselection techniques, and demonstrating the improved approximation quality
associated with using the latter.

###### keywords:

DEIM, CUR decomposition , iterative SVD, low-rank approximation, adaptive
sampling, iterative subselection

††journal: xxxxx

##  1 Introduction

In data analysis applications and machine learning, the data set is often
represented by a matrix $A\in\mathbb{R}^{m\times n}$, where $m$ is the number
of data points and $n$ is the number of features or variables. Data matrices
are usually large and in many cases, a key step in the analysis is to
approximate the data using a few features and/or a few data points so that one
can easily manipulate, understand, and interpret the data. The optimal
approximation is obtained by the truncated singular value decomposition
(TSVD). On the other hand, this approximation problem can also be reduced to
identifying a good subset of columns and rows in the data matrix, a CUR
factorization. A CUR decomposition is an alternative solution to the SVD,
motivated by the fact that in several applications, it can be challenging to
have a concrete interpretation of the singular vectors. Additionally, the
singular vectors fail to preserve properties such as nonnegativity, and
sparsity, if the original data matrix has these. Theoretical computer science
and numerical linear algebra communities have extensively studied column
subset selection (index selection) algorithms. In the former, researchers have
primarily focused on selecting good columns using randomized algorithms with
provable error bounds [4, 8, 10, 12, 17]. On the other hand, numerical linear
algebra researchers have primarily focused on deterministic algorithms, which
exploit the SVD or rank-revealing QR factorizations that select columns by
pivoting rules [3, 6, 16, 23]. This paper focuses on the deterministic
algorithms that exploit the SVD for index selection; in particular, the
discrete empirical interpolation method (DEIM) [2, 7, 11, 20].

Throughout the paper, we denote the spectral (2-norm) and the Frobenius norm
by $\lVert\cdot\rVert$ and $\lVert\cdot\rVert_{F}$, respectively, and the
infinity-norm by $\lVert\cdot\rVert_{\infty}$. The notations $C^{+}$ and
$C^{\top}$ denote the Moore–Penrose pseudoinverse and the transpose of $C$,
respectively. We use MATLAB notation to index vectors and matrices; thus,
$A(:,\mathbf{p})$ denotes the $k$ columns of $A$ whose corresponding indices
are in vector $\mathbf{p}\in\mathbb{N}_{+}^{k}$.

The DEIM algorithm [2, 7] is a technique used to select some important column
and row indices from an $m\times n$ data matrix $A$ [20], where without loss
of generality $m\geq n$. We wish to select $k\ll n$ relevant row and column
indices. The first step in the DEIM procedure is to compute an (reduced) SVD,
$A=U\Sigma V^{\top}$. The associated cost when a direct method is used is
$\mathcal{O}(mn^{2})$, independent of the value of $k$. The paper also
considers iterative methods to approximate the leading $k$ singular vectors.
Having the left and right singular vectors contained in $U$ and $V$,
respectively, we are interested in selecting distinct row indices
$s_{1},\dots,s_{k}$ from the set $\\{1,\dots,m\\}$ and column indices
$p_{1},\dots,p_{k}$ from the set $\\{1,\dots,n\\}$. The result of the method
may also be represented by an $m\times k$ row selection matrix $S$ and an
$n\times k$ column selection matrix $P$, whose columns are the standard basis
vectors indexed by the selected indices. The corresponding CUR factorization
is of the form (In line with [22], we will use the letter $M$ rather than $U$
for the middle matrix because $U$ is used to denote the matrix containing left
singular vectors)

$\begin{array}[]{ccccc}A&\approx&C&M&R~{},\\\\[-1.42262pt] {\scriptstyle
m\times n}&&{\scriptstyle m\times k}&{\scriptstyle k\times k}&{\scriptstyle
k\times n{}}\end{array}$ (1)

where matrices $C=AP$ and $R=S^{\top}A$ (both of full rank) are a subset of
the columns and rows of A, respectively, and the middle matrix $M$ of full
rank is computed such that the decomposition is as close to $A$ as possible.
Given $C$ and $R$, a standard procedure to determine $M$ is (see, e.g., [20,
Sec. 2], where also an alternative is presented) by two consecutive least
squares problems:

1: | Solve the least squares problem $CX\approx A$ for $X\in\mathbb{R}^{k\times n}$
---|---
| with solution $X=(C^{\top}C)^{-1}C^{\top}A$.
2: | Solve the least squares problem $R^{\top}M^{\top}\approx X^{\top}$ for $M\in\mathbb{R}^{k\times k}$
| with solution $M=XR^{\top}(RR^{\top})^{-1}$.

Both steps are optimal with respect to the spectral and Frobenius norm. It is
important to note that the solution in the spectral norm may not be unique. In
many applications, one cares primarily about key columns or rows of $A$,
rather than an explicit $A\approx CUR$ factorization. Thus, an interpolative
decomposition of the form $A=C\widehat{M}$ or $A=\widehat{M}R$ [23].

We will now describe how the DEIM scheme selects the row indices, and the
selection of the column indices follows similarly. The DEIM algorithm starts
from the dominant left singular vector $\mathbf{u}_{1}$, and the first index
$s_{1}$ corresponds to the largest magnitude entry in $\mathbf{u}_{1}$, that
is $|\mathbf{u}_{1}(s_{1})|=\lVert\mathbf{u}_{1}\rVert_{\infty}$. Given that
$I$ is the identity matrix, let $\mathbf{s}=[s_{1}]$, $S_{1}=I(:,s_{1})$,
$U_{1}=[\mathbf{u}_{1}]$, and define an oblique projection operator as
$\mathbb{S}=\mathbf{u}_{1}(S_{1}^{\top}\mathbf{u}_{1})^{-1}S_{1}^{\top}$.

Suppose we have $j-1$ indices, so that

$\mathbf{s}_{j-1}=\begin{bmatrix}s_{1}\\\ \vdots\\\ s_{j-1}\end{bmatrix},\quad
S_{j-1}=I(:,\mathbf{s}_{j-1}),\quad
U_{j-1}=[\mathbf{u}_{1},\dots,\mathbf{u}_{j-1}],$

and

$\mathbb{S}_{j-1}=U_{j-1}(S_{j-1}^{\top}U_{j-1})^{-1}S_{j-1}^{\top}.$

Compute the residual vector
$\mathbf{r}_{j}=\mathbf{u}_{j}-\mathbb{S}_{j-1}\mathbf{u}_{j}$, and select the
next index $s_{j}$ such that
$|\mathbf{r}_{j}(s_{j})|=\lVert\mathbf{r}_{j}\rVert_{\infty}$. It is worth
noting that, using this oblique projection operator $\mathbb{S}_{j-1}$ on
$\mathbf{u}_{j}$ ensures that the $s_{j-1}$ entry in $\mathbf{r}_{j}$ is 0,
which guarantees nonrepeating indices. We also note that for the projector
$\mathbb{S}_{j-1}$ to exist at the $j$th step, $S_{j-1}^{\top}U_{j-1}$ must be
nonsingular. The linear independence of the columns of $U$ guarantees this.
The DEIM index selection procedure is summarized in Algorithm 1 222Note that
the backslash operator used in the algorithms is a Matlab type notation for
solving linear systems and least-squares problems.. A variant of the DEIM
scheme proposed by [11] called QDEIM involves computing a column-pivoted QR
factorization on the transpose of the singular vectors to obtain the column
and row indices, which correspond to the indices of the first $k$ pivoted
columns.

Data: $U\in\mathbb{R}^{m\times k}$ with $k\leq m$ (of full rank)

Result: Indices $\mathbf{s}\in\mathbb{N}_{+}^{k}$ with non-repeating entries

1

2$\mathbf{s}(1)$ = $\operatorname*{\arg\\!\max}_{1\leq i\leq
m}~{}|(U(:,\,1))_{i}|$

3for  _ $j=2,\dots,k$_ do

4 $U(:,\,j)=U(:,\,j)-U(:,\,1:j-1)\cdot(U(\mathbf{s},\,1:j-1)\ \backslash\
U(\mathbf{s},\,j))$

5 $\mathbf{s}(j)$ = $\operatorname*{\arg\\!\max}_{1\leq i\leq
m}~{}|(U(:,\,j))_{i}|$

6 end for

Algorithm 1 DEIM index selection scheme [20]

In this paper, we explore iterative subselection variants of the DEIM scheme.
Our contribution aims to improve the approximation quality of the DEIM scheme
by iteratively invoking it, in the sense that we select subsequent columns and
rows based on the previously selected ones. That is, we modify $A$ after each
iteration by removing the information that has been captured by the previously
selected columns and rows. We show this by adapting an existing iterative
subselection technique for the DEIM index selection scheme and also propose a
new strategy. We also discuss how iterative procedures for computing a few
singular vectors of large data matrices can be used with our proposed methods.
To the best of our knowledge, this is the first deterministic DEIM-type
algorithm for large-scale data sets.

##  2 Adaptive sampling for column subset selection problem

The iterative subselection strategies proposed in this work are related to the
so-called _adaptive sampling_ for column subset selection. In this section, we
provide an overview of the adaptive sampling technique proposed by Deshpande
et al. [8]. The authors introduce a probabilistic method that iteratively
selects a subset of columns in multiple rounds to construct a rank-$k$
approximation of a matrix. This approach has been demonstrated to provide
improved accuracy and flexibility compared to _one-round sampling_ methods.
One-round sampling methods refer to selection schemes that obtain all $k$
columns in a single round.

Data: $A\in\mathbb{R}^{m\times n}$, target rank $k$, # rounds $t$, columns per
round $c$

1

Result: $C\in\mathbb{R}^{m\times tc}$

2

3$\mathbf{p}=[\ ];\quad E_{0}=A$

4for  _$i=1,\dots,t$_ do

5

6 for  _$j=1,\dots,n$_ do

7

8 if $j\in\mathbf{p}$ then $\text{pr}_{j}^{(i)}=0$ (sample without
replacement)

9 else $\text{pr}_{i}^{(j)}={\lVert E_{i-1}^{(j)}\rVert^{2}}\,/\,{\lVert
E_{i-1}\rVert_{F}^{2}}$

10 end for

11

12 $\mathbf{p}_{i}$ = set of $c$ indices sampled according to $\text{pr}_{i}$

13 $\mathbf{p}=[\mathbf{p}\ \ \mathbf{p}_{i}]$

14 $C=A(\,:,\mathbf{p})$; $E_{i}=A-CC^{+}A$

15 end for

16

Algorithm 2 Adaptive sampling for column subset selection [8]

The adaptive sampling method of [8] as summarized in Algorithm 2 involves
alternating between two steps in each round: selecting a subset of columns and
updating the probability distribution over all columns. The selection of
columns in each round is influenced by the columns picked in previous rounds.
Suppose we aim to select $k$ subsets of columns from matrix $A$. The process
begins with an initial probability distribution and randomly selects $c<k$
columns to form a matrix $C$. The selection of columns is based on the norms
of the columns, as described in [8, 12]. Each column $j$ is chosen with a
probability $\text{pr}_{i}^{(j)}={\lVert E_{i-1}^{(j)}\rVert^{2}}\,/\,{\lVert
E_{i-1}\rVert_{F}^{2}}$ (as in Algorithm 2 of Algorithm 2). After selecting
$c$ columns, the probabilities are updated based on the chosen columns, and
$c$ new columns are sampled and added to the matrix $C$. This iterative
process continues until all $k$ columns are selected. Note that assigning zero
probability to previously selected indices, as described in Algorithm 2,
represents a sampling without replacement strategy.

The authors present a detailed explanation and theoretical analysis of this
adaptive sampling technique, emphasizing its advantages and diverse
applications [8]. The algorithm improves the accuracy of a CUR decomposition
compared to one-round sampling methods as demonstrated in [8, 9, 19, 24].
Moreover, it allows for flexibility by accommodating different criteria for
selecting column and row subsets based on specific problem requirements, which
we will also discuss in Section 3. In their approach (Algorithm 2), a constant
number of columns is selected per iteration, and the residual is computed as
$E=A-CC^{+}A$.

Nevertheless, one notable drawback of this adaptive sampling technique is its
probabilistic nature, introducing randomness into the index selection process.
This randomness may sacrifice predictability to some extent. Applying the
algorithm multiple times to the same data may lead to inconsistent
factorizations, yielding varying CUR decompositions across runs. This lack of
predictability may hamper the reproducibility of the decomposition. Moreover,
the unpredictable nature may impact downstream tasks reliant on the CUR
decomposition, potentially leading to inconsistent outcomes. As a result, the
practical use of a CUR decomposition by this type of probabilistic algorithm
may be compromised.

To address this limitation, in the next section, we introduce a modified
approach based on the DEIM scheme to implement the iterative subselection
algorithm proposed by Deshpande et al. [8] and also propose a new iterative
subselection strategy. By incorporating the DEIM scheme, our method provides a
deterministic technique for iterative subselection of column indices, in
contrast to the original methods that utilized a probabilistic approach [4, 8,
9, 19, 24]. This deterministic nature enhances the reproducibility of the CUR
decomposition, making it more suitable for practical applications.

##  3 Small-scale DEIM type CUR with iterative SVDs

In this section, we introduce new index-picking schemes for constructing a CUR
decomposition. The standard DEIM scheme computes an SVD of $A$ once, after
which the indices are picked iteratively “locally optimal”. The new methods
that we present now compute an SVD in every iteration. The algorithms
adaptively select columns and rows of $A$ in several rounds. In each
iteration, we modify $A$ by removing the information that has been captured by
the previously selected columns and rows. The time complexities of the various
proposed methods after $t$ rounds are summarized in Table 1. This includes the
computational time for computing an SVD and an updated $A$ (residual) matrix
in every round.

Table 1: Summary of the dominant work of the different algorithms after $t$ rounds. The time complexity column excludes the computational cost of the DEIM scheme as it is approximately the same for all algorithms. Method | Matrix | SVD | Time
---|---|---|---
|  |  | svd |  $X$ or $M$ |  Residual ($E$)
CADP-CX (Algorithm 3) | Small | Full | $\mathcal{O}(tmn^{2})$ | $\mathcal{O}(tmnk)$ | $\mathcal{O}(tmnk)$
DADP-CX (Algorithm 4)
DADP-CUR (Algorithm 5)
Large: DADP-CX (Algorithm 7) | Large | Few | $\mathcal{O}(mn\cdot\text{nr}_{\text{in}})$ | $\mathcal{O}(mnk)$ | –

The complexity estimates and the other descriptions are similar for the first
three algorithms. It is important to highlight that when constructing a CUR
factorization using Algorithms 3 and 4, one needs to perform almost twice the
number of SVDs, $X$, and $E$, compared to what is required by Algorithm 5. For
the large-scale algorithm, estimating the precise time complexity of computing
a low-rank SVD using an iterative method (refer to Algorithm 6), as done in
Algorithm 7, can be challenging. The iterative approach involves a series of
matrix-vector multiplications and orthogonalization steps. The computational
cost of the Lanczos algorithm is typically dominated by matrix-vector
multiplications, which is $\mathcal{O}(mn)$ times the total number of inner
iterations (denoted by $\text{nr}_{\text{in}}$ in the table) needed. It is
worth noting that we do not explicitly compute the residual in Algorithm 7 as
done in Algorithms 3, 4, and 5.

###  3.1 A DEIM based adaptive sampling for column subset selection

We present a deterministic variant of the iterative subselection scheme
discussed in Section 2. The newly proposed algorithm (CADP-CX) builds upon the
original adaptive sampling algorithm [8] by leveraging the benefits of the
DEIM technique. The procedure is summarized in Algorithm 3. The method
involves iteratively selecting a constant number of column indices, denoted as
$c$, from $A$ in multiple rounds. We start by computing the leading $c<k$
singular vectors of $A$. Next, we apply the DEIM scheme (Algorithm 1) to these
singular vectors, resulting in the first set of $c$ indices. We then update
$A$ by computing the residual matrix $E$ using the interpolative decomposition
(as described in Algorithm 3 of Algorithm 3). Next, we compute the leading $c$
singular vectors of $E$ and apply the DEIM procedure again to obtain the next
set of $c$ indices. This process is repeated until we have selected all $k$
required indices.

Data: $A\in\mathbb{R}^{m\times n}$, target rank $k$, columns per round $c$
(with $c\,|\,k$)

1

Result: $C\in\mathbb{R}^{m\times k}$

2

3$\mathbf{p}=[\ ];\quad E_{0}=A$

4for  _$i=1,\dots,k/c$_ do

5

6 Compute $[\sim,\sim,V]$ = svd($E_{i-1}$)

7 $V(\mathbf{p},:\ )=0$

8 $\mathbf{p}_{i}={\sf deim}(V(:,1:c))$ (Iteratively pick $c$ indices)

9 $\mathbf{p}=[\mathbf{p}\ \ \mathbf{p}_{i}]$

10 $C=A(\,:,\mathbf{p})$; $X=C\backslash A$; $E_{i}=A-CX$

11 end for

12

Algorithm 3 DEIM based adaptive sampling for column subset selection

One consequence of utilizing the DEIM scheme here is that each column of $A$
has a chance of being selected in subsequent iterations, even if it was
selected in a previous iteration. This could result in the selection of
previously chosen columns in subsequent iterations. Algorithm 3 of Algorithm 3
is a possible sample without-replacement strategy that can alleviate this
problem. Since the DEIM procedure selects the index corresponding to the entry
of the largest magnitude in a given vector, when these indices are zeroed out
after being chosen, it guarantees that they will not be selected again.

With regards to the memory and computational complexity, computing the
residual in LABEL:{line:res} involves a full iteration over the matrix, which
has a space complexity of $\mathcal{O}(mn)$. Given that we use the DEIM
procedure and select $c$ columns per iteration, in terms of computational
complexity, a full SVD requires $\mathcal{O}(mn^{2})$, one run of the DEIM
algorithm requires $\mathcal{O}(mc^{2})$, and computing the residual in
Algorithm 3 costs $\mathcal{O}(mnk)$. The overall time complexity after $t$
rounds is $\mathcal{O}(tmn^{2}+tmc^{2}+tmnk)$.

###  3.2 A new iterative subselection method

In Algorithm 4, we introduce a new iterative subselection strategy (DADP-CX)
for a CUR factorization, which differs from the method employed in our
proposed Algorithm 3 and the adaptive sampling procedures in previous works
[4, 8, 19, 24]. In contrast to the previous strategy, which selects a fixed
number of columns or rows in each iteration, this new strategy dynamically
adjusts the selection schedule based on the decay of the singular values of
the data (the relative magnitudes of the singular values).

The motivation behind this new approach is to adapt the subselection process
according to the significance of the singular values. By considering the decay
pattern of the singular values, we can prioritize the selection of columns or
rows that contribute the most to the data’s overall structure and information.
The decay of singular values provides valuable information about the
significance of different components in the data. By leveraging this
information, the iterative subselection strategy can adapt to the specific
characteristics of the data and prioritize the selection of columns or rows
that contribute the most to its structure. This adaptability allows for a more
data-driven selection process.

1

Data: $A\in\mathbb{R}^{m\times n}$, desired rank $k$, threshold parameter
$\delta\in(0,1]$, upper limit $\ell$

2

Result: Low-rank CUR decomposition
$A_{k}=A(:,\mathbf{p})\,\cdot\,M\,\cdot\,A(\mathbf{s},:)$

3 Set $E=A$; $\mathbf{p}=[\ ]$; $\mathbf{s}=[\ ]$

4while  _length $(\mathbf{p})<k$_ do

5

6 Compute $[\sim,\Sigma,V]$ = svd($E$)

7 $V(\mathbf{p},:\ )=0$

8 Let $b$ be the last index $i\leq k-{\sf length}(\mathbf{p})$ with
$\sigma_{i}>\delta\,\sigma_{1}$

9 $c=\min(b,\ell)$; $\mathbf{p}_{c}={\sf deim}(V(:,1:c))$

10 $\mathbf{p}=[\mathbf{p}\ \ \mathbf{p}_{c}]$; $C=A(:,\mathbf{p})$;
$X=C\backslash A$

11 $E=A-CX$ (Update matrix)

12 end while

13

14Repeat steps 1–9 on $A^{\top}$ to find the row indices $\mathbf{s}$

15$M=A(:,\mathbf{p})\ \backslash\ (A\ /\ A(\mathbf{s},:))$

Algorithm 4 Singular value decay-based iterative DEIM with one-sided projected
residual

With a user-defined threshold $\delta\in(0,1]$, the small-scale version of our
method begins by computing the leading singular vectors corresponding to the
singular values greater than the threshold multiplied by the largest singular
value of $A$, i.e., all $\sigma_{i}>\delta\cdot\sigma_{1}$. Let $b$ denote the
number of singular values satisfying this condition. Additionally, we
introduce an extra parameter $\ell$ to establish an upper limit on the number
of indices per round, taking into account the number of singular values that
exceed the threshold. Consequently, we select the first $c=\min(b,\ell)$
column indices denoted by $\mathbf{p}_{c}$ by applying the DEIM scheme to the
leading $c$ right singular vectors ($V_{c}$).

Subsequently, we proceed to construct an interpolative decomposition using the
chosen column indices and compute the residual matrix $E$ by subtracting this
approximation from $A$, i.e., $E=A-CC^{+}A$. To determine the next set of
indices, we repeat the aforementioned process on $E$. Thus, we compute the
leading singular vectors of $E$ corresponding to singular values greater than
$\delta$ times the largest singular value of $E$ and repeat the procedure
mentioned earlier. As previously mentioned, since the DEIM scheme selects
indices corresponding to entries with the largest magnitude, we set the
entries in the right singular vectors that correspond to the previously
selected indices to zero. This prevents the selection of duplicate indices. We
continue this procedure until all $k$ indices are selected. We expect that the
multiple passes through $A$ would lead to a reduced approximation error. It is
worth mentioning that in Algorithm 4, there is no need to compute the initial
SVD of $A^{\top}$ since we can store the initial left singular vectors from
the SVD of $A$.

In addition to the new selection strategy described in Algorithm 4, we also
define an alternative way to compute the residual in the index selection
process, which is presented in Algorithm 5. The newly proposed iterative
subselection algorithms (Algorithms 3 and 4) and existing adaptive sampling
procedures such as those outlined in [4, 8, 9, 19, 24] define the residual as
the error incurred by projecting the matrix $A$ onto either the column space
of $C$ or the row space of $R$, i.e., $E=A-CC^{+}A$ or $E=A-AR^{+}R$,
respectively. In contrast, this new method (DADP-CUR) defines the residual as
the error incurred by simultaneously projecting $A$ onto both the column space
of $C$ and the row space of $R$. This means computing a CUR factorization at
each step using only the selected columns and rows.

Data: $A\in\mathbb{R}^{m\times n}$, desired rank $k$, threshold parameter
$\delta=(0,1]$, upper limit $\ell$

Result: Low-rank CUR decomposition
$A_{k}=A(:,\mathbf{p})\,\cdot\,M\,\cdot\,A(\mathbf{s},:)$

1

2Set $E=A$; $\mathbf{p}=[\ ]$; $\mathbf{s}=[\ ]$

3while  _length $(\mathbf{p})<k$_ do

4

5 $[U,\Sigma,V]$ = svd($E$)

6 Set $V(\mathbf{p},:\ )=0$, $U(\mathbf{s},:\ )=0$

7 Let $b$ be the last index $i\leq k-{\sf length}(\mathbf{p})$ with
$\sigma_{i}>\delta\,\sigma_{1}$

8 $c=\min(b,\ell)$; $\mathbf{p}_{c}={\sf deim}(V(:,1:c))$

9 $\mathbf{s}_{c}={\sf deim}(U(:,1:c))$

10 $\mathbf{p}=[\mathbf{p}\ \ \mathbf{p}_{c}]$, $\mathbf{s}=[\mathbf{s}\ \
\mathbf{s}_{c}]$

11 $M=A(:,\mathbf{p})\ \backslash\ (A\ /\ A(\mathbf{s},:))$

12 $E=A-A(:,\mathbf{p})\cdot M\cdot A(\mathbf{s},:)$ (Update matrix)

13 end while

Algorithm 5 Singular value decay-based iterative DEIM with two-sided projected
residual

By considering the simultaneous projection onto the column and row spaces, we
aim to use a residual that provides a more accurate representation of the
error in the CUR factorization. It takes into account the combined effect of
selecting specific columns and rows on capturing the underlying structure and
information in the data. This approach offers several potential advantages. It
allows for a more comprehensive assessment of the error in the index selection
process, considering the contributions from both the columns and rows.
Furthermore, it ensures that the residual accurately reflects the
approximation quality obtained by a CUR factorization using the selected
columns and rows. Additionally, it has the potential to reduce computational
costs compared to Algorithm 4, as the latter approach involves performing
nearly twice the number of SVDs required by Algorithm 5. Without showing
specifics, it is worth mentioning that Algorithm 3 can be adapted using the
newly defined residual.

Note that when $\delta=0$, both Algorithms 4 and 5 are equivalent to the DEIM-
type CUR factorization. In terms of time complexity, suppose we need $t$
iterations in Algorithms 4 and 5 to select all $k$ columns and rows. The cost
of solving $E$ is $\mathcal{O}(tmnk)$. The cost of an SVD and one run of the
DEIM scheme are $\mathcal{O}(mn^{2})$ and $\mathcal{O}(nc^{2})$, respectively,
where $c$ is the maximum number of columns selected per iteration. Therefore,
the overall cost of the algorithms is $\mathcal{O}(tmn^{2}+t(m+n)c^{2}+tmnk)$.
However, constructing $C$ and $R$ using Algorithm 4 requires two runs of it.
Thus, its cost is almost twice that of Algorithm 5.

For small matrices, these iterative subselection techniques may be worthwhile
as the costs are modest and the quality of the approximations may increase.
However, these schemes may be especially interesting for large matrices, for
which an SVD may be too expensive, and iterative methods are used to compute
the left and right singular vectors. We will study this situation in the next
section.

##  4 Large-scale DEIM type CUR with iterative SVDs

For large-scale matrices, taking an SVD every round in Algorithms 3, 4, and 5
will usually be prohibitively expensive. Indeed, even one (reduced) SVD will
be too costly, which means that the standard DEIM-type CUR decomposition is
generally not affordable. However, the proposed algorithm is suitable for
large-scale data, as approximating the largest singular vectors by iterative
(Krylov) methods is usually a relatively easy task. Additionally, here, we do
not explicitly compute the residual matrix as done in the proposed algorithms;
this is done implicitly in the computation of the approximate singular
vectors. Furthermore, instead of computing the full SVD as we do in Algorithms
3, 4, and 5, we now carry out:
1: M Approximate $\widehat{U}$ and $\widehat{V}$ of $E$.

This can efficiently be carried out by an implicitly restarted version of
Lanczos bidiagonalization (see, e.g., [1]). The idea is as follows. Let
$k<{\widehat{k}}$ be the minimal and maximal dimension of the subspaces. We
first carry out ${\widehat{k}}$ steps of Lanczos bidiagonalization summarized
by the matrix equations

$E\widehat{V}_{\widehat{k}}=\widehat{U}_{\widehat{k}}B_{\widehat{k}},\quad
E^{\top}\widehat{U}_{\widehat{k}}=\widehat{V}_{\widehat{k}}B_{\widehat{k}}^{\top}+\beta_{\widehat{k}}\widehat{\mathbf{v}}_{{\widehat{k}}+1}\mathbf{e}_{\widehat{k}}^{\top},$

where $B_{\widehat{k}}$ is bidiagonal. The singular values of
$B_{\widehat{k}}$ are approximations to those of $E$, and the singular vectors
lead to approximations to those of $E$. With the SVD
$B_{\widehat{k}}=W\widehat{\Sigma}Z^{\top}$, we get

$E(\widehat{V}_{\widehat{k}}Z)=(\widehat{U}_{\widehat{k}}W)\widehat{\Sigma},\quad
E^{\top}(\widehat{U}_{\widehat{k}}W)=(\widehat{V}_{\widehat{k}}Z)\widehat{\Sigma}+\beta_{\widehat{k}}\widehat{\mathbf{v}}_{{\widehat{k}}+1}(W^{\top}\mathbf{e}_{\widehat{k}})^{\top},$

For any upper triangular matrix $\widehat{\Sigma}$ an elegant implicit restart
procedure is possible; here $\widehat{\Sigma}$ is even diagonal. Order the
singular values in the desired way; in this case nonincreasingly. Partition
the transformed basis, redefining $\widehat{U}_{k}$ and $\widehat{V}_{k}$:

$(\widehat{U}_{\widehat{k}}W)=:[\widehat{U}_{k}\
\,\widehat{U}_{{\widehat{k}}-k}],\quad(\widehat{V}_{\widehat{k}}Z)=:[\widehat{V}_{k}\
\,\widehat{V}_{{\widehat{k}}-k}],\quad\widehat{\Sigma}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}\widehat{\Sigma}_{k}\\\
&\widehat{\Sigma}_{{\widehat{k}}-k}\end{array}\\!\\!\right]$}},$ (2)

and redefine $B_{k}=\widehat{\Sigma}_{k}$,
$\beta_{k+1}:=\beta_{{\widehat{k}}+1}$,
$\widehat{\mathbf{v}}_{k+1}:=\widehat{\mathbf{v}}_{{\widehat{k}}+1}$, and
$\mathbf{f}_{k}:=W^{\top}\mathbf{e}_{\widehat{k}}$. We can now conveniently
restart from the decomposition

$E\widehat{V}_{k}=\widehat{U}_{k}B_{k},\quad
E^{\top}\widehat{U}_{k}=\widehat{V}_{k}B_{k}^{\top}+\beta_{k}\widehat{\mathbf{v}}_{k+1}\mathbf{f}_{k}^{\top}.$

The pair $(\widehat{U}_{k},\widehat{V}_{k})$ may be viewed as a pair of
approximate invariant spaces with error $\|\mathbf{f}_{k}\|$. The spaces are
expanded with Lanczos bidiagonalization to dimension ${\widehat{k}}$, after
which the selection procedure is carried out again. This scheme is repeated
until the quantify $\|\mathbf{f}_{k}\|$ is sufficiently small. We summarize
the method in Algorithm 6. Note that MATLAB built-in function svds is a
different implementation of a related technique.

Data: $E\in\mathbb{R}^{m\times n}$, desired rank $k$, initial vector
$\mathbf{v}_{1}$, minimum and maximum dimension $k<{\widehat{k}}$, tolerance
tol

1

Result: Approximation to $k$ largest singular triplets
$(\sigma_{i},\mathbf{u}_{i},\mathbf{v}_{i})$, giving best low-rank
approximation
$E_{k}=\widehat{U}_{k}\widehat{\Sigma}_{k}\widehat{V}_{k}^{\top}$

2

3Generate $E\widehat{V}_{k}=\widehat{U}_{k}B_{k}$,
$E^{\top}\widehat{U}_{k}=\widehat{V}_{k}B_{k}^{\top}+\beta_{k+1}\widehat{\mathbf{v}}_{k+1}\mathbf{f}_{k}^{\top}$

4for  _$i=1,2,\dots$_ do

5 Expand to
$E\widehat{V}_{\widehat{k}}=\widehat{U}_{\widehat{k}}B_{\widehat{k}}$,
$E^{\top}\widehat{U}_{\widehat{k}}=\widehat{V}_{\widehat{k}}B_{\widehat{k}}^{\top}+\beta_{{\widehat{k}}+1}\widehat{\mathbf{v}}_{{\widehat{k}}+1}\mathbf{f}_{\widehat{k}}^{\top}$

6 Determine SVD $B_{\widehat{k}}=W\widehat{\Sigma}Z^{\top}$

7 Partition according to (2), restart with $\widehat{U}_{k}$,
$\widehat{V}_{k}$,

8 redefining $B_{k}:=\widehat{\Sigma}_{k}$,
$\beta_{k+1}:=\beta_{{\widehat{k}}+1}$,
$\widehat{\mathbf{v}}_{k+1}:=\widehat{\mathbf{v}}_{{\widehat{k}}+1}$,
$\mathbf{f}_{k}:=W^{\top}\mathbf{e}_{\widehat{k}}$

9 Stop if $\|\mathbf{f}_{k}\|\leq$ tol

10 end for

11

Algorithm 6 Implictly restarted Lanczos bidiagonalization [1]

Data: $A\in\mathbb{R}^{m\times n}$, desired rank $k$, threshold parameter
$\delta\in(0,1]$, upper limit $\ell$

Result: Low-rank CUR decomposition
$A_{k}=A(:,\mathbf{p})\,\cdot\,M\,\cdot\,A(\mathbf{s},:)$

1 Set $E=A$; $\mathbf{p}=[\ ]$;

2 while  _length $(\mathbf{p})<k$_ do

3 Compute $[\sim,\Sigma,V]$ = svds($E$) (Algorithm 6)

4 finding $\sigma_{b}$, the last index $i\leq k-{\sf length}(\mathbf{p})$ with
$\sigma_{i}>\delta\,\sigma_{1}$ or

5 at most $\sigma_{\ell}$. Let $c=\min(b,\ell)$

6 $V(\mathbf{p},:\ )=0$; $\mathbf{p}_{c}={\sf deim}(V(:,1:c))$

7 $\mathbf{p}=[\mathbf{p}\ \ \mathbf{p}_{c}]$; $C=A(:,\mathbf{p})$

8 Update an incremental QR decomposition $C=QT$

9 $E@(\mathbf{x})\ \mathbf{y}=A\mathbf{x};\
\mathbf{y}=\mathbf{y}-Q(Q^{\top}\mathbf{y});$ (Update matrix implicitly)

10 with transpose $E^{\top}@(\mathbf{x})\
\mathbf{y}=\mathbf{x}-Q(Q^{\top}\mathbf{x});\ \mathbf{y}=A^{\top}\mathbf{y}$

11

12 end while

13Repeat steps 1–9 on $A^{\top}$ to find the row indices $\mathbf{s}$

$M=A(:,\mathbf{p})\ \backslash\ (A\ /\ A(\mathbf{s},:))$

Algorithm 7 Large-scale: Singular value decay-based iterative DEIM with one-
sided projected residual

In Algorithm 7 we provide the large-scale version of Algorithm 4 by employing
Algorithm 6. Without showing details, it is worth noting that this can also be
adapted for Algorithms 3 and 5. It is important to note that the threshold
parameter $\delta$ and the upper limit $\ell$ on the number of indices to be
selected per round are incorporated within the implementation of Algorithm 6.
The efficiency of this method is because for the procedure in Algorithm 6, we
do not need the matrix $E$ in explicit form; only matrix-vector products (MVs)
with $E$ and $E^{\top}$ are necessary (see Algorithm 7). The routine of
Algorithm 6 also takes several MVs that depend on the distribution of the
singular value and the starting vector. The cost of computing the singular
values and vectors in Algorithm 7 depends on the total number of inner
iterations of Algorithm 6. The number of iterations required by the Lanczos
algorithm depends on the size of the matrix. In a matrix-vector product
$E\mathbf{x}$ for a vector $\mathbf{x}$, the component $A\mathbf{x}$ costs
$\mathcal{O}(mn)$ for a full matrix. In the case of a sparse matrix with $d$
nonzeros entries per row, the cost reduces to $\mathcal{O}(md)$. The
aggregated cost of Algorithm 7 is only $\mathcal{O}(mk^{2})$. In Algorithm 7,
the computation of $Q(Q^{\top}\mathbf{x})$ requires $\mathcal{O}(mk)$
operations. The cost of solving the least squares problem $M=C^{+}AR^{+}$
would be $\mathcal{O}(mnk)$, which is relatively expensive. Nevertheless, it
is important to highlight that this step is necessary for all CUR methods as
the final step.

As previously mentioned, it is not necessary to compute the initial SVD of
$A^{\top}$ in this case, as we can simply retain the initial left singular
vectors obtained from the SVD of $A$. The value of $\delta$ will typically
depend on the data set. A value close to $1$ may be favorable for the
approximation result but is more expensive since Algorithm 6 needs to be
carried out approximately $k$ times. However, having $\delta=1$ implies that
we select one index per iteration, and thus, we need just the first right and
left singular vectors of $E$ corresponding to the largest singular value. We
can reduce the computational cost by specifying an earlier convergence
criterion for finding the approximate leading right and left singular vectors.
We use Wedin’s theorem for this. The theorem bounds the distance between
subspaces and the proof is in (cf., e.g., [21, pp. 260–262]).

######  Theorem 4.1.

(Wedin’s Theorem) Given $E\in\mathbb{R}^{m\times n}$, let

$[U_{1}\ U_{2}\ U_{3}]^{\top}\ E\ [V_{1}\
V_{2}]=\left[\begin{array}[]{cc}\Sigma_{1}&0\\\ 0&\Sigma_{2}\\\
0&0\end{array}\right],$

be the SVD of $E$ (where the singular values are not necessarily
nonincreasing). The singular subspaces of interest are in the column spaces of
$U_{1}$ and $V_{1}$. Let the inexact/approximate singular subspaces be in the
column spaces of $\widehat{U}_{1}$ and $\widehat{V}_{1}$ in the decomposition

$[\widehat{U}_{1}\ \widehat{U}_{2}\ \widehat{U}_{3}]^{\top}\ \widehat{E}\
[\widehat{V}_{1}\
\widehat{V}_{2}]=\left[\begin{array}[]{cc}\widehat{\Sigma}_{1}&0\\\
0&\widehat{\Sigma}_{2}\\\ 0&0\end{array}\right].$

Now let $\Phi$ be the matrix of canonical angles between
$\operatorname{Range}(U_{1})$ and $\operatorname{Range}(\widehat{U}_{1})$, and
$\Theta$ be the matrix of canonical angles between
$\operatorname{Range}(V_{1})$ and $\operatorname{Range}(\widehat{V}_{1})$.
Given the residuals
$F_{1}=E\widehat{V}_{1}-\widehat{U}_{1}\widehat{\Sigma}_{1}$,
$F_{2}=E^{\top}\widehat{U}_{1}-\widehat{V}_{1}\widehat{\Sigma}_{1}$, suppose
that there is a number $\alpha>0$ such that

$\min|\sigma(\widehat{\Sigma}_{1})-\sigma(\Sigma_{2})|\geq\alpha\quad\text{and}\quad\sigma_{\min}(\widehat{\Sigma}_{1})\geq\alpha.$

Then

$\sqrt{\lVert\sin{\Phi}\rVert_{F}^{2}+\lVert\sin{\Theta}\rVert_{F}^{2}}\leq\frac{\sqrt{\lVert
F_{1}\rVert_{F}^{2}+\lVert F_{2}\rVert_{F}^{2}}}{\alpha}.$

Theorem 4.1 shows that the computed singular vectors extracted by the
projection method are optimal up to the factor in the right-hand side of the
above inequality. This implies that any change in the entries of the computed
singular vectors is bounded by this factor. Note that $\Sigma_{2}$ is unknown.
For our context, we use $\widehat{\Sigma}_{2}$ as an approximation to
$\Sigma_{2}$. Since we are only concerned with approximating the first leading
right singular vector $\widehat{\mathbf{v}}_{1}$, we approximate
$\alpha\approx\widehat{\sigma}_{1}-\widehat{\sigma}_{2}$. Let
$m_{1}(\widehat{\mathbf{v}}_{1})$ and $m_{2}(\widehat{\mathbf{v}}_{1})$ denote
the largest and second-largest entries in $\widehat{\mathbf{v}}_{1}$,
respectively, and let
$\mathbf{f}_{2}=E^{\top}\widehat{\mathbf{u}}_{1}-\widehat{\sigma}_{1}\widehat{\mathbf{v}}_{1}$
be the residual vector (associated with residual matrix $F_{2}$). The above
svds routine results in
$\mathbf{f}_{1}=E\widehat{\mathbf{v}}_{1}-\widehat{\sigma}_{1}\widehat{\mathbf{u}}_{1}=0$
(associated with residual matrix $F_{1}$). The DEIM algorithm selects the
index corresponding to the largest element in the magnitude of a vector.
Therefore, when $\delta=1$, one can set an early convergence criterion to find
the first singular vector that corresponds to the largest singular value,
using the following approximate bound:

$m_{1}(\widehat{\mathbf{v}}_{1})-m_{2}(\widehat{\mathbf{v}}_{1})\lesssim
2\,\frac{\lVert\mathbf{f}_{2}\rVert}{\widehat{\sigma}_{1}-\widehat{\sigma}_{2}}.$

##  5 Error Analysis

In this section, we provide a known theoretical error bound applicable to a
general class of CUR factorizations. Thus, this bound also holds for our
proposed methods. A detailed constructive proof is in [20], but we provide the
necessary details here for the reader’s convenience.

Let $P\in\mathbb{R}^{n\times k}$ and $S\in\mathbb{R}^{m\times k}$ be selection
matrices with some columns of the identity indexed by the indices selected by
employing the index selection techniques described in this paper.

######  Theorem 5.1.

[20, Thm. 4.1] Given $A\in\mathbb{R}^{m\times n}$ and a target rank $k$, let
$V\in\mathbb{R}^{n\times k}$ and $U\in\mathbb{R}^{m\times k}$ be the leading
$k$ right and left singular vectors of $A$. Suppose $C=AP$ and
$R=S^{\top}\\!A$ are of full rank, and $V^{\top}\\!P$ and $S^{\top}U$ are
nonsingular, then, with $M=C^{+}\\!AR^{+}$, a rank-$k$ CUR decomposition
constructed by the proposed techniques satisfies

$\lVert
A-CMR\rVert\leq(\eta_{\mathbf{s}}+\eta_{\mathbf{p}})\,\sigma_{k+1}\quad
with\quad\eta_{\mathbf{s}}<\sqrt{\tfrac{nk}{3}}\,2^{k}~{},\quad\eta_{\mathbf{p}}<\sqrt{\tfrac{mk}{3}}\,2^{k},$

where $\eta_{\mathbf{p}}=\lVert(V^{\top}\\!P)^{-1}\rVert$,
$\eta_{\mathbf{s}}=\lVert(S^{\top}U)^{-1}\rVert$.

Proof  Let $\mathbb{P}=P(V^{\top}P)^{-1}V^{\top}$ and
$\mathbb{S}=U(S^{\top}U)^{-1}S^{\top}$ be oblique projectors. Note that
$V^{\top}\mathbb{P}=V^{\top}$ and $\mathbb{S}U=U$, implying that
$V^{\top}(I-\mathbb{P})=0$ and $(I-\mathbb{S})U=0$. Using $M=C^{+}\\!AR^{+}$,
we have [18, 20]

$\displaystyle\lVert A-CMR\rVert$ $\displaystyle=\lVert
A-CC^{+}AR^{+}R\rVert=\lVert(I-CC^{+})A+CC^{+}A(I-R^{+}R)\rVert$
$\displaystyle\leq\lVert(I-CC^{+})A\rVert+\lVert CC^{+}\rVert\,\lVert
A(I-R^{+}R)\rVert.$

Leveraging the fact that $CC^{+}$ is an orthogonal projector, $\lVert
CC^{+}\rVert=1$ and [20, Lemmas 4.1 and 4.2]

$\displaystyle\lVert(I-CC^{+})A\rVert$ $\displaystyle\leq\lVert
A(I-\mathbb{P})\rVert=\lVert A(I-VV^{\top})(I-\mathbb{P})\rVert$
$\displaystyle\leq\lVert(V^{\top}\\!P)^{-1}\rVert\,\lVert
A(I-VV^{\top})\rVert,$ $\displaystyle\lVert A(I-R^{+}\\!R)\rVert$
$\displaystyle\leq\lVert(I-\mathbb{S})A\rVert=\lVert(I-\mathbb{S})(I-UU^{\top})A\rVert$
$\displaystyle\leq\lVert(S^{\top}U)^{-1}\rVert\,\lVert(I-UU^{\top})A\rVert,$

we have that

$\lVert A-CMR\rVert\leq\lVert(V^{\top}\\!P)^{-1}\rVert\,\lVert
A(I-VV^{\top})\rVert+\lVert(S^{\top}U)^{-1}\rVert\,\lVert(I-UU^{\top})A\rVert.$

Since $U$ and $V$ contain the leading $k$ left and right singular vectors,
respectively, $\lVert(I-UU^{\top})A\rVert=\lVert
A(I-VV^{\top})\rVert=\sigma_{k+1}$. Hence

$\lVert
A-CMR\rVert\leq(\lVert(V^{\top}\\!P)^{-1}\rVert+\lVert(S^{\top}U)^{-1}\rVert)\cdot\sigma_{k+1}.$

##  6 Experiments

We conduct numerical experiments to evaluate the empirical performance of the
DEIM scheme [20], the QDEIM procedure [11], the MaxVol method [15], and the
iterative subselection techniques discussed in this paper. The following are
the iterative subselection methods we evaluate:

  1. CADP-CX: refers to Algorithm 3.

  2. DADP-CX: represents Algorithm 4.

  3. DADP-CUR: corresponds to Algorithm 5.

  4. CADP-CUR: denotes the adapted version of Algorithm 3 with the residual defined as $E=A-CMR$.

To assess the effectiveness of our algorithms, we test them on various data
analysis tasks in different application domains, such as analyzing
recommendation systems, categorizing text and retrieving information, and
image compression. Our evaluation includes synthetic and real-world data
matrices, both sparse and dense, with varying sizes ranging from small to
large scale, which we summarize in Table 2. In the implementation, we perform
the column-pivoted QR factorization and the reduced SVD using the MATLAB
built-in functions qr and svd, respectively. For Algorithms 6 and 7, we use
our implementation of the Lanczos bidiagonalization method of [1] by
incorporating the threshold parameter $\delta$ and upper limit $\ell$ on the
number of singular vectors to be computed. Unless otherwise stated, in all the
experiments we use as default the number of rounds $t=10$, the parameter
$\delta=0.8$, and upper limit $\ell=k/10$.

Table 2: Various examples and dimensions considered. Exp. | Domain | Matrix | $m$ | $n$
---|---|---|---|---
1 | Synthetic | Sparse | 100000 |  12300
2 | Text categorization | Sparse |  123139 | 15210
3 | Text categorization | Sparse |  128293 | 18933
4 | Image compression | Dense |  113500 |  15000
5 | Recommendation system | Dense |  114116 |  12100
6 | Economic modeling | Sparse |  129610 | 29610
7 | Computational fluid dynamics problem | Sparse |  130412 | 30412
8 | Image compression | Dense |  150000 |  13072

######  Experiment 6.1.

In our first experiment, we investigate how different choices of $\delta$ and
the number of rounds $t$ affect the approximation accuracy of the various
iterative subselection strategies. We use the relative approximation error
${\lVert A-CMR\rVert}\,/\,{\lVert A\rVert}$ as the evaluation metric. For this
experiment just as in [20], we generate a sparse, nonnegative matrix
$A\in\mathbb{R}^{m\times n}$, with $m=100000$ and $n=300$, of the form

$A=\sum_{j=1}^{10}\frac{2}{j}\,\mathbf{x}_{j}\
\mathbf{y}_{j}^{\top}+\sum_{j=11}^{300}\frac{1}{j}\,\mathbf{x}_{j}\
\mathbf{y}_{j}^{\top},$

where $\mathbf{x}_{j}\in\mathbb{R}^{m}$ and $\mathbf{y}_{j}\in\mathbb{R}^{n}$
are sparse vectors with random nonnegative entries (i.e., $\mathbf{x}_{j}={\sf
sprand}(m,1,0.025)$ and $\mathbf{y}_{j}={\sf sprand}(n,1,0.025)$.

From Fig. 1 we observe that increasing the number of rounds $t$ or $\delta$
does not necessarily lead to a monotonic decrease in the approximation errors.
The result implies that one needs to carefully choose the parameter $\delta$
or the number of rounds to get the optimum advantage of using the iterative
subselection strategies. For this experiment, we also observe that using the
residual $E=A-CMR$ instead of $A-CC^{+}A$ for the iterative subselection
yields better approximation errors in the delta strategy while it produces
worse approximation errors in the constant number of columns strategy.

$2$$4$$6$$8$$10$$2.15$$2.2$$2.25$$2.3$$2.35$$2.4$$\cdot 10^{-2}$# of
rounds$\|A-CMR\|\ /\ \|A\|$CADP-CURCADP-
CX$0.4$$0.5$$0.6$$0.7$$0.8$$2.15$$2.2$$2.25$$2.3$$\cdot 10^{-2}$Deltas
($\delta$s)DADP-CURDADP-CX Figure 1: Relative approximation errors for the
various iterative subselection DEIM CUR approximation algorithms for $k=30$.
The right figure represents selecting a constant number of columns and rows
per iteration; the left is the delta strategy. In all cases, increasing the
number of rounds or delta does not lead to a monotonic decrease in the
approximation errors.

######  Experiment 6.2.

Our next experiment is to demonstrate that the iterative subselection
techniques yield better approximation results than one-round sampling. We
perform the experiment using four real data sets and report the relative
approximation error ${\lVert A-CMR\rVert}\,/\,{\lVert A\rVert}$ of each
algorithm on each data set.

The first two data sets are relevant to text categorization and information
retrieval applications. In such data analysis problems, a “bag of words”
approach is commonly employed to represent documents. We opt for the
Reuters-21578 text categorization collection, which comprises documents that
were featured on Reuters’ newswire in $1987$. This data set is extensively
used as a benchmark in the text classification community, consisting of
$21578$ documents categorized into $135$ categories. For our experiment, we
use the preprocessed data set, which has $18933$ unique terms and $8293$
documents [5]. We normalize the rows of the sparse matrix, which has
dimensions $8293\times 18933$, to have a unit length. The second data set, the
Internet term document data, is from the Technion Repository of Text
Categorization Datasets (TechTC) [13]. We use the test set $26$, which
consists of a collection of $139$ documents on two topics with $15210$ terms
describing each document333http://gabrilovich.com/resources/data/techtc/. As
in [20], the $139\times 15210$ TechTC matrix rows are scaled to have a unit
2-norm.

The third data set is the Gisette data. Gisette is a handwritten digit
recognition problem. The problem is to separate the highly confusable digits
‘4’ and ‘9’. The digits have been size-normalized and centered in a fixed-size
image of dimension $28\times 28$. The resulting data set is of dimension
$13500\times 5000$.

The next data set pertains to the recommendation system analysis field, where
the primary objective is to provide service or purchase suggestions to users.
Collaborative filtering is a commonly used technique in recommendation
systems, which involves recommending items to users that were previously liked
by customers with comparable preferences. The Jester data set is frequently
employed as a benchmark in recommendation system research [14]. The data set
comprises $73421$ users and their ratings for $100$ jokes. We limit our
analysis to users who have provided ratings for all $100$ jokes, resulting in
a $14116\times 100$ matrix. We normalize the matrix by subtracting the mean of
each column from all the entries in that column.

The final two data sets are standard test matrices specifically designed for
sparse matrix problems. These data matrices are sourced from the publicly
available SuiteSparse Matrix Collection. The data set referred to as g7jac100,
is derived from the “Overlapping Generations Model” used to study the social
security systems of the G7 nations. It is a sparse matrix with dimensions
$29610\times 29610$ and contains $335972$ numerically nonzero entries with a
low rank of $21971$. We utilize the invextr1-new matrix, associated with a
computational fluid dynamics problem, which has dimensions $30412\times 30412$
and represents a rank-$27502$ structure. It contains $1793881$ nonzero
entries.

From Fig. 2, we can see that in all cases our iterative subselection-based CUR
algorithm has a lower approximation error than all the _one-round_
deterministic index selection algorithms considered. We also observe that the
approximation error of the QDEIM and the MaxVol techniques do not always
decrease monotonically with increasing $k$ values. We choose the number of
rounds for the CADP-CUR and CADP-CX algorithms to be $t=10$, and the parameter
$\delta=0.8$ and upper limit $\ell=k/10$ for the DADP-CUR and DADP-CX
algorithms. The results of all four proposed iterative subselection algorithms
are comparable.

$10$$20$$30$$40$$50$$40$$50$$60$$70$$\cdot 10^{-2}$$\|A-CMR\|\ /\
\|A\|$TechTC$10$$20$$30$$40$$50$$20$$30$$40$$\cdot
10^{-2}$Reuters$10$$20$$30$$40$$50$$20$$30$$40$$\cdot 10^{-2}$$\|A-CMR\|\ /\
\|A\|$Gisette$10$$20$$30$$40$$50$$20$$40$$60$$\cdot
10^{-2}$Jester$100$$200$$300$$400$$500$$20$$40$$\cdot
10^{-2}$rank-$k$$\|A-CMR\|\ /\
\|A\|$g7jac100$500$$600$$700$$800$$900$$1{,}000$$6$$8$$10$$12$$14$$16$$\cdot
10^{-2}$rank-$k$invextr1-newDEIMMAXVOLQDEIMCADP-CURDADP-CURCADP-CXDADP-CX
Figure 2: Relative approximation errors as a function of $k$ for the various
iterative subselection DEIM CUR approximation algorithms compared with some
standard CUR approximation algorithms using real data sets.

######  Experiment 6.3.

As stated in Section 4, when dealing with large-scale matrices, performing a
full (even reduced) SVD in each iteration of algorithms 3, 4, and 5 can often
become excessively costly. To evaluate the efficiency of the proposed
algorithms, which compute the full SVD, compared to their respective large-
scale versions that utilize Lanczos bidiagonalization to find a limited number
of singular vectors (refer to Algorithm 7), we conduct experiments using a
large-scale data set. Specifically, we utilized the cifar-10 data set, which
consists of $60000$ color images sized $32\times 32$ pixels, divided into ten
different classes with $6000$ images per class. The data set is divided into
$50000$ training set images and $10000$ test set images. We focus on the
training data set containing $50000$ images. Each image in the data set has
been reshaped into a 1D array of length $3072$, resulting in a dense matrix of
size $50000\times 3072$. For our analysis, we approximate this matrix using a
rank-100 approximation.

Table 3 presents the results obtained from running the various algorithms. We
observe that the large-scale variants (i.e., the various adaptations of
Algorithm 7), which utilize an iterative method for computing a few SVDs,
demonstrate higher efficiency while maintaining comparable approximation
quality compared to the algorithms that compute the full SVD. Notably, both
for the full SVD and the iterative SVD routines, the algorithms with the
residual defined as $E=A-CMR$ exhibit greater efficiency than those with the
residual computed as $E=A-CC^{+}A$. Therefore, our new approach to computing
the residual for the iterative subselection proves to be more efficient than
the existing method while maintaining similar approximation accuracy for this
experiment.

Table 3: Comparison of large-scale iterative subselection algorithms (iterative method for computing few SVDs) and small-scale iterative subselection algorithms (full SVD computation) on the cifar-10 data set approximation. Method | Full SVD |  | Iterative SVD
---|---|---|---
| Relative error | Runtime (s) |  | Relative error | Runtime (s)
CADP-CUR | $0.028$ | $3.24\cdot 10^{2}$ |  | $0.028$ | $1.61\cdot 10^{2}$
DADP-CUR | $0.027$ | $5.06\cdot 10^{2}$ |  | $0.027$ | $1.66\cdot 10^{2}$
CADP-CX | $0.026$ | $6.21\cdot 10^{2}$ |  | $0.026$ | $5.40\cdot 10^{2}$
DADP-CX | $0.026$ | $1.15\cdot 10^{3}$ |  | $0.026$ | $5.78\cdot 10^{2}$

##  7 Conclusions

New approaches for selecting subsets of columns and rows using iterative
subselection strategies have been presented. The first one is a DEIM
adaptation of the so-called adaptive sampling [8] for column subset selection.
This procedure follows a fixed selection schedule, choosing a predetermined
number of columns or rows in each iteration. In contrast, the second proposed
iterative subselection strategy dynamically adjusts the selection schedule
based on the decay of the singular values of the data. This approach aims to
prioritize the selection of columns or rows that contribute the most to the
overall structure and information of the data. By considering the significance
of singular values and leveraging their decay pattern, the algorithm can adapt
to the unique characteristics of the data, resulting in a more data-driven
selection process.

Additionally, we also introduce an alternative approach for computing the
residual in the index selection process. The first two iterative subselection
algorithms we propose, i.e., Algorithm 3 and Algorithm 4, as well as existing
adaptive sampling procedures [4, 8, 9, 19, 24], define the residual as the
error resulting from projecting the matrix $A$ onto either the column space of
$C$ or the row space of $R$, i.e., $E=A-CC^{+}A$ or $E=A-AR^{+}R$,
respectively. In contrast, our new method defines the residual as the error
incurred by simultaneously projecting $A$ onto both the column space of $C$
and the row space of $R$. This entails computing a CUR factorization at each
step using only the selected columns and rows.

We have also discussed how iterative procedures for computing a few singular
vectors of large data matrices can be used with the newly proposed strategies.
We have presented an adaptation of Algorithm 4 for the large-scale case in
Algorithm 7, which can straightforwardly be adapted for Algorithm 3 and
Algorithm 5. For each of the iterative subselection strategies proposed in
this paper, we invoke the DEIM index selection method. However, we note that
other deterministic index selection schemes such as the QDEIM technique [11]
and the MaxVol procedure [15] may be employed. We have demonstrated through
empirical analysis that the proposed methods in this work can produce better
approximation results than the traditional method of _one-round sampling_ of
all columns and rows.

Overall, the proposed techniques may be useful for improving the accuracy of a
CUR decomposition, but may also introduce additional complexities that need to
be carefully addressed. The choice of whether to use the proposed iterative
subselection methods or not may depend on the specific problem or application,
as well as the trade-offs between accuracy, complexity, and computational
resources.

## Acknowledgements

This work has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
812912.

## References

  * Baglama and Reichel [2005] J. Baglama, L. Reichel,  Augmented implicitly restarted Lanczos bidiagonalization methods,  SIAM J. Sci. Comput. 27 (2005) 19–42.
  * Barrault et al. [2004] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera,  An ‘empirical interpolation’ method: Application to efficient reduced-basis discretization of partial differential equations,  Comptes Rendus Math. 339 (2004) 667–672.
  * Bischof and Hansen [1991] C. H. Bischof, P. C. Hansen,  Structure-preserving and rank-revealing QR-factorizations,  SIAM J. Sci. Comput. 12 (1991) 1332–1350.
  * Boutsidis and Woodruff [2017] C. Boutsidis, D. P. Woodruff,  Optimal CUR matrix decompositions,  SIAM J. Comput. 46 (2017) 543–589.
  * Cai et al. [2005] D. Cai, X. He, J. Han,  Document clustering using locality preserving indexing,  IEEE Trans. Knowl. Data Eng. 17 (2005) 1624–1637.
  * Chandrasekaran and Ipsen [1994] S. Chandrasekaran, I. C. F. Ipsen,  On rank-revealing factorisations,  SIAM J. Matrix Anal. Appl. 15 (1994) 592–622.
  * Chaturantabut and Sorensen [2010] S. Chaturantabut, D. C. Sorensen,  Nonlinear model reduction via discrete empirical interpolation,  SIAM J. Sci. Comput. 32 (2010) 2737–2764.
  * Deshpande et al. [2006] A. Deshpande, L. Rademacher, S. S. Vempala, G. Wang,  Matrix approximation and projective clustering via volume sampling,  Theory Comput. 2 (2006) 225–247.
  * Deshpande and Vempala [2006] A. Deshpande, S. Vempala,  Adaptive sampling and fast low-rank matrix approximation,  in: Proceedings of the 10th RANDOM APPROX, 2006, pp. 292–303.
  * Drineas et al. [2008] P. Drineas, M. W. Mahoney, S. Muthukrishnan,  Relative-error CUR matrix decompositions,  SIAM J. Matrix Anal. Appl. 30 (2008) 844–881.
  * Drmac and Gugercin [2016] Z. Drmac, S. Gugercin,  A new selection operator for the discrete empirical interpolation method—Improved a priori error bound and extensions,  SIAM J. Sci. Comput. 38 (2016) A631–A648.
  * Frieze et al. [2004] A. Frieze, R. Kannan, S. Vempala,  Fast monte-carlo algorithms for finding low-rank approximations,  J. ACM 51 (2004) 1025–1041.
  * Gabrilovich and Markovitch [2004] E. Gabrilovich, S. Markovitch,  Text categorization with many redundant features: Using aggressive feature selection to make svms competitive with c4.5,  in: The 21st International Conference on Machine Learning, 2004, p. 41.
  * Goldberg et al. [2001] K. Goldberg, T. Roeder, D. Gupta, C. Perkins,  Eigentaste: A constant time collaborative filtering algorithm,  Inform. Retrieval 4 (2001) 133–151.
  * Goreinov et al. [2010] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, N. L. Zamarashkin,  How to find a good submatrix,  in: Matrix Methods: Theory, Algorithms And Applications, World Scientific, Singapore, 2010, pp. 247–256.
  * Gu and Eisenstat [1996] M. Gu, S. C. Eisenstat,  Efficient algorithms for computing a strong rank-revealing QR factorization,  SIAM J. Sci. Comput. 17 (1996) 848–869.
  * Guruswami and Sinop [2012] V. Guruswami, A. K. Sinop,  Optimal column-based low-rank matrix reconstruction,  in: Proc. Annu. ACM-SIAM Symp., 2012, pp. 1207–1214.
  * Mahoney and Drineas [2009] M. W. Mahoney, P. Drineas,  CUR matrix decompositions for improved data analysis,  Proc. Natl. Acad. Sci. USA 106 (2009) 697–702.
  * Paul et al. [2015] S. Paul, M. Magdon-Ismail, P. Drineas,  Column selection via adaptive sampling,  Adv. Neural Inf. Process Syst. 28 (2015).
  * Sorensen and Embree [2016] D. C. Sorensen, M. Embree,  A DEIM induced CUR factorization,  SIAM J. Sci. Comput. 33 (2016) A1454–A1482.
  * Stewart and Sun [1990] G. W. Stewart, J. G. Sun, Matrix perturbation theory, Academic press, 1990.
  * Strang [2019] Strang, G., Linear Algebra and Learning from Data, SIAM, Philadelphia, 2019.
  * Voronin and Martinsson [2017] S. Voronin, P. G. Martinsson,  Efficient algorithms for CUR and interpolative matrix decompositions,  Adv. Comput. Math. 43 (2017) 495–516.
  * Wang and Zhang [2013] S. Wang, Z. Zhang,  Improving CUR matrix decomposition and the Nyström approximation via adaptive sampling,  J. Mach. Learn. Res. 14 (2013) 2729–2769.