1

Block CUR : Decomposing Large Distributed
Matrices

arXiv:1703.06065v1 [stat.ML] 17 Mar 2017

Urvashi Oswal, Swayambhoo Jain, Kevin S. Xu, and Brian Eriksson

Abstract
A common problem in large-scale data analysis is to approximate a matrix using a combination of specifically
sampled rows and columns, known as CUR decomposition. Unfortunately, in many real-world environments, the
ability to sample specific individual rows or columns of the matrix is limited by either system constraints or cost. In
this paper, we consider matrix approximation by sampling predefined blocks of columns (or rows) from the matrix.
This regime is commonly found when data is distributed across multiple nodes in a compute cluster, where such
blocks correspond to columns (or rows) of the matrix stored on the same node, which can be retrieved with much less
overhead than retrieving individual columns stored across different nodes. We propose a novel algorithm for sampling
useful column blocks and provide guarantees for the quality of the approximation. We demonstrate the practical utility
of this algorithm for computing the block CUR decomposition of large matrices in a distributed setting using Apache
Spark. Using our proposed block CUR algorithms, we can achieve a significant speed-up compared to a regular CUR
decomposition with the same quality of approximation.

I. I NTRODUCTION
The ability to perform large-scale data analysis is often limited by two opposing forces. The first force is the
need to store data in a matrix format for the purpose of analysis techniques such as regression or classification.
The second force is the inability to store the data matrix completely in memory due to the size of the matrix
in many application settings. This conflict gives rise to storing factorized matrix forms, such as SVD or CUR
decompositions [1].
We consider a matrix A with m rows and n columns, i.e., A ∈ Rm×n . Using a truncated k number of
singular vectors (e.g., where k < min {m, n}), the singular value decomposition (SVD) provides the best rank-k
approximation to the original matrix. The singular vectors of a SVD often do not preserve the structure in original
data. Preserving the original structure in the data may be desirable due to many reasons including interpret-ability
in case of biometric data or for storage efficiency in case of sparse matrices. This has led to the introduction of
the CUR decomposition, where the factorization is performed with respect to a subset of the rows and columns of
the matrix itself. This specific decomposition describes the matrix A as the product of a subset of the matrix rows
R and a subset of the matrix columns C (along with a matrix U that fits to A).
UO is with the Department of Electrical and Computer Engineering, University of Wisconsin – Madison, SJ is with the Department of Electrical
and Computer Engineering, University of Minnesota – Twin Cities, KSX is with the Electrical Engineering and Computer Science Department,
University of Toledo, BE is with Technicolor Research – Los Altos. Author emails: uoswal@wisc.edu, jainx174@umn.edu,
kevin.xu@utoledo.edu, brian.eriksson@technicolor.com

2

Fig. 1: (Left) Distributed storage of a large matrix across multiple nodes in a cluster. Blocks of columns are allocated to each of
the G nodes. (Right) Example block CUR decomposition, where C t ∈ Rm×s for t ∈ [g] is sampled from {A(jt ) : jt ∈ [G]}.

Significant prior work has examined how to efficiently choose the rows and columns in the CUR decomposition
and has derived worst-case error bounds (e.g., [2]). These methods have been applied successfully to many real-world
problems including genetics [3], astronomy [4], and mass spectrometry imaging [5]. Unfortunately, these large-scale
datasets require distributed storage, a regime where there can be substantial overhead involved in querying individual
rows or columns of a matrix. In these regimes, it is more efficient to retrieve predefined blocks of rows or columns
at one time corresponding to the rows or columns stored on the same node, as shown in Figure 1(a), in order to
minimize the overhead. In doing so, one forms a block CUR decomposition, as shown in Figure 1(b). Current CUR
decomposition techniques do not take advantage of this predefined block structure. Another setting where block
CUR decompositions are useful involves processing data where the ordering of rows or columns is meaningful,
such as images, video, or speech data matrices. In this setting, sampling contiguous blocks of columns is necessary
for interpretability of the factorized representation.
Using these insights into real-world applications of CUR decomposition, this paper makes a series of contributions.
We propose a randomized Block CUR algorithm for subset selection of rows and blocks of columns and derive
worst-case error bounds for this randomized algorithm. This algorithm performs fast block sampling taking advantage
of the natural storage of matrices in distributed environments. Additionally, this algorithm does not require a full
SVD calculation, saving significant computation time in the large-scale setting. Finally, we present worst-case error
bounds for this algorithm using theoretical results related to approximating matrix multiplication and generalized
`2 -regression in this sampling blocks setting.
Additionally, we demonstrate empirically that the proposed block CUR algorithms can achieve a significant speedup when used to decompose large matrices in a distributed data setting. We conduct a series of CUR decomposition
experiments in Amazon AWS using both synthetic and real-world data. In this distributed environment, we find
that our Block CUR approach achieves a speed-up of 2x to 6x for matrices larger than 12000 × 12000. This is
compared with previous CUR approaches that sample individual rows and columns and while achieving the same
matrix approximation error rate.
The rest of the paper is organized as follows. The CUR decomposition model is introduced in Section II. We detail
the specific steps of our block CUR algorithm and corresponding worst-case relative error bounds in Section III.

3

We present experiments using our proposed block CUR algorithms in a distributed setting using Apache Spark in
Section IV. Finally, Section V introduces theoretical results used in the proofs of the relative error bounds.
II. S ETUP AND BACKGROUND
The need to factorize a matrix using a collection of rows and columns of that matrix has motivated the CUR
decomposition literature. CUR decomposition is focused on sampling rows and columns of the matrix to provide
a factorization that is close to the best rank-k approximation of the matrix. One of the most fundamental results
for a CUR decomposition of a given matrix A ∈ Rm×n was obtained in [1]. This relative error bound result is
summarized in the following theorem.
Theorem 1. (Theorem 2 from [1] applied to AT ) Given A ∈ Rm×n and an integer k ≤ min{m, n}, let r =
2

2

O( kε2 ln( 1δ )) and c = O( rε2 ln( 1δ )). There exist randomized algorithms such that, if c columns are chosen to construct
C and r rows are chosen to construct R, then with probability ≥ 1 − δ, the following holds:
kA − CU RkF ≤ (1 + ε)kA − Ak kF
where ε, δ ∈ (0, 1), U = W † and W is the scaled intersection of C and R.
This theorem states that as long as enough rows and columns of the matrix are acquired (r and c, respectively), then
the CUR decomposition will be within a constant factor of the error associated with the best rank-k approximation
of that matrix.
Central to the proposed randomized algorithm was the concept of sampling columns of the matrix based on a
leverage score. The leverage score measures the contribution of each column to the approximation of A.
Definition 1. The leverage score of a column is defined as the squared row norm of the top-k right singular vectors
of A corresponding to the column:
`j = kV TA,k ej k22 , j ∈ [n],
where V A,k consists of the top-k right singular vectors of A as its rows, and ej picks the j-th column of V TA,k .
The CUR algorithm involves randomly sampling r rows using probabilities generated by the calculated leverage
scores to obtain the matrix R, and thereafter sampling c columns of A based on leverage scores of the R matrix
to obtain C.
The key technical insight in [1] is that the leverage score of a column measures “how much” of the column lies
in the subspace spanned by the top-k left singular vectors of A. By sampling columns that lie in this subspace
more often, we get a relative-error low rank approximation of the matrix.
Unfortunately, these prior results require both the calculation of SVD of A and the sampling of arbitrary rows
and columns of the matrix A. Both requirements may be either unrealistic or inefficient. In this paper, we focus on
the problem of efficiently sampling pre-defined blocks of columns (or rows) of the matrix to provide a factorization

