EFFICIENT RANDOMIZED ALGORITHMS FOR THE
FIXED-PRECISION LOW-RANK MATRIX APPROXIMATION

arXiv:1606.09402v3 [math.NA] 10 Feb 2018

WENJIAN YU∗ , YU GU† , AND YAOHANG LI‡
Abstract. Randomized algorithms for low-rank matrix approximation are investigated, with the
emphasis on the fixed-precision problem and computational efficiency for handling large matrices.
The algorithms are based on the so-called QB factorization, where Q is an orthonormal matrix.
Firstly, a mechanism for calculating the approximation error in Frobenius norm is proposed, which
enables efficient adaptive rank determination for large and/or sparse matrix. It can be combined with
any QB-form factorization algorithm in which B’s rows are incrementally generated. Based on the
blocked randQB algorithm by P.-G. Martinsson and S. Voronin, this results in an algorithm called
randQB EI. Then, we further revise the algorithm to obtain a pass-efficient algorithm, randQB FP,
which is mathematically equivalent to the existing randQB algorithms and also suitable for the
fixed-precision problem. Especially, randQB FP can serve as a single-pass algorithm for calculating
leading singular values, under certain condition. With large and/or sparse test matrices, we have
empirically validated the merits of the proposed techniques, which exhibit remarkable speedup and
memory saving over the blocked randQB algorithm. We have also demonstrated that the single-pass
algorithm derived by randQB FP is much more accurate than an existing single-pass algorithm. And
with data from a scenic image and an information retrieval application, we have shown the advantages
of the proposed algorithms over the adaptive range finder algorithm for solving the fixed-precision
problem.
Key words. adaptive rank determination, randomized algorithm, low-rank matrix approximation, pass-efficient algorithm, fixed-precision problem.
AMS subject classifications. 15A18, 65F30, 65F15, 68W20, 60B20