4

that is close to the best rank-k approximation of the matrix. In the following section, we propose and analyze a
randomized algorithm for sampling blocks of the matrix based on block leverage scores.
III. T HE B LOCK CUR A LGORITHM
A block is defined as a collection of s columns1 . Let there be G = dn/se possible blocks in A. We consider the
blocks to be predefined due to natural constraints or cost, such as data partitioning in a distributed compute cluster.
The goal of the block CUR algorithm is to approximate the underlying matrix A using g blocks of columns and
r rows, as represented in Figure 1.
Given the new regime of submatrix blocks, we must first define new block-specific variables for leverage scores,
matrix incoherence, and matrix stable rank. Informed by the prior CUR literature, we begin by defining a block
leverage score for each block of columns.
Definition 2. The block leverage score of a group of columns is defined as the sum of the squared row norms of
the top-k right singular vectors of A corresponding to the columns in the block:
`g (A, k) = kV TA,k E g k2F , g ∈ [G],
where V A,k consists of the top-k right singular vectors of A, and E g picks the columns of V TA,k corresponding
to the elements in block g.
Much like the individual column leverage scores defined in [1], the block leverage scores measure how much a
particular column block contributes to the approximation of the matrix A.
We next define a property of matrix rank relative to the collection of matrix blocks. Specifically, we focus on
the concept of matrix stable rank from [6]. With respect to our matrix blocks, we define the block stable rank as
the minimum stable rank across all matrix blocks.
Definition 3. Let V A,k consist of the top-k right singular vectors of A then the block stable rank is defined as
αA = min kV TA,k E g k2F /kV TA,k E g k22 .
g∈[G]

where E g picks the columns of V TA,k corresponding to the elements in block g.
Intuitively, this gives a measure of how informative the worst matrix column block is.
Finally, when we sample rows uniformly at random, we can give relative error approximation guarantees when
the matrix A satisfies a notion of column space incoherence. This avoids pathological constructions of rows of A
that cannot be sampled at random.
Definition 4. The top-k column space incoherence is defined as
µ := µ(U TA,k ) =

m
max kU TA,k ei k22
k i

1 To keep notation simple, we will consider just sampling column blocks, although the techniques and derivations also hold for row blocks
by applying the same to the transpose.

5

where ei picks the i-th column of U TA,k .
The column space incoherence is used to provide a guarantee for approximation without computing the SVD of
the entire matrix A.
A. Algorithm Details
The Block CUR Algorithm, detailed in Algorithm 1, takes as input the matrix A and returns as output an r × n
matrix R consisting of a small number of rows of A and an m × c matrix C consisting of a small number of
column blocks from A.
Algorithm 1: Block CUR
Input : A, target rank k, size of each block s, error parameter ε, positive integers r, g
b = CU R
Output: C, R, A
1) Row subset selection: Sample r rows uniformly from A according to pi = 1/m for i ∈ [m] and compute
R = S TR A.
2) Column block subset selection: For t = {1, 2, . . . , g}, select a block of columns jt ∈ [G] independently
kV T

E i k2

R,r
F
=
for i ∈ [G] and update S, where V R,r consists of
with probability pjt = P r[jt = i] = `i (R,r)
r
r
the top-r right singular vectors of R, and E i picks the columns V TR,r corresponding to the elements in
block i. Compute C = AS.
b = CU R where U = W † , and W = RS is the scaled intersection of R and C.
3) CUR approximation: A

In Algorithm 1, C = AS, where S ∈ Rn×gs is the block scaling and sampling matrix and the (jt , t)-th non-zero
s × s block of S is defined as
S jt ,t = √

1
Is
gpjt

where g = c/s is the number of blocks picked by the algorithm. An example of the sampling matrix S with blocks
chosen in the order [1, 3, 2] is as follows:


√1 I
 gp1 s

S n×gs




=





0

0

0

0

√1 Is
gp3

0

0

0





√ 1 I s

gp2
.

0 


0

This sampling matrix picks the blocks of columns and scales each block to compute C = AS. A similar sampling
and scaling matrix S R is defined to pick the blocks of rows and scale each block to compute R = S TR A.
In addition to considering block sampling of columns, another key advantage of this algorithm is not requiring
the computation of a full SVD of A. In many applications, it may not feasible to compute the SVD of the entire
matrix A. In these cases, algorithms requiring knowledge of the leverage scores cannot be used. Instead, we use an
estimate of the block leverage scores called the approximate block leverage scores. A subset of the rows are chosen
uniformly at random, and the block scores are calculated using the top-k right singular vectors of this row matrix
instead of the entire A matrix. This step can also be replaced with other fast approximate calculations of leverage

6

r

Method
Traditional CUR

O



Block CUR

O



k2
ε2
k2
ε2

log

1
δ




log( 1δ )

 4 c

O kε6 log3 ( 1δ )

 4
3 1
log
(
)
O αsk
6
δ
Rε

No.
sampling
ops. 
 of

3 1
k2
k4
O
+
log
(δ)
2
ε6


 2ε
4
3 1
k
k
+
log
(
)
O
2
6
ε
αR ε
δ

TABLE I: Table showing the number of rows r, columns c, and total number of sampling operation needed for given  in
traditional CUR and block CUR. For ease of comparison, we show the block CUR results from Corollary 1 with full SVD
computation.

scores involving sketching or additional sampling [7], [8]. The advantage of using our approximate leverage scores
is that the same set of rows is used to approximate the scores and also to compute the CUR approximation. Hence
no additional sampling or sketching steps are required.
B. Theoretical Results and Discussion
The main technical contribution of the paper is a relative-error bound on the quality of approximation using
blocks of columns or rows to approximate a matrix A ∈ Rm×n . The following result provides a relative-error
guarantee for the Block CUR approximation using Algorithm 1,
2

Theorem 2. Given A ∈ Rm×n with incoherent top-k column space, i.e. µ ≤ µ0 , let r = O(µ0 kε2 ln( 1δ )) and
2

g = O( αRr ε2 ln( 1δ )). There exist randomized algorithms such that, if r rows and g column blocks are chosen to
construct R and C, respectively, then with probability ≥ 1 − δ, the following holds:
kA − CU RkF ≤ (1 + ε)kA − Ak kF
where ε, δ ∈ (0, 1) and U = W † is the pseudoinverse of the scaled intersection of C and R.
We provide a sketch of the proof in Section V and provide the details in the appendix.
We now detail the differences between our technique and prior CUR algorithms. This includes additional
assumptions required, algorithmic trade-offs, and discussion of sampling and computational complexity.
1) Block Stable Rank: The theorem tells us that the number of blocks required to achieve an ε relative error
depends on the structure of the blocks (through αR ). Intuitively, this is saying the groups that provide more
information improve the approximation faster than less informative groups. The αR term depends on the stable
or numerical rank (a stable relaxation of exact rank) of the blocks. The stable rank r = r(A) = kAk2F /kAk22
is a relaxation of the rank of the matrix, in fact it is stable under small perturbations of the matrix A [6]. For
instance, the stable rank of an approximately low rank matrix tends to be low. The αR term defined in Theorem
2 is the minimum stable rank of the column blocks. Thus, the αR term gives a dependence of the block sampling
complexity on the stable ranks of the blocks.
It is easy to check that 1 ≤ αR ≤ s. In the best case, when all the groups have full stable rank with equal singular
values, αR achieves it’s maximum. The worst case αR = 1 is achieved when a group or block is rank-1. That is,
sampling groups of rank s gives us a lot more information than groups of rank 1, which leads to a reduction in the
final group sampling complexity.

7

2) Incoherence: The column space incoherence (Definition 4) is used to provide a guarantee for approximation
without computing the SVD of the entire matrix A. However, if it is possible to compute the SVD of the entire
matrix, then the rows can be sampled using row leverage scores, and the incoherence assumption can be dropped.
The relative error guarantee for the full SVD Block CUR approximation is stated in the following corollary.
2

2

Corollary 1. Given A ∈ Rm×n , let r = O( kε2 ln( 1δ )) and g = O( αRr ε2 ln( 1δ )). There exist randomized algorithms
such that, if r rows and g column blocks are chosen to construct R and C, respectively, then with probability
≥ 1 − δ, the following holds:
kA − CU RkF ≤ (1 + ε)kA − Ak kF
where ε, δ ∈ (0, 1), and U = W † is the pseudoinverse of the scaled intersection of C and R.
The corollary follows by similar analysis as Theorem 2 so we defer the proof to the appendix. Other than block
sampling, the setup of this result is equivalent to the traditional column sampling result stated in Theorem 1. Next,
we will use this result to compare block sampling with traditional column sampling.
3) Complexity: The running time of Algorithm 1 is essentially driven by the time required to compute the SVD
of R and the construction of R, C and U . Algorithm 1 requires O(SV D(R)) time, construction of R requires
O(rn) time, construction of C takes O(mc) time, construction of W requires O(rc) time and construction of U
takes O(r2 c) time.
The sampling complexity of the Block CUR derived in Theorem 2 and the full SVD Block CUR in Corollary
1 tells us the number of sampling operations or queries that need to be made to the memory in order to construct
the R and C matrices. In order to compare complexity of block sampling with individual column sampling we
focus our attention on the full SVD Block CUR result in Corollary 1. As shown in Table I the total number of
individual columns queried by traditional CUR is always less than or equal to those required by block CUR because
1 ≤ αR ≤ s. However, recall that in the distributed setting the key advantage of sampling column blocks is the
efficieny of retrieving blocks of rows or columns at one time corresponding to the rows or columns stored on the
same node. Thus, in Table I we see that the number of queries made by the Block CUR algorithm are much less
than traditional CUR. In the best case scenario, we can get a speedup proportional to the block size (when αR = s).
In the corner case when each block contains only one important column, block sampling may sample s times
more columns in total than column sampling to reach the same error (since αR = 1 and individual sampling has
the flexibility to pick only the important columns). However, in this case (with αR = 1) the number of queries or
sampling operations for Block CUR will be equal to those made by traditional CUR thus achieving the same time
complexity.
In the next section, we demonstrate the speedup achieved by computing the block CUR decomposition of large
matrices in a distributed setting using Apache Spark.
IV. N UMERICAL RESULTS
In this section we provide empirical comparisons between the block CUR and traditional CUR algorithms on
both synthetic and real-world data. We report the relative-error of the decomposition (i.e., kA − CU RkF /kAkF )