1. Introduction. Low-rank matrix factorizations, like the partial singular value
decomposition (SVD) and the rank-revealing QR factorization, play a crucial role in
data analysis and scientific computing. In recent years, techniques based on randomization have been investigated for performing the computation and low-rank factorization of large matrices [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. They involve the same or fewer
floating-point operations (f lops) than classical algorithms, and are more efficient by
exploiting modern computing architectures.
A basic idea of the randomized techniques is using random projection to approximate the dominant subspace of a matrix. For an m × n matrix A, suppose the
orthogonal basis vectors of this approximate subspace form an m × k orthonormal
matrix Q. Then, we have [1, 2]:
(1)

A ≈ QB,

where B is a k × n matrix, and
(2)

B = QT A.

Standard factorizations, e.g., SVD, can be further performed on the smaller matrix
B, to obtain the low-rank factorizations of A.
∗ TNList, Department of Computer Science and Technology, Tsinghua University, Beijing 100084,
China (yu-wj@tsinghua.edu.cn).
† Department of Computer Science and Technology and Institute for Interdisciplinary Information
Sciences, Tsinghua University, Beijing 100084, China (guyu13@mails.tsinghua.edu.cn).
‡ Department of Computer Science, Old Dominion University, Norfolk, VA 23529, USA
(yaohang@cs.odu.edu).

1

The basic randQB algorithm
Input: A, k, s.
Output: Q, B.
(1) Ω= randn(n, k + s)
(2) Q = orth(AΩ)
(3) B = QT A

The blocked randQB algorithm
Input: A, ε, b.
Output: Q, B .
(1) for i = 1, 2, 3, · · ·
(2) Ωi = randn(n, b)
(3) Qi = orth(AΩi )

T
(4) Qi = orth(Qi − i−1
j=1 Qj Qj Qi )
T
(5) B i = Qi A
(6) A = A − Qi B i
(7) if A < ε then stop
(8) end for
(9) Q = [Q1 · · · Qi ]; B = [B T1 · · · B Ti ]T .

(a) The randQB algorithm

(b) The randQB b algorithm

Fig. 1. The randomized algorithms for the QB factorization.

The approximation presented by (1) and (2) can also be regarded as a kind of lowrank factorization of A, called QB factorization or QB approximation in this work. In
[2], a basic randomized scheme for computing the QB approximation was presented, as
shown in Figure 1(a). For producing close to optimal rank-k approximation, the oversampling scheme using a random Gaussian matrix Ω with k + s columns is employed,
where s is a small integer. We use randQB to denote this algorithm.
Usually, the problem of low-rank matrix approximation falls into two categories:
• The fixed-rank problem, where the rank parameter k is given.
• The fixed-precision problem, where we seek Q and B with as small as possible
size such that kA − QBk < ε, where ε is a given accuracy tolerance.
A blocked variant of the randQB algorithm proposed in [1], i.e. the randQB b
algorithm in Figure 1(b), is suitable for the fixed-precision problem. It incrementally
builds the factors Q and B based on the combination of the randQB algorithm
and the blocked Gram-Schmidt scheme, and measures the approximation error by
explicitly maintaining the residual matrix. However, it is inefficient or even fails for
handling large matrix, because maintaining the residual matrix is costly in runtime
and memory usage.
In this work, the randomized algorithms for the fixed-precision problem are investigated considering their adaptability to large and/or sparse matrices. Firstly, a
mechanism is proposed for calculating the error of QB approximation in Frobenius
norm during the iterative process of building Q and B. It does not require maintaining the residual matrix (or updating matrix A), and thus avoids fill-in while handling a
sparse A. This mechanism is also applicable to other iterative computing procedures,
e.g., the trQRCP algorithm [14]. Secondly, the algorithm is further revised to largely
reduce the number of passes over matrix A, in order to adapt the scenarios where the
cost of accessing matrix A is expensive. These techniques result in two algorithms
called randQB EI and randQB FP, which inherit the merits of randQB/randQB b
algorithms and have extra benefits. Numerical experiments are carried out on a multicore computer to validate the efficiency and accuracy of the proposed algorithms for
handling large or sparse matrices in practical scenarios. The results show that our
randQB EI and randQB FP algorithms have up to 3X speedup and 3X memory saving over an implementation of the randQB b algorithm for dense matrices. They
2

also exhibit up to 22X speedup over the implementation of randQB b algorithm for
sparse matrices. Compared with the single-pass algorithm and the adaptive randomized range finder in [2], the proposed algorithms exhibit much better accuracy
or avoid the large overestimation of the sizes of Q and B. For reproducibility, we
have shared the Matlab codes of the proposed algorithms and experimental data on
https://github.com/WenjianYu/randQB auto.
2. Technical preliminaries. This section summarizes the background we need
for presenting the proposed techniques. Throughout the paper, we measure vectors
with their Euclidean norm. Two kinds of matrix norm are usually considered: Frobenius norm and spectral norm (l2 -norm). The spectral norm of a matrix is relatively
difficult to calculate, though
P it is often more informative for noisy data [24]. The
Frobenius norm kAkF = ( i,j |A(i, j)|2 )1/2 is easier for calculation, and thus more
widely used in data analysis and machine learning applications [10]. We measure
matrices with their Frobenius norm by default. We also assume that all matrices are
real valued, although the generalization to complex matrices is of no difficulty.
2.1. Randomized algorithms. To produce a rank-k factorization for an m × n
matrix A, with the basic randQB algorithm in Figure 1(a) one obtains an m × l
orthonormal matrix Q and an l × n matrix B, where l > k due to over-sampling.
With this QB approximation, the standard factorizations can be efficiently computed.
For example, the standard SVD algorithm can be performed on B, which results in
T
B = Ũ Σ̃Ṽ . Then,
(3)

T

A ≈ QB = QŨ Σ̃Ṽ .

The first k columns of matrices QŨ and Ṽ and the k × k upper-left submatrix of Σ̃
approximate the rank-k SVD factors of A. Similarly, by changing the factorizations
made on B, one obtains the approximate QR factorization and CUR factorization,
etc [1, 4]. Notice that the accurate truncated SVD provides the optimal low-rank
approximation [15]. However, computing accurate SVD of a large matrix is costly, and
in many applications the optimality is not necessary. It is thus acceptable to compute
the approximate low-rank factorizations for gains in computational efficiency.
The approximation error of the randomized algorithm is a random variable. The
authors of [2] have studied the properties of the error in term of spectral and Frobenius
norms, and given the bounds on their expectation and variance.
The truncated QR factorization with column pivoting can also be used for lowrank matrix approximation [11, 12]. Classical pivoted QR factorization has the disadvantage that it is hard to be parallelized or to take usage of BLAS-3 operation.
Recently, the pivoted QR factorization was largely accelerated through utilizing a
randomized technique, which achieves the efficiency comparable to the unpivoted QR
decomposition [13, 14]. This makes the truncated pivoted QR factorization competitive for low-rank approximation. Notice that QR factorization can be regarded as a
special case of the QB factorization. Therefore, one of the techniques proposed here
(c.f. Sect. 3) could benefit the solution of the fixed-precision problem based on QR
factorizations, as well.
The major computation of the randQB algorithm lies at the multiplication of
A and Ω. If A is sparse or a structured matrix (often implicitly defined) for which
matrix-vector products can be rapidly evaluated, the cost of multiplication can be
largely reduced (even to O(m + n) flops). The implicitly-defined structured matrix
often arises from physical problems, such as a discretized integral operator applied via
3

the fast multipole method, and is sometimes referred to as an implicit sparse matrix.
For the fixed-precision problem, an adaptive randomized range finder was proposed in [2]. It employs the incremental sampling approach with a probabilistic error
estimator to determine the size of Q and B. It is based on the statement that
r
2
T
(4)
kA − QQ Ak2 ≤ 10
max k(A − QQT A)ω (i) k,
π i=1,··· ,r
with probability at least 1 − 10−r [2]. Here k · k2 stands for the spectral norm, ω (i) is
a random vector, and r is a small integer, e.g., r = 10. However, the error estimator
often overestimates the approximation error, yielding much larger output matrices
than what is necessary.
The randQB b algorithm in [1] is based on the single-vector version of randQB
algorithm with Gram-Schmidt procedure, which allows to construct the QB factorization and to evaluate its error step by step. In order to exploit blocking to attain high
performance of linear algebraic computation, the algorithm is then converted to the
blocked randQB algorithm in Figure 1(b), where step (6) is for calculating the residual
matrix. Therefore, the randQB b algorithm allows a precise error calculation for the
fixed-precision problem. Notice that step (4) there is a re-orthogonalization operation
which eases the accumulation of round-off error under floating point arithmetic.
The blocked randQB algorithm is able to produce Q and B in smaller sizes than
the adaptive randomized range finder. However, explicitly maintaining the residual
matrix brings extra time and memory cost if a large matrix A is handled. This
disadvantage becomes more serious if A is also sparse, because the fill-in phenomena
leads to a dense residual matrix.
In order to reveal the difference among the relevant randomized algorithms, we
give a brief comparison of them for the fixed-rank problem, presented as Table 1.
svds denotes the Matlab built-in command for truncated SVD [22], which is based
on a Krylov subspace iterative method. trQRCP denotes the randomized pivoted QR
factorization [14]. randQB FP is one of the contributions in this paper, which is a
pass-efficient algorithm (c.f. Sect. 4). To depict the performance of the algorithms
in the situation where the cost of accessing matrix entries is expensive (e.g. A are
stored in slow memory) [2], we include the number of passes over matrix A in Table
1. From the table, we see that the randQB FP algorithm inherits the merits of the
randQB algorithm, and is also suitable for the fixed-precision problem.
Table 1
Properties of the randomized algorithms for the fixed-rank matrix approximation (P is a small
integer for power scheme, and b is the block size).

svds
Computational efficiency
low
Adaptive rank determination no
Number of passes over matrix αk *
*
†
‡

randQB randQB b trQRCP randQB FP
high
high
high
high
no
yes
yes†
yes
2 or 2+2P
4k/b
3k/b + 1 1 or 1+2P ‡

α is a number larger than 1.
This was not mentioned in [14]. For more detail, please see Remark 3.2.
Suppose matrix A is stored in the row-major format.

3. An efficient Frobenius-norm error indicator and its application. In
this section, we first propose an error indicator for measuring the approximation error
in Frobenius norm and an efficient framework for solving the fixed-precision problem
of QB factorization. Then, the randQB EI algorithm is derived based on the blocked
4

randQB algorithm. Finally, the accuracy and validity of the error indicator in floating
point arithmetic is discussed.
3.1. An error indicator. We first give a theorem regarding Frobenius norm of
the error of the QB approximation.
theorem 1. Let A be an m × n matrix. Q denotes an m × k orthonormal matrix
(k < m), and B = QT A. Then,
(5)

kA − QBk2F = kAk2F − kBk2F .
Proof. Due to the property of Frobenius norm, we know for any matrix M ,

(6)

kM k2F = tr(M T M ),

where tr(·) calculates the trace of a matrix. Because Q is orthonormal and B = QT A,
(A − QB)T (A − QB) = (A − QQT A)T (A − QQT A)
(7)

= AT A − 2AT QQT A + AT QQT QQT A
= AT A − AT QQT A
= AT A − B T B.

Now, applying the trace operation to both sides of (7), and according to (6), we obtain
(5).
Theorem 1 suggests that, if we have access to B and kAkF is known a prior,
we can calculate the error of QB approximation without referring to the residual
matrix. This leads to a framework for solving the fixed-precision problem, presented
as Algorithm 1. It suits to any algorithm that incrementally generates the rows of
B, including randQB b, trQRCP and the randQB FP algorithm presented in Sect.
4. Below we prove the correctness of Algorithm 1.
Algorithm 1 A framework for solving the fixed-precision QB factorization problem
Input: an m × n matrix A; desired accuracy tolerance ε.
Output: Q, B, such that kA − QBkF < ε.
1: Q = [ ]; B = [ ];
# empty matrices
2: E = kAk2F
# initialization of the error indicator
3: for i = 1, 2, 3, · · · do
4:
Generate Qi and B i , s.t. [Q, Qi ] is orthonormal and B i = QTi A.
5:
Q = [Q,
 Qi]
B
6:
B=
Bi
7:
E = E − kB i k2F
# update the error indicator
8:
if E < ε2 then stop
9: end for
theorem 2. After the i-th iteration of loop in Algorithm 1 is executed,
(8)

E (i) = kA − Q(i) B (i) k2F ,

where E (i) , Q(i) and B (i) denote the values of E, Q and B after the i-th iteration of
the loop is executed, respectively.
5

Proof. Based on step 7 of Algorithm 1 and the property of Frobenius norm,
(9)

E (i) = kAk2F −

i
X
j=1

kB j k2F = kAk2F − kB (i) k2F .

Because Q(i) is an orthonormal matrix and B (i) = Q(i)T A, (8) can be obtained by
applying Theorem 1.
From Theorem 2, we see that E in Algorithm 1 equals to the square of the approximation error. It is an error indicator updated through calculating the Frobenius
norm of B i . This yields two benefits: we no longer need to maintain the residual
A − QB, and the approximation error can be calculated with very small cost.
Steps 7 and 8 in Algorithm 1 can be replaced by a row-by-row calculation scheme.
7a: for j = 1, 2, · · · mi do
# mi denotes the number of rows of B i
7b: E = E − kB i (j, :)k2F
8a: if E < ε2 then
8b:
remove B i (j + 1 : mi , :) from B; remove Qi (:, j + 1 : mi ) from Q
8c:
stop
8d: end if
8e: end for
This adds negligible cost, but allows us to determine the certain row of B i where the
accuracy tolerance is just attained. It makes the column (row) number of outputted Q
(B) an arbitrary integer, instead of a multiple of block size b in randQB b algorithm.
Remark 3.1. If A is an implicit sparse matrix, the proposed framework needs
more effort for calculating its Frobenius norm. The columns of A can be solved by
multiplying A with the canonical basis vectors. Therefore, the cost for calculating
kAkF will be O(n(m + n)) flops, if each matrix-vector product costs O(m + n) flops.
This might be affordable, as it is executed just once.
3.2. The randQB EI algorithm. The combination of Algorithm 1 and the
randQB b algorithm results in Algorithm 2 (called randQB EI), whose steps 4∼7
replace step 4 in Algorithm 1. Notice that step 7 in Algorithm 2 looks the same as
step (5) in randQB b algorithm, but is actually different. And, step (3) in randQB b
algorithm becomes step 5 in Algorithm 2, which is the blocked Gram-Schmidt orthogonalization of AΩi (notice B = QT A).
Due to Theorem 2 and the orthogonality of Q, we have the following proposition.
proposition 1. The randQB EI algorithm (Algorithm 2) is equivalent to the
randQB b algorithm, when executed in exact arithmetic.
Assuming that multiplying two dense matrices of sizes m × n and n × l costs
Cmm mnl flops, and performing an economic QR factorization of an m×n dense matrix
costs Cqr mn min(m, n) flops, we can analyze the flop counts of the relevant algorithms
and compare their performance for handling a dense A. If TrandQB and TrandQB b
denote the runtime of the algorithms randQB and randQB b respectively [1],
(10)
(11)

TrandQB ∼ 2Cmm mnl + Cqr ml2 ,
TrandQB

b

2
∼ 3Cmm mnl + Cmm ml2 + Cqr ml2 ,
t

where l is the number of columns in the resulting matrix Q, and t satisfies l = tb.
Note that b is practically much smaller than l (say, b =10 or 20), although the optimal
choice of block size depends strongly on what hardware is used [1].
6

Algorithm 2 The randQB EI algorithm for the fixed-precision problem
Input: an m × n matrix A; desired accuracy tolerance ε; block size b.
Output: Q, B, such that kA − QBkF < ε.
1: Q = [ ]; B = [ ];
2: E = kAk2F
3: for i = 1, 2, 3, · · · do
4:
Ωi = randn(n, b)
5:
Qi = orth(AΩi − Q(BΩi ))
6:
Qi = orth(Qi − Q(QT Qi ))
# re-orthogonalization
7:
B i = QTi A
# no need to calculate QTi (A − QB)
8:
Q = [Q,
 Qi]
B
9:
B=
Bi
10:
E = E − kB i k2F
11:
if E < ε2 then stop
12: end for

For the randQB EI algorithm, the runtime can be similarly depicted:
(12)

TrandQB

EI

2
∼ 2Cmm mnl + Cmm (2m + n)l2 + Cqr ml2 .
t

Because l is usually much smaller than m and n, we see that the flop count of
randQB EI is about 2/3 of that of randQB b, and is comparable to that of the
basic randQB algorithm. Notice that Cqr is several times larger than Cmm . So, the
flop count of randQB EI could be smaller than that of randQB in the situation where
t is a large number. If A is sparse, this would more possibly happen, because the
randQB algorithm loses the benefit brought by the BLAS-3 operation.
The advantage of randQB EI over randQB b becomes more prominent if A is
a sparse matrix. With the proposed error indicator, we no longer need the residual
matrix. In contrast, it is always a dense matrix in the randQB b algorithm, and costs
much larger memory and induces much more computations.
Remark 3.2. Algorithm 1 can also be combined with the trQRCP algorithm
[14]. Although trQRCP generates a column permutation matrix as well, it does not
affect the Frobenius norm of each partial or the entire matrix of B. Therefore, this
will produce another efficient algorithm for adaptive low-rank matrix approximation,
which also adapts to sparse matrices.
3.3. Floating point arithmetic. Below we discuss the accuracy of the error
indicator in floating point arithmetic. We use E(·) and Er (·) to denote the functions
of error and relative error, respectively.
As the error indicator E = kAk2F − kBk2F , its calculated value Ê cannot be
accurate when E is very small, due to the cancellation in calculation. In floatingpoint arithmetic, the machine precision mach characterizes the maximum relative
error of converting a real number to its floating-point representation, i.e.
(13)

∀x, |Er (x)| ≤ mach .

According to the definition of Frobenius norm, kAk2F is the summation of squares of
matrix entries. Therefore, the relative error of kAk2F is bounded by 2mach . The same
7

thing applies to kBk2F . So,
|E(E)| = |E(kAk2F − kBk2F )| ≤ |E(kAk2F )| + |E(kBk2F )|
≤ 2mach (kAk2F + kBk2F )

(14)

< 4mach kAk2F .

If we want to guarantee that E has a relative error no more than δ, i.e., |E(E)| ≤ δE,
we shall enforce
4mach kAk2F ≤ δE.

(15)

√
This means the preset accuracy tolerance ε, which is larger than E at the termination of the algorithm, should satisfy:
r
r
√
4mach kAk2F
4mach
(16)
ε> E≥
=
kAkF .
δ
δ
So, we obtain the following Theorem.
theorem 3. Suppose matrix
q A and accuracy tolerance ε are the input to the

randQB EI algorithm. If ε >

4mach
kAkF ,
δ

the relative error of the calculated error

indicator E must be no more than δ. E.g., if ε > 2.1 × 10−7 kAkF , the error of E is
within 1% in the double-precision floating arithmetic, where mach ≈ 1.11 × 10−16 .
Notice that, as an error indicator for the fixed-precision problem, E should have
sufficient accuracy (e.g., with relative error ∼ 1% or less). Otherwise, the outputted
QB factorization would not satisfy the preset accuracy tolerance.
Besides, the orthogonality of Q also affects the accuracy of E. As the number
of columns in Q increases, its orthogonality gradually degrades. This issue occurs
for Q produced either by a single run of QR factorization (based on Householder
transformation) or by a Gram-Schmidt procedure followed by the re-orthogonalization
step. We will investigate its effect in the following experiment.
An n × n matrix A is constructed to have singular values according to a decaying
exponential. Two instances are tested, with singular value σj = e−j/20 and σj =
e−j/200 , j = 1, 2, · · · , respectively. The results obtained from executing randQB EI
algorithm (b = 10) and randQB algorithm with different values of rank parameter l
are shown in Figure 2. Note that for some large value of l, the error indicator can
be of negative value, such that it cannot be drawn in the log-scale plot. From the
figure, we can validate the correctness of Theorem 3. Providing that the square of
error kA − QBk2F /kAk2F > (2.1 × 10−7 )2 = 4.4 × 10−14 , the error indicator matches
the square of error very well. This holds even when l is larger than 3, 000, which
corresponds to the situation with larger accumulated round-off error. In Figure 2, the
value of kQT Q − Ik∞ 1 is also plotted, which reveals the loss of orthogonality of Q.
The results show that this issue is not severe in both the randQB and randQB EI
algorithms, although it gradually increases as the columns of Q are increased.
Remark 3.3. Theorem 3 suggests the limitation of the error indicator and the
proposed algorithms for the fixed-precision problem. It means that the efficient framework for adaptive rank determination would not work, in double-precision floating
arithmetic, for the problem with the accuracy tolerance ε less than 2.1 × 10−7 kAkF .
1 This

measure of loss of orthogonality follows Cleve Moler’s blog with title “Compare GramSchmidt and Householder Orthogonalization Algorithms” posted on Oct. 17, 2016.
8

error indicator(randQB)
error indicator(randQB_EI)

error indicator(randQB)
error indicator(randQB_EI)

10 -6

error 2 (randQB)
error (randQB_EI)
||Q*Q-I|| (randQB)
||Q*Q-I|| (randQB_EI)

10 -8

error 2 (randQB)

10 -12

2

error 2 (randQB_EI)
||Q*Q-I|| (randQB)
||Q*Q-I|| (randQB_EI)

10 -13

10 -10
10 -14

10 -12

10 -15

10 -14

10 -16

10 -16

10 -18
150

200

250

300

350

400

450

3200

3400

3600

l

3800

4000

l

(a) n = 3, 000, σj = e−j/20

(b) n = 20, 000, σj = e−j/200

Fig. 2. Normalized error indicator E/kAk2F , error square kA − QBk2F /kAk2F and the loss of
orthogonality of Q v.s. Q’s column number l in randQB algorithm and randQB EI algorithm.

4. A pass-efficient algorithm for the QB factorization. The technique in
last section efficiently solves the fixed-precision problem measured in Frobenius norm.
However, the randQB EI algorithm is not suitable for the scenarios where accessing
matrix entries is very expensive (e.g., A is too large and has to be stored on hard
disk), as it visits A for considerable times (see Table 1). In this section, we propose
a pass-efficient algorithm which largely reduce the number of passes over A.
4.1. The version without re-orthogonalization. We first consider the fixedrank problem where the rank parameter k is given. A preliminary pass-efficient algorithm (presented as Algorithm 3) can be derived from the randQB b or randQB EI
algorithm. The steps correspond to those of randQB EI in Algorithm 2, one by
one, except that the applications of A are moved out of the loop, and step 6 for
re-orthogonalization is ignored. In Algorithm 3, step 6 is the same as step 4 of
randQB EI algorithm, and steps 7 and 8 correspond to step 5 of randQB EI. In step
9, R−T
is the inverse of the transpose of an upper triangular matrix Ri . This step
i
can be regarded as solving linear equations with the coefficient matrix RTi , which is
implemented by “\” operator in Matlab.
Because Y i = Qi Ri ,
Qi = Y i R−1
i .

(17)

Substituting it into step 7 of the randQB EI algorithm, we have:
T
B i = (Y i R−1
i ) A
T
= R−T
i Yi A

(18)

T
T T
T
= R−T
i (G(:, (i − 1)b + 1 : ib) − Ωi B Q )A

T
T
T T
T
= R−T
i ((A G(:, (i − 1)b + 1 : ib)) − Ωi B Q A)
T
T
T
= R−T
i (H(:, (i − 1)b + 1 : ib) − Ωi B B).

9

Algorithm 3 A preliminary pass-efficient algorithm for fixed-rank QB factorization
Input: A ∈ Rm×n , k, s, block size b.
Output: Q, B.
1: Q = [ ]; B = [ ]; l = k + s;
2: Ω = randn(n, l)
3: G = AΩ
T
4: H = A G
5: for i = 1, 2, 3, · · · , l/b do
6:
Ωi = Ω(:, (i − 1)b + 1 : ib)
7:
Y i = G(:, (i − 1)b + 1 : ib) − Q(BΩi )
8:
[Qi , Ri ] = qr(Y i )
T
T
T
9:
B i = R−T
i (H(:, (i − 1)b + 1 : ib) − Ωi B B)
10:
Q = [Q,
 Qi]
B
11:
B=
Bi
12: end for

This means that step 9 of Algorithm 3 is equivalent to step 7 of randQB EI. Therefore,
we obtain the following proposition.
proposition 2. Algorithm 3 is mathematically equivalent to the fixed-rank version of randQB EI algorithm without re-orthogonalization.
4.2. The version with re-orthogonalization. In reality, the loss of orthogonality among the columns of {Q1 , Q2 , · · · } occurs due to the accumulation of round-off
error. This means the re-orthogonalization step, i.e. step 6 in Algorithm 2, cannot be
ignored. Below we derive the revisions to Algorithm 3 to address this issue.
The re-orthogonalization step can be expressed as
Q̃i R̃i = Qi − QQT Qi ,

(19)

where Q̃i 6= Qi and R̃i 6= I due to the loss of orthogonality. Q̃i is better orthogonal
to the previously generated {Q1 , Q2 , · · · , Qi−1 } than Qi . Now, we need to derive a
formula for calculating B i which does not involve A explicitly. Based on (17),
−1

Q̃i = (I − QQT )Y i R−1
i R̃i .

(20)

So, the corresponding formula for B i is:
T

B̃ i =Q̃i A
(21)

=(R̃i Ri )−T Y Ti (I − QQT )A

=(R̃i Ri )−T (ΩTi AT − ΩTi B T QT )(A − QB),

where the formula of Y i and the equality B = QT A is taken usage of. The product
of the last two brackets can be further simplified.

(22)

(ΩTi AT − ΩTi B T QT )(A − QB)

= H Ti − GTi QB − ΩTi B T B + ΩTi B T QT QB
= H Ti − Y Ti QB − ΩTi B T B.
10

In the deduction, Gi and H i denote G(:, (i − 1)b + 1 : ib) and H(:, (i − 1)b + 1 : ib)
in Algorithm 3, respectively. Therefore,
(23)

B̃ i = (R̃i Ri )−T (H Ti − Y Ti QB − ΩTi B T B) .

Based on (19) and (23), we can derive the version with re-orthogonalization for
Algorithm 3. We just need to replace the 9th step with the following steps.
9a: [Qi , R̃i ] = qr(Qi − Q(QT Qi ))
9b: Ri = R̃i Ri
T
T
T
T
9c: B i = R−T
i (H(:, (i − 1)b + 1 : ib) − Y i QB − Ωi B B)

Notice that Qi and B i are overwritten to stand for Q̃i and B̃ i , respectively. Based
on Proposition 2 and the above deduction, we see that the pass-efficient algorithm
with re-orthogonalization is also mathematically equivalent to the fixed-rank versions
of randQB EI and randQB b algorithms.
This algorithm with fewer passes over A is called randQB FP. Based on the
notations in Section 3, its flop count analysis is as follows.
(24)

TrandQB

FP

2
∼ 2Cmm mnl + 2Cmm (m + n)l2 + Cqr ml2 ,
t

where t satisfies l = tb. Compared with the randQB EI algorithm, the randQB FP
algorithm has slightly larger flop count. However, while handling a dense A its actual
runtime may be shorter because it lumps the multiplications with A.
Remark 4.1. The round-off error may affect the accuracy of B i , and it increases
as the number of iterations increases. However, this may not be an issue for practical
low-rank approximation problems. In Section 5, we will present numerical experiments
to validate the effectiveness of the randQB FP algorithm, which shows it works very
well for many applications with the rank parameter up to several thousands or the
relative Frobenius-norm error of approximation as small as 10−7 .
Remark 4.2. The randQB FP algorithm can derive a single-pass algorithm, if
matrix A is stored in the row-major format or is revealed row(s) by row(s). Suppose
Ai,: denotes the i-th row of A. With it we have the i-th row of G, Gi,: = Ai,: Ω. Then,
P
because H = AT G = i ATi,: Gi,: , the i-th term in this summation can be obtained.
With all rows of A, in this way we can accomplish steps 3 and 4 in the randQB FP
algorithm with only one pass over A. It should be pointed out that this algorithm is
not a general single-pass algorithm, as it has the restriction of the matrix. For more
general single-pass algorithms for low-rank matrix approximation, please refer to the
recent work [23].
4.3. The inclusion of power iteration scheme. The error of randomized
QB factorization could be large for the matrix whose singular value decays slowly [2].
So, the power iteration scheme has been proposed to relieve this weakness [1, 2, 7].
Conceptually, the power iteration means replacing A with (AAT )P A, where P is
an integer. However, in floating-point computation any singular components smaller
1/(2P +1)
than σ1 mach
will be lost. This makes the orthonormalization steps after the
applications of A and AT necessary, and P should not be set to a large number.
Incorporating the power iteration, we have the randQB FP algorithm for the fixedprecision problem presented as Algorithm 4, where the error indicator E is utilized.
In Algorithm 4, a sufficient large value of ˜l should be set according to problemspecific experience and the concern of computing time. If the set ˜l is not large enough
for attaining the specified accuracy criterion, we need to re-generate the Ω matrix and
11

Algorithm 4 The randQB FP with power scheme for the fixed-precision problem
Input: A ∈ Rm×n , desired accuracy tolerance ε, block size b, power parameter P .
Output: Q, B, such that kA − QBkF ≤ ε.
1: Q = [ ]; B = [ ];
2: Ω = randn(n, ˜
l), where ˜l is a sufficiently large number.
3: for i = 1 : P do
4:
G = orth(AΩ)
5:
Ω = orth(AT G)
6: end for
7: G = AΩ
T
8: H = A G
9: E = kAk2F
10: for i = 1, 2, 3, · · · do
11:
Ωi = Ω(:, (i − 1)b + 1 : ib)
12:
Y i = G(:, (i − 1)b + 1 : ib) − Q(BΩi )
13:
[Qi , Ri ] = qr(Y i )
14:
[Qi , R̃i ] = qr(Qi − Q(QT Qi ))
15:
Ri = R̃i Ri
T
T
T
T
16:
B i = R−T
i (H(:, (i − 1)b + 1 : ib) − Y i QB − Ωi B B)
17:
Q = [Q,
 Qi]
B
18:
B=
Bi
2
19:
E=
√ E − kB i kF
20:
if E < ε then stop
21: end for
rerun the algorithm to collect additional columns/rows of Q and B. This situation
and the power scheme both increase the number of passes over A. But compared
to other algorithms for the fixed-precision problem, this fixed-precision randQB FP
algorithm involves much fewer passes over A.
5. Numerical results. In this section we compare the proposed algorithms
against several existing algorithms in terms of execution time, memory usage and
accuracy. All experiments are carried out on a Linux server with two 12-core Intel
Xeon E5-2630 CPUs @ 2.30 GHz, and 32GB RAM. For comparison of speed, the
proposed algorithms have been implemented in C based on the codes shared by the
authors of [1, 16]. The program is coded with OpenMP derivatives, and compiled
with the Intel ICC compiler with MKL libraries [17], to take full advantage of the
multi-core CPUs. The QR factorization and other basic linear algebra operations are
implemented through LAPACK routines which are automatically executed in parallel.
5.1. Comparison of speed. We compute the QB factorization of an n × n
matrix A. Notice the singular value distribution of matrix is immaterial for this
runtime comparison. Four algorithms are compared:
• The randQB algorithm in Figure 1(a);
• The randQB b algorithm in Figure 1(b), obtained from [16];
• The randQB EI algorithm presented in Section 3;
• The randQB FP algorithm presented in Section 4.
We compare their speed using both dense and sparse matrices, both as a function of
the dimension of the matrix and the parameter l denoting the number of the output
12

Q’s columns. The block size is b = 20 for the randQB b, randQB EI and randQB FP
algorithms. For each runtime measurement, the average time over 20 runs is reported.
Notice that the compared randQB b algorithm is an efficient parallel implementation
open-sourced in [16], also based on Intel MKL libraries.
In the first experiment we test the algorithms on dense matrices of varying size.
n ranges from 2,000 to 40,000. The value of l is always 200. The results are shown
in Figure 3 for the situations without and with the power scheme. The data of the
blocked randQB algorithm for the matrix with n = 40, 000 are not available due to
unreasonably long runtime of the program from [16]. From the results in Figure 3(a),
we see that the randQB EI and randQB FP algorithms are 2.4X (13.47s vs. 31.78s)
and 4.0X (8.01s vs. 31.78s) faster than the implementation of randQB b algorithm

1

1

10

Time(s)

Time(s)

10

0

10

randQB
randQB_EI
randQB_FP
randQB_b [1]

−1

10

2000

4000

10000
n

20000

0

10

randQB(P=1)
randQB_EI(P=1)
randQB_FP(P=1)
randQB_b(P=1) [1]

−1

10

40000

2000

(a) Without the power scheme

4000

10000
n

20000

40000

(b) With the power scheme

Fig. 3. Runtime of the randomized QB-factorization algorithms for dense matrices (l = 200).

respectively, when n = 32, 000. If the power scheme is imposed, the acceleration ratios
decrease to 1.8X and 3.0X respectively, which are still remarkable. The randQB
algorithm has the fastest computational speed, but its advantage over randQB FP
algorithm becomes marginal when the matrix size is large.
Here we only show the runtime results with the power parameter P = 1, as for
many applications this already achieves sufficient accuracy.
The memory costs for some large matrices are listed in Table 2. For the randQB,
randQB EI and randQB FP algorithms, the memory cost is mainly due to storing
matrix A. For the randQB b algorithm, it needs additional memory to store the
residual matrix and the product of QB. So, the proposed algorithms consume about
1/3 of that used by the blocked randQB algorithm. If we allow that A can be
overwritten by the residual matrix, the memory cost of the randQB b algorithm
[1] can be reduced, but still 2X larger than the proposed algorithms. For the largest
Table 2
The memory usage of the randomized QB-factorization algorithms for dense matrices (l = 200)

n

randQB

16,000
24,000
32,000
40,000

2,308 MB
4,792 MB
8,253 MB
12,694 MB

randQB EI randQB FP randQB b [1]
2,303 MB
4,796 MB
8,253 MB
12,581 MB
13

2,357 MB
4,873 MB
8,356 MB
12,714 MB

6,237 MB
13,618 MB
23,931 MB
N.A.

case with n = 40, 000, the randQB b algorithm actually requests more memory than
the size of RAM (∼ 32 GB), which explains the aforementioned long runtime of
randQB b.
The second experiment is about the algorithms’ efficiency for sparse matrices. We
generate sparse matrices with roughly 0.3% non-zero elements. They are stored in
CSR (compressed sparse row) format [18]. The runtimes of the algorithms are shown
in Figure 4. The results of the randQB b algorithm for the matrices with n ≥ 40, 000
are not available due to unreasonably long runtime. In contrast, it only takes a couple
of seconds for the other algorithms to process the largest matrix with n = 48, 000. We
see that the proposed algorithms take usage of the sparsity, while the blocked randQB
algorithm cannot. The speedup ratios of the former to the latter increase as the matrix size increases. For n = 32, 000, the randQB EI and randQB FP algorithms are
more than 22X and 14X faster than the implementation of randQB b algorithm, re-

1

1

10

Time(s)

Time(s)

10

0

10

randQB
randQB_EI
randQB_FP(CSR)
randQB_FP(COO)
randQB_b [1]

−1

10

−2

10
2000

4000

10000
n

20000

0

10

randQB(P=1)
randQB_EI(P=1)
randQB_FP(CSR,P=1)
randQB_FP(COO,P=1)
randQB_b(P=1) [1]

−1

10

−2

10
2000

50000

(a) Without the power scheme

4000

10000
n

20000

50000

(b) With the power scheme

Fig. 4. Runtime of the randomized QB-factorization algorithms for sparse matrices (l = 200).

spectively. Different from the situation for dense matrices, the randQB EI algorithm
becomes faster than randQB FP. This implies that lumping the multiplications of
a sparse matrix all together brings less benefit. And, randQB EI could run faster
than the randQB algorithm for matrix size over 48,000. This can be explained by
the comparison of (10) and (12), the inefficiency of orthogonalizing the whole matrix
of l columns, and that the sparse matrix removes the benefit of BLAS-3 operation to
randQB. Another interesting phenomenon is that if we instead store a large sparse
matrix with the COO (coordinate) format, the runtime of randQB FP algorithm can
be reduced by 30%. This means that the COO format is more adaptive to parallel
computing. If we set P = 1, similar observations regarding the experimental results
can be drawn, as shown in Figure 4(b).
The memory cost of these algorithms are listed in Table 3, from which we see more
prominent memory saving of the proposed algorithms over the randQB b algorithm.
While compared with randQB, the proposed algorithms consume comparable memory.
Lastly, we test a dense matrix with size n = 8, 000, and vary the value of l. The
trends of the runtime are plotted in Figure 5. It shows that the randQB EI and
randQB FP algorithms without the power scheme are about 1.9X and 2.5X faster
than randQB b, respectively. If the power scheme is imposed, the speedup ratios to
randQB b decrease, but randQB FP is still more than 2X faster than randQB b.
14

Table 3
The memory usage of the randomized QB-factorization algorithms for sparse matrices (l = 200)

n

randQB randQB EI randQB FP randQB b [1]

20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0

162
232
293
338
405

MB
MB
MB
MB
MB

174
239
303
359
426

MB
MB
MB
MB
MB

223
312
402
488
582

randQB
randQB_EI
randQB_FP
randQB_b [1]
Time(s)

Time(s)

16,000
24,000
32,000
40,000
48,000

200

400

600

800

1000

l

20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0

MB
MB
MB
MB
MB

6,153 MB
13,572 MB
23,917 MB
N.A.
N.A.

randQB(P=1)
randQB_EI(P=1)
randQB_FP(P=1)
randQB_b(P=1) [1]

200

400

600

800

1000

l

(a) Without the power scheme

(b) With the power scheme

Fig. 5. Runtime of the randomized QB-factorization algorithms for a dense matrix (n = 8, 000)
with varying value of l.

5.2. Comparison of accuracy. Three kinds of matrices are tested standing for
different distribution patterns of singular values:
• Matrix 1 (slow decay): A = U ΣV , where U and V are randomly drawn
matrices with orthonormal columns, and the diagonal matrix Σ has diagonal
elements σj = 1/j 2 .
• Matrix 2 (fast decay): A is formed just like Matrix 1, but the diagonal
elements of Σ is given by σj = e−j/7 . It reflects a fast decay of singular
values.
• Matrix 3 (S-shape decay): A is built in the same manner as Matrix 1
and Matrix 2, but the diagonal elements of Σ are given by σj = 0.0001 +
(1 + ej−30 )−1 . It makes the singular values first hover around 1, then decay
rapidly, and finally level out at about 0.0001.
For each kind, we generate a 2, 000 × 2, 000 matrix, for which we compare the errors
of the proposed techniques and the blocked randQB scheme [1] for varying l values.
The results are shown in Figure 6, where we see that the proposed techniques have
just the same accuracy as the blocked randQB algorithm. If we use the power scheme,
even with a power parameter as small as P = 1, the errors of the randQB EI and
randQB FP algorithms are remarkably reduced. And, the power schemes with P =
1 and P = 2 produce indistinguishable results for the tested matrices. Both are
extremely close to the optimal results from SVD.
15

0

10

SVD
randQB_EI(P=0)
randQB_EI(P=1)
randQB_EI(P=2)
randQB_FP(P=0)
randQB_FP(P=1)
randQB_FP(P=2)
randQB_b(P=0)
randQB_b(P=1)
randQB_b(P=2)

−1

||A−QB||/||A||

10

−2

10

−3

10

−4

10

0

20

40

60

80

100
l

120

140

160

180

200

(a) Errors of Matrix 1 whose singular values decay slowly
0

10

−2

||A−QB||/||A||

10

SVD
randQB_EI(P=0)
randQB_EI(P=1)
randQB_FP(P=0)
randQB_FP(P=1)
randQB_b(P=0)
randQB_b(P=1)

−4

10

−6

10

0

10

20

30

40

50
l

60

70

80

90

100

(b) Errors of Matrix 2 whose singular values decay rapidly
0

10

SVD
randQB_EI(P=0)
randQB_EI(P=1)
randQB_FP(P=0)
randQB_FP(P=1)
randQB_b(P=0)
randQB_b(P=1)

−1

||A−QB||/||A||

10

−2

10

−3

10

−4

10

0

10

20

30

40

50
l

60

70

80

90

100

(c) Errors of Matrix 3 with S-shape decay of singular values
Fig. 6. Errors of different algorithms for approximating the test matrices (n = 2, 000, b = 10).

5.3. Performance of the single-pass algorithm. Without the power scheme,
the randQB FP algorithm is a single-pass algorithm (see Algorithm 3). This is because
G = AΩ and H = AT G can be executed through one pass over matrix A, providing
that A is in the row-major format. Another single-pass algorithm was proposed in
[2], as a remedy to the randQB algorithm. It is shown in Figure 7, whose step (3)
T

produces matrices Q and Q̃ such that A ≈ QQT AQ̃Q̃ . Then, a small matrix
16

T

T

B̂ = QT AQ̃ is approximately solved in step (4), because Q̃ Ỹ ≈ Q̃ AT QQT Ω̃ =
T

B̂ QT Ω̃. This single-pass algorithm corresponds to the low-rank factorization in form
T
of A ≈ QB̂ Q̃ . Obviously, it includes more approximations and is not equivalent
to the randQB algorithm. In contrast, Algorithm 3 is mathematically equivalent to
randQB (see Proposition 2), and is supposed to be more accurate. With Matrix 1
and Matrix 2 from Sect. 5.2, we can compare the accuracy of the both algorithms,
whose results are plotted in Figure 8.
A single-pass algorithm for low-rank approximation
Input: A, l.
T

Output: Q, Q̃, B̂ such that A ≈ QB̂ Q̃ .
(1) Ω = randn(n, l); Ω̃ = randn(m, l);
(2) Y = AΩ; Ỹ = AT Ω̃;
(3) Q = orth(Y ); Q̃ = orth(Ỹ );
T

T

(4) B̂ = (Ω̃ Q)−1 Ỹ Q̃ .
Fig. 7. An existing single-pass algorithm for low-rank matrix factorization [2].
0

2

10

10

0

10

−1

||A−QB||/||A||

||A−QB||/||A||

10

−2

10

−2

10

−4

10

−3

10

SVD
single−pass[2]
randQB_FP

−6

10

−4

10

−8

0

50

100
l

150

10

200

0

20

40

60

80

100

l

Fig. 8. The approximation errors of the two single-pass algorithms and truncated SVD for
Matrix 1 (left) and Matrix 2 (right).

From Figure 8 we see that the approximation error of the single-pass algorithm in
[2] is often one order of magnitude larger than that of our randQB FP based algorithm.
Actually, it does not even decrease as the rank of the approximation matrix increases.
We also calculate the top 50 singular values, and the over-sampling with s = 10 is
applied to the both algorithms. The results are shown in Figure 9, along with those
obtained from the randQB algorithm, where the results of randQB FP and randQB
are indistinguishable. For the matrix with slow decay of singular value the result
from randQB FP shows moderate accuracy on the top singular values, whose error is
usually orders of magnitude smaller than that of the single-pass algorithm in [2].
5.4. Results of solving the fixed-precision problems. In this subsection
we test the proposed algorithms with some fixed-precision problems. The optimal
solution is the factorization with the smallest sizes of Q and B, which corresponds to
less amount of subsequent computation. It can be achieved by first calculating SVD
Pmin(m,n)
of the input matrix A, and then checking ( j=k+1 σj2 )1/2 , where σj is A’s j-th
17

0

0

10

10
single−pass [2]
randQB_FP
randQB
SVD

−1

−1

10

j

−2

10

σ

σ

j

10

−3

single−pass [2]
randQB_FP
randQB
SVD

−3

10

10

−4

10

−2

10

−4

10

20

30

40

10

50

10

20

30

40

50

Fig. 9. The computational results of top singular values for Matrix 1 (left) and Matrix 2 (right).

singular value, to determine the smallest k satisfying the accuracy criterion. Here, we
always consider the accuracy criterion with a relative tolerance: kA−QBkF < kAkF .
The proposed algorithms are compared with the SVD based method and the
adaptive randomized range finder (Algorithm 4.2 in [2]). The row-by-row calculation
scheme mentioned in Section 3 is implemented into our algorithms. The experiments
are carried out with Matlab on the aforementioned Linux server. The built-in commands like “svd”, “qr”, etc. are employed, which naturally take advantage of parallel
computing.
For the randQB FP algorithm, we empirically set ˜l = 50b. With a suitable value of
b, this produces a large enough ˜l for attaining the accuracy criteria in the experiments.
A more sophisticated approach for setting ˜l and b can be investigated in the future.
We first construct the three kinds of matrices in Sect. 5.2, each of 8, 000 × 8, 000
size, and test them with the four methods. Their results are shown in Table 4. For
randQB EI and randQB FP, the power scheme with P = 1 is used. The block size is
set to b = 10 in all tests, except the last one for which b = 40. In Table 4, “” stands
for the threshold for relative error, and “error” means the relative Frobenius-norm
error of the produced QB factorization. From the table we see that the results of
randQB EI and randQB FP algorithms all satisfy the set accuracy demands. And,
the corresponding ranks (i.e., the number of columns in Q) are very close to the
the optimal values from the SVD based approach. As for the runtime, the proposed
algorithms are usually several tens times faster than SVD. Notice that our Matlab
programs are less optimized than the built-in svd command. So, more significant
speedup could be expected for the implementation in C. Although the adaptive range
Table 4
The results on solving the fixed-precision problems based on three test matrices

Matrix

randQB EI



type

rank time(s)

error

randQB FP truncated SVD RangeFinder [2]
rank time(s) rank

time(s) rank

time(s)

Matrix 1

1e-2
1e-4

15
327

1.19
8.29

9.3e-3
9.98e-5

15
328

3.13
3.71

15
313

123

115
2,084

1.1
29.1

Matrix 2

1e-4
1e-5

66
82

2.16
2.68

8.37e-5
8.86e-6

66
82

2.73
3.17

65
81

115

101
113

1.1
1.1

1.56
18.7

4.1e-3
33
1.499e-3 1,587

3.62
15.8

32
1,587

126

3,618
7,916

87.8
379

Matrix 3

1e-2
33
1.5e-3 1,588

18

finder is built on a theory with spectral norm of matrix, in our experiments it always produces a QB factorization satisfying the accuracy demand in Frobenius norm.
However, its results (factor matrices) are much larger than necessary.
We then test the algorithms with two real data. One is from a scenic image [19],
and the other is from an information retrieval application “AMiner” [20]. The colored
image is represented by a 9, 504 × 4, 752 matrix. The other is an 8, 130 × 100, 000
keyword-person matrix produced with the term frequency and inverse document frequency (TF-IDF) model [21]. This sparse matrix has about 0.2% nonzero elements.
The computational results are listed in Table 5, with different power parameters and
block sizes. They again validate that the proposed algorithms can automatically
satisfy the accuracy criterion. And, with P = 2 the result of rank is substantially
reduced, approaching the optimal value. For the same power scheme, setting larger
block size b we can reduce the runtime of randQB EI. In contrast, the runtime of
randQB FP increases with the block size, as we have set ˜l = 50b. Notice that with the
relative error tolerance  = 0.1, the image is largely compressed (∼ 7X size reduction),
with little loss of quality (see Figure 10). And, the singular value of “AMiner” matrix
decays very slowly, but even with large approximation error its low-rank approximation could bring improved performance of information retrieval (c.f. [21], Sect. 11.3).
Table 5
The results on solving the fixed-precision problems based on two real data

Matrix  parameters

randQB EI
rank time(s) error

image 0.1

AMiner 0.5

randQB FP truncated SVD RangeFinder [2]
rank time(s) rank

P =1, b=10 468
P =1, b=20 468

8.1 0.0999 471
4.23 0.0999 472

3.25
4.44

P =2, b=10 441
P =2, b=20 441

9.98 0.0999 443
5.76 0.0999 443

3.47
7.26

426

P =1, b=50 2440 108 0.4999 2,449 143
2,115
P =2, b=50 2,229 134 0.4999 2,242 205

(a) original image

time(s) rank

time(s)

44.2

2,913

79.0

1,049

8,018

399

(b) compressed image (10% error)

Fig. 10. The original image and the compressed image obtained with the proposed algorithm.

For the second data, which is a large sparse matrix, the Krylov subspace iterative
method “svds” [18, 21] is also tested. However, it costs 2,281 seconds for computing
the first 1,000 singular values/vectors. It is much slower than executing “svd” to the
19

matrix’s dense version. Besides, “svd” requests more than 20 GB memory, while the
proposed randomized algorithms only costs 3 GB memory or so for this case.
6. Conclusions. Efficient techniques are proposed for the fixed-precision lowrank approximation of large matrices. Our contributions are as follows.
• A simple and accurate error indicator in Frobenius norm is proposed, which
enables efficient rank determination and can be used in the blocked randQB
algorithm [1] and other incremental QB-form factorization algorithms (like
that in [14]). We have proved its accuracy and validity for the problems with
relative accuracy tolerance larger than 2.1 × 10−7 . Numerical experiments
on large dense and sparse matrices have shown that the proposed rank determination scheme brings several to several tens times speedup and memory
saving to the blocked randQB algorithm, without loss of accuracy.
• Base on the blocked randQB algorithm, we propose a pass-efficient algorithm
called randQB FP. It is mathematically equivalent to the blocked randQB algorithm, but reduces the passes over matrix A to the fewest. The randQB FP
algorithm also suits to the fixed-precision problem, and can derive a singlepass algorithm under certain condition. Numerical results have validated the
efficiency and accuracy of the randQB FP algorithm, and shown that the
derived single-pass algorithm is much more accurate than an existing counterpart.
• Real data are tested to demonstrate the effectiveness of the proposed algorithms for the fixed-precision problem. Compared with the adaptive range
finder approach [2], the proposed algorithms run faster and produce much
smaller factor matrices while attaining the accuracy criterion.
Future work includes extending and applying the proposed algorithms to more
practical data mining and machine learning scenarios.
REFERENCES
[1] P.-G. Martinsson and S. Voronin, A randomized blocked algorithm for efficiently computing
rank-revealing factorizations of matrices, SIAM J. Sci. Comput., 38(2016), no. 5, pp. S485
- S507.
[2] N. Halko, P.-G. Martinsson and J. A. Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Review, 53
(2011), no. 2, pp. 217–288.
[3] P. Drineas and M. W. Mahoney, RandNLA: Randomized numerical linear algebra, Communications of the ACM, 59(2016), no. 6, pp. 80–90.
[4] S. Voronin and P.-G. Martinsson, RSVDPACK: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and GPU
architectures, arXiv preprint, arXiv:1502.05366v3 [math.NA], 2016.
[5] S. Eriksson-Bique, M. Solbrig, M. Stefanelli, S. Warkentin, R. Abbey, and I. C. F.
Ipsen, Importance sampling for a Monte Carlo matrix multiplication algorithm, with application to information retrieval, SIAM J. Sci. Comput., 33 (2011), no. 4, pp. 1689–1706.
[6] P. Drineas, R. Kannan, and M. W. Mahoney, Fast Monte Carlo algorithms for matrices
II: Computing a low-rank approximation to a matrix,, SIAM J. Sci. Comput., 36 (2006),
pp. 158–183.
[7] V. Rokhlin, A. Szlam, and M. Tygert, A randomized algorithm for principal component
analysis,, SIAM J. Matrix Anal. Appl., 31 (2009), no. 3, pp. 1100–1124.
[8] H. Ji and Y. Li, GPU accelerated randomized singular value decomposition and its application
in image compression, In Proceedings of Modeling, Simulation, and Visualization Capstone
Conference, Sulfolk, VA, 2014, 7 pages.
[9] T. Mary, I. Yamazaki, J. Kurzak, P. Luszczek, S. Tomov, and J. Dongarra, Performance
of random sampling for computing low-rank approximations of a dense matrix on GPUs,
In Proc. SC’2015, Austin, TX, Nov. 2015, 11 pages.
20

[10] M. W. Mahoney, Randomized algorithms for matrices and data, Foundations and Trends in
Machine Learning 3 (2011), no. 2, pp. 123–224.
[11] G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd ed., Johns Hopkins University
Press, Baltimore, MD, 1996.
[12] M. Gu and S. C. Eisenstat, Efficient algorithms for computing a strong rank-revealing QR
factorization, SIAM J. Sci. Comput., 17 (1996), no. 4, pp. 848–869.
[13] P.G. Martinsson, G. Quintana-Orti, N. Heavner, and R. van de Geijn. Householder QR
factorization with randomization for column pivoting (HQRRP), SIAM J. Sci. Comput.,
39 (2017), no. 2, pp. C96–C115.
[14] J. A. Duersch, and M. Gu, Randomized QR with column pivoting, SIAM J. Sci. Comput., 39
(2017), no. 4, pp. C263–C291.
[15] C. Eckart and G. Young , The approximation of one matrix by another of lower rank,
Psychometrika, 1 (1936), pp. 211–218.
[16] S. Voronin and P.G. Martinsson, RandQR, http://amath.colorado.edu/faculty/martinss/
main codes.html, 2015.
[17] Intel Parallel Studio XE Cluster Edition for Linux, https://software.intel.com/en-us/
intel-parallel-studio-xe, 2017.
[18] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM Press, Philadelphia, MD, 2003.
[19] http://flickr.com/photos/38459790@N03/4398318416.
[20] http://www.aminer.org.
[21] L. Elden, Matrix Methods in Data Mining and Pattern Recognition, SIAM Press, Philadelphia,
MD, 2007.
[22] Richard Lehoucq, Daniel Sorensen, and Chao Yang., ARPACK Users Guide: Solution of
Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM Press,
Philadelphia, MD, 1998.
[23] Joel A. Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher, Practical sketching algorithms for low-rank matrix approximation, SIAM J. Matrix Anal. Appl., 38 (2017),
no. 4, pp. 1454–1485.
[24] H. Li, G. C. Linderman, A. Szlam, K. P. Stanton, Y. Kluger, and M. Tygert, Algorithm
971: An implementation of a randomized algorithm for principal component analysis,
ACM Trans Math Softw., 43 (2017), no. 3, art. no. 28.

21