8

1

1

Block CUR
Individual CUR

0.8

0.6

0.4

0.6

0.4

0.2

0.2

0

5000

10000

15000

20000

25000

30000

0.4

0

0

35000

0.6

0.2

0

0

Block CUR
Individual CUR

0.8

Relative error

Relative error

Relative error

0.8

1

Block CUR
Individual CUR

10000

20000

30000

40000

50000

0

50000 100000 150000 200000 250000 300000 350000

Computation Time ( in ms )

Computation Time ( in ms )

Computation Time ( in ms )

Fig. 2: (Left) Performance on 8000 × 8000 matrix with rank 800. (Center) Performance on 12000 × 12000 matrix with rank
1200, and (Right) Performance on 20000 × 20000 matrix with rank 2000.

6
5

Speedup

4
3
2
1
0
0

50

100

150

200

250

Block size

Fig. 3: Runtime speed-up from block sampling when compared to individual column sampling for varying block sizes.

and the running time of each algorithm on different data-sets.
We implemented the algorithms in Scala 2.10 and Apache Spark 2.11 on Amazon Elastic Map-Reduce (Amazon
EMR). The compute cluster was constructed using four Amazon m4.4xlarge instances, with each compute node
having a Intel Xeon E5-2676 v3 processor and 64 GB of RAM. Using Spark, we store the data sets as resilient
distributed dataset (RDD), a collection of elements partitioned across the nodes of the cluster (see Figure 1(a)). In
other words, Spark partitions the data into many blocks and distributes these blocks across multiple nodes in the
cluster. Using block sampling, we can approximate the matrix by sampling only a subset of the important blocks.
Meanwhile, individual column sampling would require looking up all the partitions containing specific columns of
interest. Our experiments examine the runtime speed-up from our block sampling CUR that exploits the partitioning
of data.
A. Synthetic Experiments
The synthetic data is generated by A = U V where U ∈ Rm×k and V ∈ Rk×n are random matrices with i.i.d.
Gaussian random entries, resulting in a low rank matrix A. We perform CUR decomposition on matrices of size
m × n with m = n, target rank k, and number of blocks G (set here across all experiments to be 100). The leverage
scores are calculated by computing the SVD of the rows sampled uniformly with R ∈ Rr×n . We sample one-sixth
of the rows in all experiments.
Figure 2 shows the plots for the relative error achieved with respect to the runtime required to sample the C
and R matrices for both the block CUR and traditional CUR algorithms. To focus on the speed-up achieved by

9

0.22
Block CUR
Individual CUR

0.2

Relative error

0.18
0.16
0.14
0.12
0.1
0.08
0

5000

10000

15000

20000

25000

30000

35000

Computation Time ( in ms )

Fig. 4: Performance on 900 × 10000 Arcene dataset with block size 12.

taking into account the block storage of the data we compare the running time of only the sampling operations of
the algorithms (which excludes the time required to compute the SVD). We note that other steps in both algorithms
can be updated to include faster variants such as the approximation of the leverage scores by sketching or sampling
[7].
We vary g, the number of blocks chosen, from 1 to 6. The number of columns chosen is thus c = gs, where
s denotes the number of columns in a block and varies from 50 to 200. We repeat each algorithm (block CUR
and traditional CUR) twice for the specified number of columns, with each realization as a point in the plot. The
proposed block CUR algorithm samples the c columns in g blocks, while the traditional CUR algorithm samples
the c columns one at a time.
Consistently, these results show that block sampling achieves the relative error much faster than the individual
column sampling – with performance gains increasing as the size of the matrix grows, as shown in Figure 3. While
the same amount of data is being transmitted regardless of whether block or individual column sampling is used,
block sampling is much faster because it needs to contact fewer executors to retrieve blocks of columns rather the
same number of columns individually. In the worst case, sampling individual columns may need to communicate
with all of the executors, while block sampling only needs to communicate with g executors. Thus, by exploiting
the partitioning of the data, the block CUR approach is able to achieve roughly the same quality of approximation
as traditional column-based CUR, as measured by relative error, with significantly less computation time.
B. Real-World Experiment
We also conduct experiments on the Arcene dataset [9] which has 900 rows and 10000 columns. We compare
the running time for both block and traditional CUR decomposition. We again find consistent improvements for
the block-wise approach compared with individual column sampling. With block size s = 12, sampling up to 10
groups led to an average speed up of 11.22 over individual column sampling, as shown in Figure 4. The matrix is
very low rank and sampling a few groups gave small relative errors.
V. P ROOF OF MAIN RESULT
In this section we provide the proof sketch for the relative-error guarantee of approximation provided by
Theorem 2. In proving this results, our theoretical results reveal the following contributions:

10

•

As a consequence of the proof of Lemma 1, we show that using block leverage scores if enough blocks are
sampled with high probability
kA − AS(RS)† RkF ≤ (1 + ε)kA − AR† RkF .
This gives a guarantee on the approximate solution obtained by solving a block-sampled regression problem
min

X∈Rm×r

k(AS) − X(RS)kF

instead of the entire least squares problem.
•

In Lemma 2, we show that the multiplication of two matrices A and B can be approximated by the product
of smaller sampled and scaled blocks of the matrices.

•

As a special case of the above result, when R = A we get a bound for the block column subset selection

2
problem. If g = O( αkA ε2 log 1δ ) blocks are chosen, then with probability at least 1 − δ we have
kA − CC † AkF ≤ (1 + ε)kA − Ak kF .

A. Proof Sketch
We start with Lemma 1, which states a non-boosting approximation error result for Algorithm 1.
2

2

Lemma 1. Given A ∈ Rm×n with incoherent top-k column space, i.e. µ ≤ µ0 , let r = O(µ0 kε2 ) and g = O( αRr ε2 ).
If rows and column blocks are chosen according to Algorithm 1, then with probability at least 0.7, the following
holds:
kA − CU RkF ≤ (1 + ε)kA − Ak kF
where ε ∈ (0, 1) and U = W † is the pseudoinverse of the scaled intersection of C and R.
Proof Sketch: First, note U = (RS)† and C = AS.
kA − CU RkF = kA − AS(RS)† RkF
Recall that R ∈ Rr×n has rank no greater than r; A ∈ Rm×n ; ε ∈ (0, 1); and that the same column blocks from
R and A are picked with the following probability distribution:
pi =

kV TR,r E i k2F
`i (R, r)
=
, ∀i ∈ [G].
r
r

2

Here, we show that if g = O( αRr ε2 ) blocks are chosen, then with probability at least 0.85 we have
kA − AS(RS)† RkF ≤ (1 + ε)kA − AR† RkF .
First step of the proof is to show that V TR,r S is full rank. Suppose δ 0 = 0.05.
Step 1: Using Lemma 2 (stated after this sketch) with β = 1, if g = O(r2 ln(1/δ 0 )/αR ε21 ) we get the following
with probability ≥ 1 − δ 0 we have kV TR,r V R,r − V TR,r SS T V R,r k2 ≤

ε1
2 .

Further, noting that V TR,r V R,r = I R,r

11

gives us a bound on the singular values of V TR,r S, for all i,
|1 − σi2 (V TR,r S)| = |σi (V TR,r V R,r ) − σi (V TR,r SS T V R,r )|
≤ kV TR,r V R,r − V TR,r SS T V R,r k2 ≤

ε1
2

Thus, it follows for all singular values of V TR,r S
r
r
ε1
ε1
1−
≤ σi (V TR,r S) ≤ 1 +
2
2

(1)

If ε1 < 1, then we have, with probability ≥ 1 − δ 0 , that V TR,r S is full rank.
Step 2: Using this result, we show that kA − AS(RS)† RkF can be bounded by the sum of three terms to
separate the random part and matrix perturbation part,
kA − AS(RS)† RkF = kA − AS(V TR,r S)† V TR,r kF
⊥T
⊥
⊥T
≤ kAV ⊥
R,r V R,r kF + kAV R,r V R,r SkF kΩk2
⊥T
T
+ kAV ⊥
R,r V R,r SS V R,r kF

where Ω = (V TR,r S)† − (V TR,r S)T . From Lemma 2 and using the fact that V TR,r S is full rank along with basic
linear algebra, we bound these terms separately with probability ≥ 1 − δ 0 as follows:
kQSkF ≤ (1/δ 0 )kQkF for many Q
√
kΩk2 ≤ ε1 / 2
√
r
T
⊥T
⊥T
SS
V
k
≤
V
kAV ⊥
kAV ⊥
√
R,r
F
R,r R,r
R,r V R,r kF
0
δ αR g

(2)
(3)
(4)

The proofs are provided in the appendix. Thus, we can conclude the following with probability ≥ 1 − 3δ 0 = 0.85:
kA − AS(RS)† RkF


 
1
1
⊥T
√ +
√
≤ 1+
ε1 kAV ⊥
R,r V R,r kF
δ0 2 δ0 2 2
≤ (1 + ε0 ) kA − AR† RkF
For δ 0 = 0.15/3 = 0.05, we have ε0 = 2ε1 /0.05.
So far we have shown the following
kA − CU RkF = kA − AS(RS)† RkF
≤ (1 + ε0 )kA − AR† RkF
Next, we bound kA − AR† RkF . Since A has incoherent column space, the uniform sampling distribution
pj = 1/m satisfies eqn. (13) in [1] with β = 1/µ0 . Consequently, we can apply modified version of Theorem 1 in
[1] to get with probability at least 0.85, kA − AR† RkF ≤ (1 + ε)kA − Ak kF . Finally, we get with probability

12

at least 0.7
kA − CU RkF ≤ (1 + ε0 )2 kA − Ak kF
≤ (1 + ε00 )kA − Ak kF

letting ε00 = 3ε0

This completes the proof of Lemma 1.
The result in Theorem 2 follows by applying standard boosting methods to Lemma 1 and running Algorithm 1
t = ln( 1δ ) times. By choosing the solution with minimum error and observing that 0.3 < 1/e, we have that the
relative error bound holds with probability greater than 1 − e−t = 1 − δ. The details of all the proofs are provided
in the appendix.
B. Block Multiplication Lemma
Finally, we state an intermediate result that was used in the proof of Lemma 1. It shows that the multiplication
of two matrices A and B can be approximated by the product of the smaller sampled and scaled block matrices.
This is the key lemma in proving the main result.
Lemma 2. Let A ∈ Rm×n , B ∈ Rn×p , and ε, δ ∈ (0, 1). Suppose αA = mink



kA(k) kF
kA(k) k2

2

. Construct C and R

using the sampling probabilities pi given by
kA(i) k2F
pi ≥ β PG
(j) 2
j=1 kA kF
for all i ∈ [G] and a constant β ∈ (0, 1]. Then, with probability at least 1 − δ,
1
kAkF kBkF .
kAB − CRkF ≤ √
δ βαA g
The proof is provided in the appendix.
VI. R ELATED WORK
The concept of sampling the important columns of a matrix based on the notion of subspace sampling first
appeared in context of fast `2 regression in [10] and was refined in [1] to obtain performance error guarantees for
CUR matrix decomposition. These guarantees were subsequently improved in follow-up work [2].
Modified versions of this problem have been studied extensively for adaptive sampling [11], divide-and-conquer
algorithms for parallel computations [12], and input-sparsity algorithms [13]. The authors of [11] propose an
adaptive sampling-based algorithm which requires only c = O(k/ε) columns to be sampled. The authors of [13]
also proposed an optimal, deterministic CUR algorithm. In [14], the authors prove the lower bound of the column
selection problem; at least c = k/ε columns are selected to achieve the (1 + ε) ratio. This rate is achieved by
more complicated column sampling algorithms, which can be extended to block sampling to achieve tighter bounds
on the sampling complexity. We stick to subspace sampling algorithms for the sake of simplicity and defer the
extension of more complicated adaptive block sampling algorithms to future work.

13

Recently, there has been emphasis on settings with the matrix partially known. For example, a CUR algorithm
and active column subset selection for partially observed matrices were proposed in [8] and [15], respectively. Our
work differs from [8] due to their requirement of prior uniform sampling of the matrix to estimate the leverage
scores, while our algorithm does not require these additional samples. With respect to [15], this prior work is
focused on column subset selection and does not provide bounds for the CUR approximation.
In contrast to all known prior work, we focus on the block setting where a block of columns is sampled rather
than sampling a single column, motivated by the distributed data setting. We analyze this problem by extending the
notion of subspace sampling to the block setting and give relative error guarantees for block CUR decomposition.
To best of our knowledge, this paper is the first to analyze the problem of column block subset selection.
Recent work in [16] considers a distributed implementation of CUR matrix decomposition on Apache Spark
and demonstrate its capability on very large data sets. In contrast, in this paper we consider the more natural
environment of block sampling for distributed computation, explore the performance advantages of block sampling
over individual column sampling, and provide the first theoretical error guarantees for block CUR decomposition.
VII. C ONCLUSIONS AND F UTURE W ORK
In this paper we extended the problem of CUR matrix decomposition to the block setting which is naturally
relevant to distributed storage systems. We proposed a novel algorithm and derived its performance bounds. We
demonstrated its pratical utility on real-life distributed storage systems. Some possible future directions for this
work include calculating the leverage scores quickly or adaptively, and considering the algorithms and error bounds
when the matrix has a pre-specified structure like sparsity.
R EFERENCES
[1] Petros Drineas, Michael W Mahoney, and S Muthukrishnan. Relative-error cur matrix decompositions. SIAM Journal on Matrix Analysis
and Applications, 30(2):844–881, 2008.
[2] Michael W Mahoney and Petros Drineas. Cur matrix decompositions for improved data analysis. Proceedings of the National Academy
of Sciences, 106(3):697–702, 2009.
[3] Peristera Paschou, Elad Ziv, Esteban G Burchard, Shweta Choudhry, William Rodriguez-Cintron, Michael W Mahoney, and Petros Drineas.
Pca-correlated snps for structure identification in worldwide human populations. PLoS Genet, 3(9):e160, 2007.
[4] Ching-Wa Yip, Michael W Mahoney, Alexander S Szalay, István Csabai, Tamás Budavári, Rosemary FG Wyse, and Laszlo Dobos. Objective
identification of informative wavelength regions in galaxy spectra. The Astronomical Journal, 147(5):110, 2014.
[5] Jiyan Yang, Oliver Rbel, Michael W Mahoney, and Benjamin P Bowen. Identifying important ions and positions in mass spectrometry
imaging data using cur matrix decompositions. Analytical chemistry, 87(9):4658–4666, 2015.
[6] Mark Rudelson and Roman Vershynin. Sampling from large matrices: An approach through geometric functional analysis. Journal of the
ACM (JACM), 54(4):21, 2007.
[7] Petros Drineas, Malik Magdon-Ismail, Michael W Mahoney, and David P Woodruff. Fast approximation of matrix coherence and statistical
leverage. Journal of Machine Learning Research, 13(Dec):3475–3506, 2012.
[8] Miao Xu, Rong Jin, and Zhi-Hua Zhou. Cur algorithm for partially observed matrices. In Proceedings of the International Conference on
Machine Learning (ICML), pages 1412–1421, 2015.
[9] Isabelle Guyon, Steve R Gunn, Asa Ben-Hur, and Gideon Dror. Result analysis of the nips 2003 feature selection challenge. In NIPS,
volume 4, pages 545–552, 2004.

14

[10] Petros Drineas, Michael W Mahoney, and S Muthukrishnan. Sampling algorithms for l 2 regression and applications. In Proceedings of
the seventeenth annual ACM-SIAM symposium on Discrete algorithm, pages 1127–1136. Society for Industrial and Applied Mathematics,
2006.
[11] Shusen Wang and Zhihua Zhang. A scalable cur matrix decomposition algorithm: lower time complexity and tighter bound. In Advances
in Neural Information Processing Systems, pages 647–655, 2012.
[12] Lester W Mackey, Michael I Jordan, and Ameet Talwalkar. Divide-and-conquer matrix factorization. In Advances in Neural Information
Processing Systems, pages 1134–1142, 2011.
[13] Christos Boutsidis and David P Woodruff. Optimal cur matrix decompositions. In Proceedings of the 46th Annual ACM Symposium on
Theory of Computing, pages 353–362. ACM, 2014.
[14] Christos Boutsidis, Petros Drineas, and Malik Magdon-Ismail. Near-optimal column-based matrix reconstruction. SIAM Journal on
Computing, 43(2):687–717, 2014.
[15] Yining Wang and Aarti Singh. Provably correct active sampling algorithms for matrix column subset selection with missing data. arXiv
preprint arXiv:1505.04343, 2015.
[16] Alex Gittens, Jey Kottalam, Jiyan Yang, Michael F. Ringenburg, Jatin Chhugani, Evan Racahk, Mohitdeep Singh, Yushu Yao, Curt Fischer,
Oliver Ruebel, Benjamin Bowen, Norman G. Lewis, Michael W. Mahoney, Venkat Krishnamurthy, and Prabhat. A multi-platform evaluation
of the randomized CX low-rank matrix factorization in Spark. In International Workshop on Parallel and Distributed Computing for Large
Scale Machine Learning and Big Data Analytics (ParLearning), 2016.

A PPENDIX
A. Proof of Lemma 2
Proof. Note that,
"
E

A(jt ) B (jt )
gpjt

!

#
=

i1 i2

G
X

pk

k=1

A(k) B (k)
gpk

!
=
i1 i2

1
(AB)i1 i2
g

Since each block is picked independently we have,
! #
" g
X A(jt ) B (jt )
var[(CR)i1 i2 ] = var
gpjt
t=1
i i
! 1 2#
"
g
(jt )
X
A B (jt )
=
var
gpjt
t=1
i1 i2
 
"
!2 
! #2 
g
(jt )
(j
)
t
X
A B (jt )
A B (jt )
E 
−E

=
gp
gpjt
jt
t=1
i1 i2
i1 i2


!2
G
X
(AB)2i1 i2
A(k) B (k)

pk
= g
−
gpk
g2
k=1


=

1

g

i1 i2


2
G
A(k) B (k)
X

i1 i2

k=1

pk



− (AB)2i1 i2 

15

E[kAB − CRk2F ] =

p
m X
X

var[(CR)i1 i2 ]

i1 =1 i2 =1




2
p
G
m X
A(k) B (k)
X
X
1

i1 i2
− (AB)2i1 i2 
=

g
p
k
i =1 i =1
1

=

p
G
m
2
X
1 X X  (k)
A B (k)
gpk i =1 i =1
i1 i2

k=1

=

≤

where αA = mink



kA(k) kF
kA(k) k2

2

!
−

kABk2F
g

2

gpk

−

kABk2F
g

G
X
kA(k) k22 kB (k) k2F
k=1

≤

1

G
X
kA(k) B (k) k2F
k=1

≤

k=1

2

gpk

G
X

kA(k) k2

k=1

kA(k) kF

!2

kA(k) k2F kB (k) k2F
gpk

1
kAk2F kBk2F
βαA g

. Also, note 1 ≤ αA ≤ s.

By Jensen’s inequality,
E[kAB − CRkF ]2 ≤ E[kAB − CRk2F ]
≤

1
kAk2F kBk2F
βαA g

And by Markov’s inequality, with probability ≥ 1 − δ, we have
kAB − CRkF ≤

1
1
E[kAB − CRkF ] ≤ √
kAkF kBkF
δ
δ βαA g

B. Proof of Lemma 1
First, note U = (RS)† and C = AS.
kA − CU RkF = kA − AS(RS)† RkF
Recall that R ∈ Rr×n has rank no greater than r; A ∈ Rm×n ; ε ∈ (0, 1); and that the same column blocks from

16

R and A are picked with the following probability distribution:
kV TR,r E i k2F
, ∀i ∈ [G].
r

pi =

We can use Lemma 1.1 (stated and proved in next section) with probability at least 0.85 we have
kA − AS(RS)† RkF ≤ (1 + ε)kA − AR† RkF .
In the paper we collapse Lemma 1.1 into the proof of Lemma 1, but here we state and prove the Lemma separately
since it may be of general interest.
Next, we bound kA − AR† RkF . Since A has incoherent column space, the uniform sampling distribution
pj = 1/m satisfies eqn. (13) in [1] with β = 1/µ0 . Consequently, we can apply modified version of Theorem 1 in
[1] we get with probability at least 0.85, kA − AR† RkF ≤ (1 + ε)kA − Ak kF . Finally, we get with probability
0.7
kA − CU RkF ≤ (1 + ε0 )2 kA − Ak kF
≤ (1 + ε00 )kA − Ak kF

letting ε00 = 3ε0

This completes the proof of Lemma 1.
1) Approximating generalized `2 regression in the block setting:
In this section, we give theory for generalized least squares using block subset selection that is used to prove
the main results for the algorithms but applies to arbitrary matrices A and B. Given matrices A ∈ Rm×n and
B ∈ Rr×n , the generalized least squares problem is
min

X∈Rm×r

kA − XBkF .

c = AB † . To approximate this
It is well-known that the solution to this optimization problem is given by X
problem by a subsampled problem, we sample some blocks of columns from A and B to approximate the standard
`2 regression by the following optimization:
min

X∈Rm×r

k(AS) − X(BS)kF .

The solution of this problem is given by X̃ = AS(BS)† . In the following lemma, we give a guarantee stating
that, when enough blocks are sampled with the specified probability, the approximate solution is close to the actual
solution to the `2 regression.
Lemma 3. Suppose B ∈ Rr×n has rank no greater than k; A ∈ Rm×n ; ε, δ ∈ (0, 1); and let the same column
blocks from B and A be picked with the following probability distribution:
pi =

k(V B,k )(i) k2F
, ∀i ∈ [G].
k

17

2

If g = O( αBkδ4 ε2 ) blocks are chosen, then with probability at least 1 − δ we have
kA − AS(BS)† BkF ≤ (1 + ε)kA − AB † BkF .
Proof. Let B = U k Σk V Tk and α = maxi



kV T
k Ei k 2
kV T
k Ei k F

2

−1 2 −2 −2
We start by showing V Tk S is full rank. Using Lemma 2, if g ≥ 8αR
k δ ε1 and 0 < ε1 < 1, we get the

following with probability ≥ 1 − δ1
k
ε1
kV Tk V k − V Tk SS T V k k2 = kI k − V Tk SS T V k k2 ≤ 4 √
≤
δ αR g
2
This further gives us a bound on the singular values of V Tk S, for all i,
|1 − σi2 (V Tk S)| = |σi (V Tk V k ) − σi (V Tk SS T V k )| ≤ kI k − V Tk SS T V k k2 ≤ ε1

(5)

Thus, it follows for all singular values of V Tk S
√

1 − ε1 ≤ σi (V Tk S) ≤

√

1 + ε1

(6)

Now, consider
kΩk2 = k(V Tk S)† − (V Tk S)T k2

by definition

= kΣ−1
− ΣV Tk S k2
V TS

from SVD

k

= max σi (V Tk S) −
i

= max
i

≤q
≤p

1
σi (V Tk S)

by definition

|σi2 (V Tk S) − 1|
|σi (V Tk S)|

by simple manipulation

kV Tk V k − V Tk SS T V k k2

from equation 5

1 − kV Tk V k − V Tk SS T V k k2
ε1 /2

from lemma 2 w.h.p.

1 − ε1 /2
√
≤ ε1 / 2

if ε1 < 1 =>

Also, for any Q we have ,
" g
X
1
E[kQSk2F ] = E
Q(jt )
√
gp
j
t
t=1

2
F

#


g
X
1
Q(jt )
=
E
gp
j
t
t=1

2


=

F

p

g X
G
X
t=1 i=1

√
1 − ε1 /2 > 1/ 2

pi

1
Q(i)
gpi

2
F

= kQk2F

By Jensen’s inequality,
E[kQSkF ]2 ≤ E[kQSk2F ] = kQk2F
By applying Markov’s inequality, we get with probability ≥ 1 − δ 0 ,
kQSkF ≤

1
1
E[kQSkF ] ≤ 0 kQkF
δ0
δ

(7)

18

The following will be useful later,
T
T
T
† T
AS(BS)† B = AS(U k Σk V Tk S)† U k Σk V Tk = AS(V Tk S)† Σ−1
k U k U k Σk V k = AS(V k S) V k

We break down the left hand term into 3 manageable components:
kA − AS(BS)† BkF
= kA − AS(V Tk S)† V Tk kF

by Lemma (shown above)

⊥T
T
† T
= kA − AV k V Tk S(V Tk S)† V Tk + AV ⊥
k V k S(V k S) V k kF

⊥T
since (V k V Tk + V ⊥
kVk )=I

⊥T
T
† T
= kA − AV k V Tk + AV ⊥
k V k S(V k S) V k kF

since V Tk S is full rank

⊥T
T
† T
≤ kA − AV k V Tk kF + kAV ⊥
k V k S(V k S) V k kF

triangle inequality

⊥T
⊥ ⊥T
T
T
= kAV ⊥
k V k kF + kAV k V k S(Ω + (V k S) )kF

define Ω = (V Tk S)† − (V Tk S)T

⊥T
⊥ ⊥T
⊥ ⊥T
T
≤ kAV ⊥
k V k kF + kAV k V k SkF kΩk2 + kAV k V k SS V k kF


1
⊥T
⊥ ⊥T
⊥ ⊥T
T
≤ 1 + 0 kΩk2 kAV ⊥
k V k kF + kAV k V k V k − AV k V k SS V k kF
δ


1
1
⊥T
⊥T
kAV ⊥
≤ 1 + 0 kΩk2 kAV ⊥
√
k V k kF kV k kF
k V k kF +
δ
δ2 αR g
√


1
k
⊥T
⊥ ⊥T
kAV ⊥
≤ 1 + 0 kΩk2 kAV k V k kF + √
k V k kF
δ
δ2 αR g


1
ε1
⊥T
≤ 1 + 0 kΩk2 + √
kAV ⊥
k V k kF
δ
8δ2


1
ε1
≤ 1 + 0 kΩk2 + √
kA − AB † BkF
δ
8δ2

where (*) follows since

√
√ k
αR g

≤

δ1

√k
αR g

≤

submultiplicativity
by (7) and since V ⊥T
k Vk =0
Lemma 2 andkV k kF =

see (*) below
since AV k V Tk = AB † B

ε1
√
8

≤ (1 + ε) kA − AB † BkF
6ε1
δ .
2

Lemma 2 is used, then g ≥ 36 ∗ 8αR εk2 δ4 . Finally, note that ε1 ≤ ε < 1 by assumption.

C. Proof of Corollary 1
First, note U = (RS)† and C = AS.
kA − CU RkF = kA − AS(RS)† RkF

k

⊥T
kA − AB † BkF = kAV ⊥
k V k kF

Thus, we can conclude the following with probability ≥ 1 − (δ 0 + δ1 + δ2 ) = 1 − δ


 
1
1
kA − AS(BS)† BkF ≤ 1 + √
+
ε1 kA − AB † BkF
2δ2
2δ 0

by setting δ 0 = δ1 = δ2 = δ/3 and ε =

√

19

Similar to the proof of Lemma 1, we can use Lemma 1.1 with probability at least 0.85 we have
kA − AS(RS)† RkF ≤ (1 + ε)kA − AR† RkF .
Recall that A ∈ Rm×n ; ε ∈ (0, 1); and that the rows R are picked from A with the following probability
distribution:
pi =

keTi U A,k k22
, ∀i ∈ [m].
k

We bound kA − AR† RkF using Theorem 1 in [1] we get with probability at least 0.85, kA − AR† RkF ≤
(1 + ε)kA − Ak kF . Finally, we get with probability 0.7
kA − CU RkF ≤ (1 + ε0 )2 kA − Ak kF
≤ (1 + ε00 )kA − Ak kF
This completes the proof of Corollary 1.

letting ε00 = 3ε0

