IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

1Fast Randomized Singular Value Thresholding for
Low-rank Optimization
Tae-Hyun Oh Student Member, IEEE, Yasuyuki Matsushita, Yu-Wing Tai Senior Members, IEEE,and In So Kweon Member, IEEE
Abstract—Rank minimization can be converted into tractable surrogate problems, such as Nuclear Norm Minimization (NNM) and Weighted
NNM (WNNM). The problems related to NNM, or WNNM, can be solved iteratively by applying a closed-form proximal operator, called
Singular Value Thresholding (SVT), or Weighted SVT, but they suffer from high computational cost of Singular Value Decomposition (SVD) ateach iteration. We propose a fast and accurate approximation method for SVT, that we call fast randomized SVT (FRSVT), with which we
avoid direct computation of SVD. The key idea is to extract an approximate basis for the range of the matrix from its compressed matrix.
Given the basis, we compute partial singular values of the original matrix from the small factored matrix. In addition, by developping a range
propagation method, our method further speeds up the extraction of approximate basis at each iteration. Our theoretical analysis shows therelationship between the approximation bound of SVD and its effect to NNM via SVT. Along with the analysis, our empirical results
quantitatively and qualitatively show that our approximation rarely harms the convergence of the host algorithms. We assess the efficiency
and accuracy of the proposed method on various computer vision problems, e.g., subspace clustering, weather artifact removal, andsimultaneous multi-image alignment and rectification.
Index Terms—Singular value thresholding, rank minimization, nuclear norm minimization, robust principal component analysis, low-rankapproximation

F

1

I NTRODUCTIONM

I nimizing the rank of a matrix can be used as a versatileregularizer to derive a low-rank solution and is neededin many mathematical models in computer vision and machinelearning [5], [10], [14], [21], [23], [29], [31], [32], [36], [43], [45].
As rank minimization, min rank(X) s.t. X ∈ C , is generally
an NP-hard problem, it is typically relaxed using the nuclear
norm (i.e., k · k∗ , the sum of the singular values), which is atight convex surrogate for the rank function, i.e., rank(·).
Nuclear norm minimization (NNM) is expressed as

X∗ = argminX f (X) + τ kXk∗ ,(1)where X ∈ Rm×n , and τ > 0 is a regularization parameter. The function f (X) can be arbitrarily defined dependingon the objectives, e.g., f (X) = kO − Xk1 in robust principal component analysis (RPCA) [4], f (X) = 21 kAX − Bk2F
in multivariate regression and multi-class learning [26],f (X) = 12 kπΨ (O) − πΨ (X)k2F in matrix completion [4],f (X) = kOX − Ok2F in subspace clustering [21], where O ismeasured data, k · k1 and k · kF are l1 and Frobenius norms
respectively, and πΨ (·) is an orthogonal projection operator
setting [πΨ (X)]i,j = [X]i,j for (i, j) ∈ Ψ and 0 otherwise.
For the purpose of better approximating the rank function
rank(·), weighted nuclear norm minimization (WNNM) [10],[14], [30], which could be non-convex depending on theweights that are used, can be alternatively adopted instead of
the standard nuclear norm.
Except for the case that the loss function f (X) is a proximity term, most previous works use first-order optimization•••

T.-H. Oh (corresponding author), and I.S. Kweon are with the Department ofElectrical Engineering, KAIST, Daejeon, Republic of Korea.
E-mail: thoh.kaist.ac.kr@gmail.com and iskweon77@kaist.ac.kr
Y. Matsushita is with Osaka University, Japan.
E-mail: yasumat@ist.osaka-u.ac.jp
Y.-W. Tai is with Tencent.
E-mail: yuwing@gmail.com

approaches, e.g., dual method [8], accelerated proximal gradient [15], augmented Lagrange multiplier (ALM) [18], alternating direction method (ADM) [18], [19], and iteratively
reweighted least-squares [24]. In the intermediate step, regardless of regularization, all these approaches have an iterative
step to solve a simple NNM, or WNNM, subproblem definedas the following nuclear norm and proximity terms:
Problem (Nuclear norm minimization). For τ ≥ 0 and
A ∈ Rm×n ,1(2)
X∗ = argminX τ kXk∗ + kX − Ak2F ,2
where optimal X∗ can be obtained by the singular value
thresholding operator defined as following.
Definition 1 (Singular value thresholding (SVT)1 [2]). The
problem (2) has a closed-form solution given by the singular value
thresholding operator Sτ (·) as>X∗ = Sτ (A) = UA Sτ (ΣA )VA,(3)

where Sτ (x) = sgn(x) · max(|x| − τ, 0) is a soft shrinkage opera>
tor [11], and UA ΣA VA
is the SVD of A.
The major computational bottleneck of NNM and WNNM
problems is the necessity of solving Eq. (2) multiple times,
where SVD computation occupies the largest computation cost,i.e., O(mn min(m, n)) for each SVD [9].
This paper proposes a fast SVT technique to accelerate general NNM and WNNM computation. Our method is motivatedby the previous study of a randomized SVD proposed by
Halko et al. [12], and we extend the original general method inseveral respects for better solving NNM and WNNM problems
that we focus on in this paper. As a result, we propose analgorithm that we call fast randomized SVT (FRSVT). We present
1

A similar result for WNNM can be found in [10], called WSVT.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEARπΨ (·)Sτ (·)Sτ (·)σi (·)
Ω
Q
O
L, L
S, E
k
p
l
r
s
η
b
v

Orthogonal projection operator with a map (index set) Ψ
SVT operator [2] with the parameter τSoft shrinkage operator [11] with the parameter τi-th singular value of a matrix
A standard Gaussian random matrix
Orthonormal column matrix
Observation matrix (Input)
Low-rank optimization related matrix or tensor
Sparse optimization related matrix or tensor
Parameter for the target rank
Parameter for over-sampling rate
Sampling (predicted) rate (l = k + p)
True rank of the target matrix A
Real sampling rate (s = min(l, r))
Number of the power iteration
Parameter for the maximum rank boundParameter vector for the polynomial error bound (v = {k, p} without the power iteration, or v = {k, p, η} with it)

TABLE 1: Summary of notations

the connection between FRSVT and low-rank approximation
with both theoretical and empirical analyses, and show theeffectiveness of the proposed method via simulations and a
few computer vision applications. Table 1 summarizes thenotations used in the rest of the paper.
A preliminary version of this paper has appeared in [34].
We extend [34] by analyzing the behavior of the proposed
method in both theoretical and empirical aspects as follows.
Detailed derivations of error bounds on both Frobenius andspectral norms, inexpensive and provable accuracy check, andcomparison with a linear-time approximate SVD [6] are newly
supplemented in Sec. 3. Based on the analyses, we provide
theoretically sound implementation tips and other options
that can potentially further improve the performance of the
proposed algorithm. Additional simulation and experimental
results are included in Sec. 4. Specifically, this paper makes thefollowing contributions:• We develop a successive truncated low-rank decomposition
method that can be generally applied to NNM and WNNM
problems. Our method achieves high efficiency with rarely
degrading the accuracy of the host algorithms.
• By exploiting the proximity of the range space over iterations,
our method further accelerates NNM computation by propagating the estimated basis at the previous iteration to the next
one for acceleration. We call this technique as range propagation(RP).
• We provide a theoretical analysis between the low-rankapproximation and our FRSVT method. We show how the
low-rank approximation affects to the SVT operation as wellas the interaction with the power iteration and over-sampling.
In addition, we show the empirical stability and behavior of
our method with respect to varying parameters.
• We apply FRSVT to various computer vision applications
and show the performance gain in comparison with previous
methods.
1.1

2high dynamic range imaging [32], batch image alignment [36],photometric stereo [43], image rectification [45], and nuclear
norm regularized learning [26].
While effective, due to the high computational complexityof solving NNM, specifically, SVD used in the SVT operator, afast SVT method has always been wanted in both small- andlarge-scale problems. Liu et al. [22] efficiently solve NNM bycasting the original convex RPCA problem into a bilinear factorization form. Although this bilinear factorization introduces
non-convexity in the objective function like other factorization
methods [7], [46], the work achieved significant speed-up.
Unlike their method, our method retains the original objectivefunction and approximately optimizes it, and is applicable togeneral NNM problems once the problems can be led to theform of Eq. (2).
With retaining the advantage of convexity, Liu et al. [23]exactly solve RPCA on a small sub-sampled matrix and propagate the seed solution to other parts via `1 filtering. Since both
Liu et al. [22] and Liu et al. [23] focus only on RPCA problem but
not general NNM problems, a fast and general SVT method isstill needed as a tool for NNM to be applied to large-scale problems. Cai et al. [3] avoid explicit SVD computation using thedual of SVT. Since their method uses Newton iterations withan inverse matrix, the input matrix needs to be preprocessedby the complete orthogonal decomposition [9] (COD) to ensurea non-singular square matrix. The standard COD consistsof twice of QR decompositions with column/row pivoting,which requires O(mn min(m, n)); therefore, the reduction ofcomputation complexity is still limited.
In other thread of works, it has been shown that theexact SVD computation is unnecessary in the inner loops of
NNM [18], [26], [29]. Mu et al. [29] propose a compressed optimization by random projection. Ma et al. [26] solve NNM related problems with a linear-time approximate SVD [6]. However, these methods become occasionally unstable, because anintermediate matrix is directly approximated by sampling or
projection, where the original information is impaired, andthe randomness leads to unstable and incorrect results. Unlike
these methods, our method only approximates the subspacebases to guarantee that most spectrum information is retained.
Scope of this work This work mainly focuses on the behavior of the fundamental core module, i.e., SVT operator,
combined with low-rank approximation. All the provided
theoretical analyses in this work are mainly focused on the
FRSVT operator per se, rather than the overall optimization
procedure of low-rank optimization problems. In rank minimization literature, there are approaches [28], [33], [41] that are
tightly entangled with the optimization procedure to promote
low-rank solutions, and they are beyond SVT like modular
operation. This paper does not cover this type of approaches,
yet focuses on the SVT operation that is used broadly [5], [10],[14], [21], [26], [32], [36], [43], [45].

Related Works

Candès et al. [4] showed that, under some mild condition,
the solution of NNM is equivalent to the solution of rank
minimization in conjunction with a sparse outlier model [4]in the RPCA problem setting. Inspired by the success of thisconvex surrogate for rank minimization, low-rank structures
have been exploited in various computer vision applications,such as rain removal [5], de-noising [10], inpainting [14], motion segmentation by subspace clustering [21], structure frommotion [23], background subtraction [29], tag transduction [29],2

FAST R ANDOMIZED S INGULAR VALUE T HRESH -OLDING (FRSVT)The heart of our method shares the idea of Liu et al. [23] andMa et al. [26], in that the solution of Eq. (2) can be found by
applying SVT to a small matrix instead of the original large
matrix as illustrated in Fig. 1. Instead of sampling columns
or rows of a matrix as in [23], [26], our method extracts
a small core matrix by finding orthonormal bases with theIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

Algorithm 1 Fast Randomized Singular Value Threshold-𝕊𝜏 [ ]𝕊𝜏 [ing (FRSVT) algorithm]=

1: Input : A ∈ Rm×n , τ > 0, l = k + p > 0 and q ≥ 0. For range propagation,
the orthonormal column matrix Q̃ of the previous iteration.

Fig. 1: Basic idea of our method. Instead of applying SVT to the large
original matrix, if we can obtain the same result by applying SVT to small
matrix with few additional efforts, the complexity could be significantly drop.
As shown in Proposition 1, once a large matrix is decomposed into a corematrix and orthonormal matrices (the left and right thin matrices in theillustration), by applying SVT to the small core matrix we can obtain thesame result with the direct SVT computation on the large matrix.

unitary invariant property. Specifically, since the NNM definedin Eq. (2) consists of unitary invariant norms, the following
equality holds:
Proposition 1. Let A = QB ∈ Rm×n , where Q ∈ Rm×n has
orthonormal columns. Then,Sτ (A) = Q Sτ (B),2:
3: if not Range propagation then4:
Sample Gaussian random matrix Ω ∈ Rn×l
5:
Y = AΩ
6:
Q = QR CP(Y)
7: else8:
Sample Gaussian random matrix Ω ∈ Rn×p
9:
Y = AΩ
10:
QY = PartialOrthogonalization(Q̃, Y)
11:
Q = [Q̃, QY ]12: end if
13: repeat
14:
Q = QR(AA> Q)
15: until η times
16: [H, C] = QR(A> Q)
17: [W, P] = PolarDecomposition(C)18: [V, D] = EigenDecomposition(P)
19: Sτ (A) = (QV) Sτ (D) (HWV)>
20:
21: Output : Sτ (A), Q̃ = QV

(4)

where Sτ (·) is the SVT operator.
Proof. When m = n and Q is an orthonormal matrix, the equality obviously holds by the unitary invariant property of norms.
For m > n, we prove the equality for the orthonormal column>
matrix. Consider an arbitrary matrix Z = UZ ΣZ VZ∈ Rn×n ,>where UZ ΣZ VZ is Z’s SVD, then
1
2
Q · Sτ (B) = Q · argmin τ kZk∗ + kZ − BkF
2
Z


1
= argmin τ kQ> Xk∗ + kQ(Z − B)k2F
2
X(by letting QZ = X, and as Q> Q = I,
and the unitary invariant norm property)1= argmin τ kXk∗ + kX − QBk2F
2
X

SVD computation, we can further reduce the computation as
described in Sec. 2.2. Also, we describe a target rank prediction
technique for further improving speed in Sec. 2.3.
2.1

= Sτ (QB) = Sτ (A),which concludes the proof.
In general, SVT requires SVD computation, and its complexity2 is O(mn2 ). Based on Proposition 1, we may avoid
expensive computation by instead computing SVT on a smaller
matrix B ∈ Rn×n when Q ∈ Rm×n is available and m ≥ n.
Given Q ∈ Rm×k that best approximates A by a rank-k matrix,
the complexity of SVD for B ∈ Rk×n becomes O(nk 2 ). Therefore, when k  n, the computation speed can be significantly
improved.
Our SVT computation iterates the following two steps: 1)
Estimating an orthonormal column matrix Q, and 2) Computing SVD of B for SVT. For SVT computation, a partial (ortruncated) SVD is frequently used to reduce the complexity inmany prior arts. Our method similarly finds a rank-k approximation (k < n) of the original matrix A as A ≈ Âk = QB. It
saves the computation of the first step as well as the second
step, because the size of matrices Q and B are reduced tom × k and k × n, respectively. By exploiting the observation that the major orthonormal k bases evolve slowly overiterations, our method efficiently initializes matrix Q at each
iteration by bypassing expensive random range estimation(Range propagation in Sec. 2.1). In addition, by avoiding direct
Without loss of generality, we assume m ≥ n in this paper.

Finding Approximate Range

Inspired by Halko et al. [12], we first estimate the orthonormal bases Q = [q1 , · · · , ql ] (where l ≥ k ) such thatspan(Q) ⊆ Range(A)3 from a matrix compressed by random
projection. Intuition of their randomized range finding algorithm is as follows. By multiplying a random vector ωj , arandom linear combination yj of the column vectors of A
is generated, which encodes the partial range of A. Suppose
A = Ak + E, where Ak is the best rank-k projection of A,
whose range is the target to be captured, and E represents
small perturbation, then a sample vector yj can be obtained by(since kQ> Xk∗ = kZk∗ = kQZk∗ = kXk∗ )23

yj = Ak ωj + Eωj .

(5)

Even though unwanted E may be included in {yj } because
the action of Ak (i.e., magnitudes of spectrum) is larger than E,
the range of Ak is dominant to be captured in {yj }. However,
if only k vectors {yj } are sampled for rank-k approximation,{yj } could not span the entire Range(Ak ). Instead, by increasing the sampling rate, most of Range(Ak ) can be captured;
therefore, we oversample l sample vectors, so that Range(Ak )can be captured as much as possible.
Given the sample matrix Y = [y1 , · · · , yl ], where l = k + p,
Q can be obtained by orthonormalizing Y. Whenr = rank(A) < l, r bases are enough to span the entire
Range(A). Indeed, when rank(A) < l, the range finding algorithm is fairly accurate and close to the exact method as
we will see the theoretical analysis in Sec. 3. Thus, we canreduce the dimension of Q for almost free by estimating therank and dominant bases by QR decomposition with Column
Pivoting (QR-CP) to Y . While Halko et al. also proposedan algorithm to adaptively and approximately determine thenumber of bases by different randomization for sampling,orthogonalization, and re-orthogonalization, our method isbased on an exact method with the same complexity using
QR-CP, and we adaptively predict the sampling rate in a
simpler manner (we will see in Sec. 2.3 or Sec. 3.2). Fortunately,3

We only consider the column range space of a matrix for simplicity of
explanation, but the row range also can be used.

Normalized magnitude

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR
1ImageGaussian
Video

0.5

0
0
100
200
Singular value index

Angular Error (Degree)

Fig. 2: Illustration of singular value decaying. [Left] Gaussian random, video
and image samples. [Right] Decaying graphs of singular values. The Red,
Green, Blue lines represent the graphs of a Gaussian random matrix, video,
and image, respectively. The spectrum of visual data decays significantly
faster.
4

Computing the Singular Values (Vectors)

The NNM problem is now reduced to SVT on a smaller
matrix B. In this section, we further reduce the computation
time of SVT on B. The SVT operator can be computed by
SVD and shrinkage on its singular values. For positive semidefinite matrices, SVD can be more efficiently computed by
Eigen decomposition (ED), which is typically faster than SVD
in our empirical tests with small matrices. To apply ED to a
general matrix, we form a positive semi-definite matrix by thefollowing decomposition:
Definition 2 (Polar decomposition [13]). Let X ∈ Cm×n ,m ≥ n. There exists a matrix W ∈ Cm×n and a unique Hermitian
positive semi-definite matrix P ∈ Cn×n such that

R=10, OR=0.1
R=50, OR=0.1
R=10, OR=0.2
R=50, OR=0.2

3

2.2

4

2
1
0
0

X = WP,
5

1015
No. iteration

20

W∗ W = I,

25

Fig. 3: Angular difference of subspaces between subsequent iterations. R:
ground-truth rank of input data, OR: oversampling rate, p/k.

in LAPACK, an efficient QR-CP routine using level-3 BLAS
(dgeqp3) is available. Moreover, our method does not require
the upper triangle matrix from QR, but only the orthonormal basis Q; therefore, we can avoid extracting the whole
triangular matrix but only compute rank and Q with dgeqp3
and orgqr routines, respectively. After obtaining Q ∈ Rm×s ,where s = min(l, r), we can compute a small matrix B by
projection on Q as B = Q> A ∈ Rs×n .
Power Iteration Among the overall process in our algorithm,
since the only approximation step is the estimation of orthonormal column matrix Q, the accuracy of our algorithm depends
only on this step. In Eq. (5), if the magnitude of the action of
Ak is not dominant against E, the directions of sample vectors
are biased and may be affected by portions of Range(Ak )⊥ .
This introduces accuracy degradation to the rest of the process.
To resolve this issue, Halko et al. [12] used a power iteration
scheme, which makes the spectrum difference between Ak andE larger by estimating Q on (AA> )η A. It improves the chanceof better capturing the range of Ak from Y = (AA> )η AΩ,
while the range space (left singular vectors) remain unchanged.
Halko et al. also showed that η = 2 or 4 power iterations are
sufficient for usual data of interest, and highly accurate range
finding can be achieved. As shown in Fig. 2, decay of singularvalues of visual data is much faster than the one of a Gaussian
random matrix. Our empirical tests also show that η = 2 issufficient enough, and it is used in all our experiments.
Range Propagation (RP) for Fast Range Finding As shownin Fig. 3, we observe that Range(A(i) ) at the i-th iteration issimilar to the one at the (i − 1)-th iteration in NNM related
problems. This motivates us to use the singular vectors at the(i − 1)-th step as an initial approximation of range bases Q(i)at the i-th step. A similar observation is exploited in [20].
To capture the change of the range space, we additionallysample p number of sample vectors {y}. We append thesampled {y} to the previous singular vector matrix Q̃(i−1) as
Q(i) = [Q̃(i−1) , y1 , · · · , yp ], and apply partial orthogonalization only for newly added {y} by a modified Gram-Schmidt
procedure [9]. The number of bases can subsequently be reduced by checking the rank with QR-CP in the first step of thepower iteration or by a method described in Sec. 3.2.

where I is the identity matrix. If rank(X) = n, then P is positive
definite and W is uniquely determined.
Note that the existence of polar decomposition is equivalent
to the existence of SVD.
We use a Newton based polar decomposition suggested by
Higham et al. [13], which has a quadratic convergence behavior.
In our experiment, it converges at a small number of iterations (typically, seven) with various different data, which isconsistent with the result of [3], [13]. Due to the requirement of
the inverse operator in Newton iterations, it is only applicable
to non-singular square matrices. Since B> ∈ Rn×s is a fullcolumn rank matrix, the non-singular square matrix can be
simply obtained from B> = HC by QR decomposition, where
we call C ∈ Rs×s a core matrix that is always non-singularand square. For this step, unlike the procedure in Sec. 2.1, nocolumn pivoting is required.
We sequentially apply the polar decomposition and ED on
the core matrix; C = WP = WVDV> , where D and V are
the eigenvalue and eigenvector matrices of P, respectively.
Since the matrices H, W, and V are orthonormal column
matrices, the diagonal matrix D is equivalent to the singularvalue matrix of B. Finally, Sτ (A) can be approximated bySτ (A) ≈ Sτ (Âs ) = (QV) Sτ (D) (HWV)> .

(6)

For the range propagation, the singular vector matrix Q̃ isstored as Q̃ = QV or HWV (according to the side of randommatrix multiplication). Overall algorithm is summarized in
Algorithm 1.
2.3

Adaptive Rank Prediction (AP) Heuristic

For SVT, only singular vectors corresponding to the singularvalues that are greater than a certain threshold are needed, andfull SVD is unnecessary. Since the rank of A(i) is unknown
before SVD, predicting its rank can avoid unnecessary computation. We observe that, in many NNM related problems, therank of A(i) tends to monotonically increasing or decreasingover iterations, and the rank is stabilized as the number of
iteration increases. As we shall see in the theorem of error
bound in Sec. 3, over-sampling is always useful to reduce theexpected error bound of FRSVT. Thus, optimistically predicting rank allows to achieve both computational efficiency andstability.
The efficiency of our method may be degraded by anexcessively high sampling rate. In such a case, we resort tothe truncated SVT by upper bounding the target rank. As

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

shown in natural image statistics of Fig. 2, the rank of A(i)is generally stabilized at low-rank in many computer visionapplication. Usually, the final accuracy is not harmed, as seenin the successes of the truncated SVD in the NNM related
problems [18], [19], [20], [22].
Based on these observations, we define over-sampling ratep as:
pi+1 =a,dρne,

if ri < li ,otherwise,

5Line No.

Operation
m×l

1
2

z}|{
z}|{
Q = QR CP( Y )Fast Random Projection

We can replace the Gaussian random matrix with other lessexpensive and fast transforms according to the JohnsonLindenstrauss lemma [27], such as a subsampled random
Fourier transform, subsampled random Hadamard transform,or sparse Gaussian random matrix (in Matlab, sprandn() function). The cost of the sparse Gaussian random matrix transform
depends on the sparsity of the matrix, but it may becomeless accurate if the matrix is too sparse. The two structured
subsampled random matrices can be multiplied with the costO(mn log l) by efficient computation methodology (e.g., FFT),
while the cost of dense Gaussian matrix multiplication isO(mnl). In our implementation, we simply use dense Gaussian
random matrices drawn from the standard normal distribution
in entry-wise.
2.5

Computational Complexity

The step-by-step computational complexities are summarizedin Table 2. It does not consider the range propagation procedure for simplicity, but the range propagation only reduces thecomputational complexity. For example, while the computationof the sample matrix Y by random projection requires O(mnl),
we can reduce it to O(mn) by the range propagation technique.
Since Lines 3 and 4 (i.e., power iteration) in Table 2 iterates ηtimes, the overall algorithm in the case of m < n is much more
efficient than the other way around. We ensure A to be m < n,if necessary by transposing it, at the beginning of the algorithm.
The orthonormalization operation of Line 4 in Table 2 (Line 14
in Alg. 1) is used for enhancing the numerical stability during
the power iteration, but in practice it may not be necessaryfor every step and could be skipped except at the last step.
However, in this work, we do not omit it but consistently adopt
the orthonormalization every step as described in Alg. 1 forsimplicity of evaluation.

O(mnl)
O(ms2 )

m×s
n×m

3
4

m×n
m×s
z}|{
z}|{ z}|{ z}|{>Y = A A
Q

m×s
m×s

5

O(mns)

m×s

z}|{
z}|{
Q = QR( Y )
n×s

2.4

m×n n×lz}|{z}|{ z}|{Y = A Ω
m×l

(7)

where m ≥ n is assumed, a is a constant (set to a = 2),
ρ ∈ (0, 1] is a constant parameter to rapidly follow the realrank ri (set to ρ = 0.05 or less), and d·e denotes the ceil
operation. The prediction rule for sampling rate is defined
as li+1 = min(ri + pi+1 , b), where b = dγne is the maximum
bound of sampling rate, γ ∈ (0, 1] is the proportion parameter.
The sampling rate li can be regarded as the predicted rank
at the i-th iteration, and ri is the number of singular valuesof A(i) that are larger than the threshold, i.e., the estimatedrank of Sτ (A(i) ). Initially, we set l0 = 0.1b. In the case of therange propagation, the number of columns in Q̃(i−1) becomesri , and pi = li − ri . When ri < li , Eq. (7) slightly oversamples, otherwise it optimistically predicts the rank of thenext iteration by a larger over-sampling rate. By virtue of lowrankness of visual data shown in Fig. 2, the optimistic ruleleads to an accurate estimate.

Complexity

O(ms2 )

n×m m×sz}|{ z}|{
z}|{
>
Q
B = A
n×s s×s

O(ms2 )n×s67z}|{ z}|{z}|{[ H , C ] = QR( B )Polar & Eigen Decomposition with

8

(Lines 17, 18 in Alg. 1)
Composition Sτ (A) (Lines 19 in Alg. 1)

C|{z}

O(ns2 )O(s3 )s×s

O(mns)TABLE 2: Computational complexity of Alg. 1. For simplicity, we do notinclude the range propagation case, but it can be easily deduced by theabove table because most of basic operations remain similar. Incorporating
the range propagation obviously requires cheaper computation than theabove algorithm.

3

A NALYSISThe proposed method provides an approximate solution of
SVT; thus, a natural question is how accurate our algorithm
is. This section analyzes this respect by first establishing therelationship between SVT and low-rank approximation. Then,
we show that, under what conditions, our method can be expected to be accurate, whereby a few design tips are describedfor potentially further improving the algorithm. Besides, weshow the bound comparison between our method and a representative work, Linear-Time SVD (LTSVD) [6].
3.1 Relationship between SVT and low-rank approximation
We can see the relationship between SVT and low-rank approximation by comparing the low-rank approximation gap
and SVT operation gap for two different matrices. We begin
with reviewing some useful relationships.
The SVT operator Sτ (·) is a proximal operator [1]; thus
there exists its corresponding dual operator Pτ (·) such that
X = Sτ (X) + Pτ (X). The dual relationship is valid by thefollowing 2-norm ball Euclidean projection operator [3].
Definition 3 (2-norm ball Euclidean projection). For τ ≥ 0,>
consider SVD of Y = UY ΣY VY. Then the 2-norm ball Euclidean
projection operator is defined by>Pτ (Y) = UY Pτ (DY )VY,

where Pτ (x) = min(x, τ ) is an element-wise operation.
Given the 2-norm ball Euclidean projection operator, wehave the following lemma with proof for completeness.
Lemma 1. The 2-norm ball Euclidean projection operator is theclosed-form optimal solution of the following problem:Pτ (Y) = min kY − Xk2F .
kXk2 ≤τ

Proof. Since the proof is straightforward, we instead provide
the idea of the proof. By von Neumann’s inequality, the optimal
solution of X should have the same singular vectors with Y .
Since the l2 -norm and the Frobenius norm are both unitaryinvariant, the problem is reduced to a diagonal matrix form.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

Thus, the projection operator has an explicit form as shown in
Definition 3.

6

where Ak is the unique and exact rank-k truncated matrix ofA, and PQ = QQ∗ is an orthogonal projector. Then,

In order to derive a tight error bound, we need a contraction
inequality, called pseudo-contraction4 [37]:===Proposition 2 (Pseudo-contraction [37]). Let A, B ∈ Rm×n andSτ (·) be a proximal operator. Then, the following pseudo-contraction
is satisfied.

=kSτ (A) − Sτ (B)k2F≤ kA − Bk2F − k(Sτ (A) − A) − (Sτ (B) − B)k2F=== kA − Bk2F − kPτ (B) − Pτ (A)k2F .
The pseudo-contraction holds for all the proximal operators. For deriving a spectral norm based error bound, we derive
the spectral pseudo-contraction.
Proposition 3 (Spectral Pseudo-contraction). Let A, B ∈Rm×n and Sτ (·) be a proximal operator. There exists a constant
1 < C ≤ 2 satisfyingkPτ (A) − Pτ (Âk )k2F
kPτ (A) − Pτ (PQ A)k2F>> 2
kPτ (UA ΣA VA) − Pτ (PQ UA ΣA VA
)kF
>
> 2
kUA Pτ (ΣA )VA − PQ UA Pτ (ΣA )VA kF(by the property of the unitary transform)kUA Pτ (ΣA ) − PQ UA Pτ (ΣA )k2F
(by the unitary invariant of the norm)k(I − PQ )UA Pτ (ΣA )k2F
2
kP⊥
Q UA Pτ (ΣA )kF .

P⊥
Q = I − PQ is the orthogonal projector onto the complementary subspace Range(PQ )⊥ , which is also close toRange(Ak )⊥ . Thus, the minimum is achieved only when
Range(PQ ) = Range(Ak ) or Range(PQ )⊥ = Range(Ak )⊥by the well known variational characterization. Therefore,X2
min kP⊥min(σj (A), τ )2 .
Ak UA Pτ (ΣA )kF =j>kkSτ (A) − Sτ (B)k2≤ CkA − Bk2 − k(Sτ (A) − A) − (Sτ (B) − B)k2 .

(8)

Proof. We have the following inequalities for Sτ (·) and Pτ (·)
from non-expansiveness of proximal operator [1]:kSτ (A) − Sτ (B)k2 ≤ kA − Bk2 ,kPτ (A) − Pτ (B)k2 ≤ kA − Bk2 .
By adding two inequalities, we havekSτ (A) − Sτ (B)k2 ≤ 2kA − Bk2 − kPτ (A) − Pτ (B)k2 ,and then this implies there exists C ≤ 2. For the lowest valueof C , using A = Sτ (A) + Pτ (A) and triangle inequality, wehavekA − Bk2⇒ kA − Bk2= kSτ (A) + Pτ (A) − (Sτ (B) + Pτ (B)) k2≤ kSτ (A) − Sτ (B)k2 + kPτ (A) − Pτ (B)k2− kPτ (A) − Pτ (B)k2 ≤ kSτ (A) − Sτ (B)k2 ,which implies C > 1. We can conclude there exists a constant
C such that 1 < C ≤ 2 satisfying Eq. (8).
Since both pseudo-contractions involve the error term between two projection operators, i.e., kPτ (A) − Pτ (B)k, thefollowing lemma is useful to see the bound.
Lemma 2. Let Âk be a rank-k approximation of A. Then,min kPτ (A) − Pτ (Âk )k2F=

X

min kPτ (A) − Pτ (Âk )k2=min(σk+1 (A), τ ).

j>kmin(σj (A), τ )2 ,Proof. Let Âk = QQ∗ A = PQ A, where Q ∈ Rm×k consists ofk dominant orthonormal bases, and Range(Q) ' Range(Ak ),4

Refer to Lemma 3.3 of Pierra et al. [37]The spectral norm case can be derived similarly.
Other than the randomization step in sampling Y , ourmethod is based on exact algorithms. Hence, the accuracy
is only affected by the randomized range estimation, which
computes a rank-k approximation of a matrix. Since there is
only one approximation step, the proposed method shares thesame theorems with Halko et al. [12] as:
Theorem 1 (Average Frobenius error bounds by randomization [12]). Let Âk be a rank-k approximation matrix A ∈ Rm×n .
For a target rank k ≥ 2 and an over-sampling parameter p ≥ 2,
where k + p ≤ min(m, n), draw a standard Gaussian matrixΩ ∈ Rn×(k+p) . Then, the theoretical minimum error and the upperbound satisfyXX
σj2 (A) ≤ EkA − Âk k2F ≤ poly(v) ·σj2 (A),j>kj>kwhere E denotes expectation with respect to the random matrix,poly(v) is a function for v = {k, p} without the power iterationor v = {k, p, η} with the power iteration.
Halko et al. [12] show that the upper bound of the error isclose to the theoretical minimum error with high probability
in conjunction with some improving techniques (e.g., oversampling and power iteration), and the bound is rather pessimistic.
In Theorem 1, in the case without power iteration, poly(·)
is defined as:poly(v = {k, p}) = (1 + k/(p − 1)) .

(9)

The polynomial bound with power iteration in the Frobenius
norm representation has no simple form [12]; therefore, instead
we observe the error bound with power iteration by referringto the spectral bound in the following Theorem 2 and Corollary 1.
Theorem 2 (Average spectral error bound by randomization
with the power scheme [12]). Let Âk be a rank-k approximation
matrix A ∈ Rm×n . For a target rank k ≥ 2 and an over-samplingIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

parameter p ≥ 2, where k + p ≤ min(m, n), draw a standard Gaussian matrix Ω ∈ Rn×(k+p) . Then, the theoretical minimum error
and the upper bound satisfy
σk+1 ≤ EkA − Âk k2 ≤

√
q


e k+p
2η+1
k
 1+
σ
+
k+1
p−1
p

P

1
! 1  2η+1
22(2η+1)

,
σj

(10)

j>k

where η denotes the number of power iterations.
Corollary 1 (Loose average spectral error bound [12]). Theorem
2 is bounded as
EkA − Âk k2 ≤ [Bound of Theorem 2] ≤ poly(v) · σk+1 , (11)1hi 2η+1q√√kwhere poly(v) = 1 + p−1
h
−
k
, v =
+ e k+p·p{k, p, η}, and h = min(m, n).
Corollary 1 is only useful for understanding the boundbehavior rather than identifying the tight bound, because it
is far broader than Theorem 2. However, it shows that, as thenumber of power iterations η increases, poly(v) approaches 1
exponentially fast, and when poly(v) ∼ 1, the spectral error is
likely to be bounded by the theoretical minimum σk+1 .
For our method, now we can derive the following Frobenius error bound for the approximate SVT (i.e., FRSVT) using
Theorem 1, Proposition 2 and Lemma 2.
Theorem 3 (Average Frobenius error bound of the approximate SVT). Let Sτ (·) be the SVT operator and Âk holds Theorem 1.
The average error satisfies the following inequality:X
σj2 (A) − G(A),EkSτ (A) − Sτ (Âk )k2F ≤ poly(v) ·j>kP2
where G(A) =
min(σ
(A),
τ
)
≥
0
.
j>k
j
Proof. By Proposition 2 and Lemma 2,

kSτ (A) − Sτ (Âk )k2F ≤ kA − Âk k2F − kPτ (A) − Pτ (Âk )k2F
X
≤ kA − Âk k2F −min(σj (A), τ )2 .
j>k

By taking the expectation of both sides,

70where 1 < C ≤ 2, and G (A) = min(σk+1 (A), τ ) ≥ 0.
Proof. We can simply derive Eq. (12) in a similar manner to
Theorem 3, but instead with Theorem 2 and Proposition 3.
Corollary 2 (Loose average spectral error bound of theapproximate SVT with power iteration). Theorem 4 is boundedas
EkSτ (A) − Sτ (Âk )k2≤ [Bound of Theorem 4]≤ C · poly(v) · σk+1 − G0 (A),1hi 2η+1q√√kwhere v={k, p, η}, poly(v)= 1 + p−1h−k+ e k+p·,ph=min(m, n), 1<C≤2 and G0 (A) = min(σk+1 (A), τ ) ≥ 0.
P
12(2η+1) 2
Proof. Since σj ≥ σj+1 for any j , the termj>k σj√2η+1
in Eq. (12) is upper-bounded by h − k · σk+1 . Then, reexpressing the inequality leads the conclusion.
Consequently, both of Theorems 3 and 4 assert thatthe bounds of the approximate SVT are tighter than therespective errors by rank-k approximation in Theorems 1 and
2 on average. Thus, we have the following properties from
Theorems 3 and 4:• Once σk+1 (A) approaches a small value, thenkSτ (A) − Sτ (Âk )k may approach close to zero, i.e., whichindicates the approximation is almost exact, roughly whenrank(A) ≤ k . Detailed analyses of how we can exploit thisproperty in the algorithm are described in Sec. 3.2.
• Fortunately, the fact that the bound becomes rapidly tighteras p or η increases in both practice and theory remains
consistent with the result of Halko et al. [12].
• The bound is independent of the size of the input matrix.
3.2 Residual Singular Value Check without Singular Value
Computation

In the previous section, we state that our algorithm may produce almost exact results when k ≥ rank(A), i.e., σk+1 = 0. We
may have a broader condition for ensuring accurate estimates,
because the bounds in Theorems 3 and 4 are tighter than thoseX X
for low-rank approximation in Theorems 1 and 2, respectively.
22≤ poly(v) ·σ (A) −
min(σj (A), τ ) ,
j>k j
j>k
In this section, we explain the condition in which the proposed
method can be expected to be accurate, and how the condition(by Theorem 1)
X
can be verified in a computationally inexpensive manner.
2= poly(v) ·σ (A) − G(A)j>k jGiven the target rank k and an approximated matrix Â
by the range approximation, if σk (Â) ≤ τ , then no accuracy
Obviously, G(A) ≥ 0 by definition.
improvement is possible. This is because all the remaining
As aforementioned, poly(v) has a simple form of Eq. (9)singular values {σj (Â)}j≥k are zeroed out by thresholding
only when power iteration is not taken into account. For theoperation with τ , and these does not affect the solution any
case with power iteration, it is useful to see the spectral error
further. Unless σk (Â) ≤ τ , we can additionally estimate abound of the approximate SVT, which is our main result.
few more bases as described in Sec. 2.1, and we may check
Theorem 4 (Average spectral error bound of the approximate the condition again until a desired accuracy is achieved. This
SVT with power iteration). Let Sτ (·) be the SVT operator and procedure can be adopted with Gram-Schmidt procedure [9] inÂk holds Theorem 2. The average spectral error satisfies the following Line 10 of Alg. 1 when an accurate estimate is required.
inequality:
While the complexity of the single additional basis estimation
only affects the factor s in Table 2 by s ← s + 1, which
EkSτ (A) − Sτ (Âk )k2 ≤still
preserves
the original complexity, checking σk (Â) requires1! 12  2η+1q
additional
computation
because it is unavailable until the end√
P 2(2η+1)2η+1k − G0 (A), of the FRSVT algorithm. Instead of directly estimating σ (Â),σk+1
C  1 + p−1+ e k+pσjkpj>k
we can utilize the following useful lemma.
(12)
EkSτ (A) − Sτ (Âk )k2F
X
≤ EkA − Âk k2F −min(σj (A), τ )2j>k

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR
m×n

Lemma 3 ( [42]). Let A ∈ R, and let q be a positive integer, αbe a real number greater than 1. By drawing an independent family{ω (i) : i = 1, 2, ..., q} of standard Gaussian vectors, then,qkAk2 ≤ α π2 max kAω (i) k2(13)i=1,...qwith probability at least 1 − α−q.

of Drineas et al. and Halko et al. [12], i.e., the over-sampling rate
to be the same as the specified target rank k .5
Corollary 3. Suppose that the sampling rate√c is set twice larger
than the target rank k , i.e., c = 2k , then  ≥ 2, and LTSVD also
holds the following bounds:X√ Xmin(m,n) 2
EkA − Âk k2F ≤σj2 (A) + 2
σj (A), (15)j>kThe lemma illustrates the relationship between the truelargest singular value and the estimate by Gaussian random
projection. From this, we derive the following proposition.
Proposition 4 (Residual singular value estimate). Let Q be any
rank-k orthonormal column matrix. Given Q, the (k+1)-th singularvalue is bounded by σk+1 (A) ≤ k(I − QQ∗ )Ak2 . Then, drawing
a standard Gaussian vector ω ,qk(I − QQ∗ )Ak2 ≤ α π2 k(I − QQ∗ )Aωk2 ,(14)with probability at least 1 − α−1 .
Proof. From the well known variational characterization of
singular values, we can obtain the minimum of k(I−QQ∗ )Ak2
with Q∗ Q = I (orthonormal column matrix) iff Q = Uk ,where Uk is the singular vector matrix of A corresponding
to the top-k largest singular values. Thus,k(I − QQ∗ )Ak2 ≥ k(I − UU∗ )Ak2 = σk+1 (A).
By drawing a single standard Gaussian random vector ω , theleft side of the above inequality can be further upper-boundedby Lemma 3 as Eq. (14) with probability at least 1 − α−1 . Thisconcludes the proof.
Since y = Aω in Eq. (14) is the same by-productwith the over-sampling during range propagation procedurein Sec. 2.1, we can check the provable upper bound estimate of (kq + 1)-th singular value σk+1 (Â) (≤ σk+1 (A)) byςk+1 = α π2 ky − QQ∗ yk2 with high probability (say 95%
with α = 20), where the residual vector computation y−QQ∗ yis also a by-product of the Gram-Schmidt procedure. Noticethat all the computation for estimating the quantity ςk+1 has
been already performed in the original FRSVT procedure;
therefore, no additional
computation is required except forq2
a scaling with α π , which is almost negligible. Thus, bychecking ςk+1 ≤ τ , we are able to avoid unnecessary additional
sampling.
3.3

Bound Comparison with Linear-Time SVD

We show the error bound comparison with a relevant methodby showing the approximation error bound of Drineas et al. [6].
Theorem 5 (Average error bound of LTSVD [6]). Let
A ∈ Rm×n and Hk be a matrix created by the LTSVD algorithm [6]by sampling c columns of A with probabilities {pi }ni=1 such thatpi ≥ β|A(i) |2 /kAk2F for some positive β ≤ 1 and  > 0. If
c ≥ 4k/β2 , then the average Frobenius error is bounded by
2
2
2
2
EkA − Hk H>
k AkF = EkA − Âk kF ≤ kA − Ak kF + kAkF .

In addition, the average spectral error is bounded byq
EkA − Hk H>
Ak
=
E
kA
−
Â
k
≤
kA − Ak k22 + kAk2F .
2k 2
kFor a fair comparison and simplicity, we set the same parameters when comparing the low-rank approximation bounds

8
EkA − Âk k2 ≤

j=1

 21
√ Xmin(m,n) 2
2
σj (A) . (16)σk+1 (A) + 2
j=1

Proof. From Theorem
5, we have β ≤ 1,  > 0 and c ≥ 4k/β2 .
qkHence,  ≥ 2 βc
. The lowest bound can be achieved at β = 1;q√therefore, if c = 2k is assumed,  ≥ 2 kc = 2. In addition,
we can obtain Eqs. 15 and 16 from the following bounds byusing the above result with Theorem 5.
√
EkA − Âk k2F ≤ kA − Ak k2F + 2kAk2F .
The average spectral bound can be derived analogously.
From Corollary 3 and Theorem 1, we obtain the following
conclusion.
Proposition 5. Suppose the target rank k ≥ 4 and the oversampling rate to be k (for Halko et al. [12], p = k , and for LTSVD,c = 2k ), then on average Halko et al. [12] without power iteration
has a tighter bound than LTSVD.
Proof. We compare the right bound terms of Theorem 1 andCorollary 3 in Eq. (15).
√ Pmin(m,n) 2
P
2
2 j=1
σj (A). (17)[LTSVD]j>k σj (A) +P[Halko et al. ]poly(v) · j>k σj2 (A).
(18)
Suppose that Eq. (17) < Eq. (18), and given Eq. (9),
P
j>k

⇒
⇒

σj2 (A) +√ min(m,n)
P
P 2
σj2 (A) < poly(v) ·σj (A)
2
j=1

j>k√ min(m,n)
P
σj2 (A) + 2
σj2 (A) < 0

k
− k−1
j=1
j>k
√ P
P 2
2
2
σj (A) + C(k) ·σj (A) < 0
1≤j≤k
j>k


√
k
C(k) = 2 − k−1

P

(19)

For k √
≥ 4, the minimum of C(k) is achieved at k = 4 asC(4) = 2 − 34 > 0. Thus, since C(k) > 0 for k ≥ 4, theleft term of Eq. (19) is positive for k ≥ 4. This contradicts
the assumption that Eq. (17) < Eq. (18), which concludes theproof.
For simplicity, we make a reasonable assumption, the target
rank k ≥ 4. Since the bound comparison is performed without
the power iteration scheme, Halko et al. [12] may have a loosebound. However, as shown in Proposition 5, Halko et al. [12]still result in a better approximation bound than LTSVD evenunder a non-ideal condition. Accompanied with the power
iteration, we may be able to derive broader conditions that
Halko et al. [12] perform better than LTSVD in a similar mannerto the proof process of Proposition 5. We would like note thatLTSVD is still fascinating in that it does not require dense
matrix multiplications for Gaussian matrix sketching used in
Halko et al. [12] by computing the approximate leveraging
score {pi } which is computationally less expensive.
5

It can be generalized with some efforts.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

−4

500

1000

1500

2000

10

0

Target rank

20

1500

Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP

−1

10

15
10

−2

10

−3

10

5
0
0

−4

500

1000

1500

Target rank(g) Computation time
Full rank 2000 × 2000.

2000

10

0

500

1000

1500

2000

Target rank(h) Accuracy
Full rank 2000 × 2000.

−3

500

1000

10

4

500

Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP

−1

10

−2

10

1
−3

1000

10

Target rank(i) Computation time
Full rank 2000 × 1000.

0

500

1000

Target rank(j) Accuracy
Full rank 2000 × 1000.

LTSVD
LANCZOSFRSVTFRSVT−RP

−10

10

20
−15

500

1000

1500

2000

10

50
40

0

500

1000

1500

2000

Target rank(e) Computation time
Low rank 2000 × 2000
60

10

2

500

40

10

Target rank

0

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

−5

60

0
0

1000

(d) Accuracy
Full rank 2000 × 1000.

3

0
0

80

0

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

Target rank(c) Computation time
Full rank 2000 × 1000

5

0

100

Accuracy

10

Computation Time (sec)

Accuracy−2

Target rank

6

0

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

10

5
0
0

2000

(b) Accuracy
Full rank 2000 × 2000.

Accuracy

Computation Time (sec)

25

1000

−1

10

Target rank(a) Computation time
Full rank 2000 × 2000
30

500

120
Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP(f) Accuracy
Low rank 2000 × 2000.
0

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP10−5

Accuracy

0
0

15

0

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

Computation Time (sec)

−3

10

20

20

Accuracy−210Computation Time (sec)40−1Computation Time (sec)60

Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

Accuracy

Computation Time (sec)

10

100
80

25

0

120

9

30
20

10

LTSVD
LANCZOSFRSVTFRSVT−RP

−10

10

10
0
0

−15

500

1000

1500

Target rank(k) Computation time
Low rank 2000 × 2000.

2000

10

0

500

1000

1500

2000

Target rank(l) Accuracy
Low rank 2000 × 2000.

Fig. 4: SVT comparisons among SVD methods. [Top] Experiments on Matlab 2010a. [Bottom] Experiments on Matlab 2014a. Accuracy is measured
by normalized root mean squared error (NRMSE), kA∗ − Âk kF /kA∗ kF , where A∗ = Sτ (DGT ), Âk = Sτ (Dk ), DGT is the input data, and Dk isapproximated by each method. For the low-rank matrix test, we generate input data matrices whose rank correspond to the target rank by multiplying two
Gaussian random matrices with m × r and r × n, while in other tests
P full rank matrices are used. In (b,d,h,j), the theoretical minimum error bound by rank
truncation in Sτ (Dk ) is provided for guidance, and it is defined as j>k σi (A∗ ), where σi (A∗ ) is computed by Matlab built-in SVD. For the low-rank casein (f,l), the theoretical minimum error is zero, so we omit the theoretical bound.

4

E XPERIMENTAL R ESULTS

In this section, we first evaluate the efficiency of the proposed
method in comparison with other methods using simulation
data. The evaluation is conducted by examining the performance of a single SVT computation and also by RPCA computation, which is arguably the most relevant NNM problem
in computer vision today. We then show NNM applications
in computer vision using real-world data; subspace clustering,semi-online weather artifact removal, and simultaneous multiimage alignment and rectification. All the experiments are conducted on a PC with Intel i7-3.4GHz and 16GB RAM. The same
shared parameters were used among different algorithms.
4.1

Evaluation using Simulation Data

We quantitatively evaluate our method in comparison to othermethods with synthetic matrices sampled from a standard
Gaussian distribution. Computational times are evaluated using Matlab 2010a and 2014a 64bits. Since the recent Matlab
has been intensively optimized on Intel CPU (mainly due toimprovement of Intel MKL), the computation efficiency has
been noticeably improved. Therefore, it is worth reporting
the performance difference on these two versions of Matlab, because most related works have been assessed with
Matlab older than 2011a [3], [20], [22], [23], [29]. For a faircomparison, we turn off multi-threading functions including
maxNumCompThreads(1) in Matlab. We describe the implementation details of each method in Table 3.
Single SVT test Figure 4 compares speed and accuracy of
SVT computation using various SVT methods, such as Matlab
built-in SVD6 (baseline), Lanczos [17], FSVT [3], LTSVD [6] andour FRSVT (with / without range propagation (RP hereafter)).
All the implementation details are summarized in Table 3.
Except for the baseline SVD, the others produce the truncatedSVD of the input matrix, and it is used for SVT computation.
6

We apply either the economic size or full SVD depending on the matrix
shape and report the faster one, but we denote either as econSVD.

We test using Gaussian random matrix
pdrawn from N (0, 1)
for the operation Sτ [·] with τ = 1/ kAk2 . For a rankk approximation, we compute rank-(k + p) approximations
by setting the over-sampling rate p to 2 for Lanczos, k for
LTSVD, 5 for FRSVT.7 We apply the power iteration twice, i.e.,η = 2. While LTSVD in Fig. 4-(g), (i), (l) is faster than ours,
the approximation error is significantly higher than any othertruncated SVDs, and it results in slower convergence in RPCA
as we will see in the following test.
Robust PCA test To observe the convergence behavior of
SVT methods in RPCA [4], we compare our method withvarious SVT methods, such as SVD, LTSVD, BLWS [20],
FSVT, RSVD8 , using an inexact augmented Lagrange multipliermethod [18] (iALM, or called alternating directional multiplier method), and other variants, such as Mu et al. [29]9 ,Active Subspace [22], L1 filtering [23] in Table 5 for resultant
performance. Also, the continuous convergence behaviors areshown in Figs. 5 (over the number of iteration) and 6 (overcomputation times).10 We generate data by following [35] with
10% rank ratio and 5% corruption ratio. We set the tradeas suggested by Candès et
off parameter λ = √ 1
max(m,n)al. [4] for all the methods. All other parameters for FRSVT
are same as mentioned in Sec. 2 and with Table 4. We apply
the proposed adaptive rank prediction described in Sec. 2.3
to LTSVD, RSVD and our FRSVT. LTSVD shows convergence
7

The rationals behind the used parameters are that Lanczos is an almost
exact method, LTSVD produces quite low accuracy despite high oversampling rate, and FRSVT is an approximate method fairly more accurate
than LTSVD even with low over-sampling rate.
8
As noted in Table 3, this implementation involves 2-side random projection and respective power iterations. By this, we can see the effects of thedifferent implementation of [12].
9For Mu et al. [29], we could not reproduce their result. However, we
include the results for thoroughness. The results are consistent with theearly reported results in Liu et al. [23].
10
The convergence of L1 filtering is not compared, because the algorithm
solves RPCA by dividing several block matrices of which convergence is
not comparable to other methods.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

−2

−6

−4

10

−6

10

−8

10

20

30

40

50

60

70

10

−4

10

−2

10

20

30

40

50

60

70

10

−6

−8

0

10

20

Iterations(a) 2000 × 2000

−4

10

10

−8

0

Iterations

ALM−econSVD
ALM−LTSVD
ALM−BLWS
Mu et. al.
ALM−FSVTALM−RSVD
ALM−FRSVTALM−FRSVT−RP

10

10

−8

0

−2

−6

10

Convergence, 10000x100

010
ALM−econSVD
ALM−LTSVD
ALM−BLWS
Mu et. al.
ALM−FSVTALM−RSVD
ALM−FRSVTALM−FRSVT−RP

10

Convergence

Convergence−4

10

Convergence, 6000x6000

010
ALM−econSVD
ALM−LTSVD
ALM−BLWS
Mu et. al.
ALM−FSVTALM−RSVD
ALM−FRSVTALM−FRSVT−RP

10

Convergence

−2

10

10

Convergence, 4000x4000

010
ALM−econSVD
ALM−LTSVD
ALM−BLWS
Mu et. al.
ALM−FSVTALM−RSVD
ALM−FRSVTALM−FRSVT−RP

Convergence

Convergence, 2000x2000

0

10

10

30

40

50

60

70

10

0

10

20

Iterations(b) 4000 × 4000

30

40

50

60

70

Iterations(c) 6000 × 6000(d) 10000 × 100

Fig. 5: RPCA comparisons. X-axis: Number of iterations, Y-axis: Stopping criterion. Except Mu et al. is based on the exact ALM method, all the othermethods are based on the inexact ALM.
Built-in SVD(·, ’econ’) is based on LAPACK 3.2.1 in IntelMKL 10.2.2.

SVDMatlab 2014a

Built-in SVD(·, ’econ’) is based on LAPACK 3.4.1 in IntelMKL 11.0.5.

LTSVD*The Matlab implementation by Ma et al. [26] is used.

Lanczos

The implementation optimized by Lin et al. [17] is used. It
is implemented by Matlab and Mex to use LAPACK (ver.
3.4.1 in Intel MKL 11.0.5 is used).

BLWS

The Matlab implementation by Lin et al. [20] is used.

Mu et al.

The Matlab implementation provided by Mu et al. [29] isused. All the parameters are directly used as suggested.

Active Subspace

The Matlab implementation provided by Liu et al. [22] isused. All the parameters are directly used as suggested.

L1 filtering

The Matlab implementation provided by Liu et al. [23] isused. All the parameters are directly used as suggested.

FSVT ***

10 0
ALM-LTSVD
ALM-BLWS
Mu et. al.
ALM-FSVT
ALM-ActiveSub
ALM-RSVD
ALM-FRSVT
ALM-FRSVT-WS

10 -2
Convergence

SVDMatlab 2010a

10 -4

10 -6

10 -8

0

100

200

300
400
Time (sec)

500

600

700

(a) Matlab 2010a
10 0

We implement RSVD [12] based on LAPACK 3.4.1 in IntelMKL 11.0.5. by referring to RedSVD**.

FRSVT ***

Our Matlab implementation based on Intel MKL 11.0.5. Partially, we implement Mex functions for QR–CP (dgeqp3)
and Polar Decomposition with BLAS and LAPACK 3.4.1
routines.

Convergence

10 -4

10 -6

We implement FSVT by Matlab and Mex based on IntelMKL 11.0.5, and the speed-up gain is checked and consistent with the reported results in Cai et al. [3]. All theparameters are directly used as suggested.

RSVD ***

ALM-econSVD
ALM-LTSVD
ALM-BLWS
Mu et. al.
ALM-FSVT
ALM-ActiveSub
ALM-RSVD
ALM-FRSVT
ALM-FRSVT-WS

10 -2

10 -8

0

100

200

300
400
Time (sec)

500

600

700

(b) Matlab 2014aFig. 6: Convergence time plots of RPCA. X-axis: Elapsed time (sec.), Y-axis:
Stopping criterion. Except Mu et al. is based on the exact ALM method,
all the other methods are based on the inexact ALM. As a representative
example, matrices of size 6000 × 6000 are used, but the results can begeneralized to other size matrices.

* LTSVD: In the implementation of [26], the leverage score estimation suggested by
the original paper [6] is omitted, which is required for the importance sampling. So
we implement it with O(mn) complexity.
** RedSVD (http://code.google.com/p/redsvd/): Since RedSVD is implemented
based on Eigen 3.0-b1, without the power iteration but with 2-side random projection,
so the original implementation is not fast and inaccurate. Thus, we re-implement
RSVD [12] by referring to RedSVD with 2-side power iteration.
*** For FSVT, RSVD and FRSVT implementation, we utilize Armadillo library [38] as ahigh level wrapper for LAPACK and BLAS of Intel MKL in our Mex implementation
part. In the Mex functions of Lanczos, Lin et al. [17] directly call LAPACK routines
without any external library.

η = 2 without adopting AP. Without AP, even in the case of
p = 2, which has been used as our standard experiment setting,iALMs with both FRSVT and RSVD do not converge to theground truth. This indicates that AP is necessary for a fixedover-sampling rate scheme. Fig. 7-(b) shows the results withthe varying number of power iterations η , where we fix p = 2
with AP. The result implies that a single power iteration is
sufficient in practice.

TABLE 3: Lists of comparison methods and details for SVT tests. All theparallelization techniques are turned off for fair comparison.

4.2

iALMMethodspAPηeconSVD
–
×
–

LTSVD
+5

BLWS
–

–

–

FSVT
–
×
–

RSVD
+2

FRSVT
+2

FRSVT-RP
+2

2

2

2

TABLE 4: Parameter settings for convergence evaluation of RPCA (p: Oversampling rate, AP: Adaptive rank prediction, η : The number of power
iterations). Mu et al. [29], Active Subspace [22] and L1 -filtering [23] are
not compatible with the parameters.

degradation when achieving comparable accuracy with others
due to rough approximation on both bases and singular values.
Other methods including our FRSVT show similar numbersof iterations and accuracy, but our method has a considerably
lower computation time for a single iteration.
Fig. 7 shows iALM behavior over varying over-sampling or
power iteration settings. In Fig. 7-(a), varying over-sampling
rate p is tested, where we fix the number of power iterations

Applications

This section shows applications of our FRSVT method tovarious types of low-rank optimization problems summarizedin Table 6, such as affine constrained NNM (Type A), nonconvex truncated NNM (Type B) and NNM on tensor structure(Type C). We show them with typical computer vision applications in what follows.
Type A - Robust Subspace Clustering by Low-Rank Representation (LRR) Many visual data are often characterized bya mixture of multiple subspaces. One of the recent promising
methods is LRR, which effectively performs subspace clustering and noise correction simultaneously. While both noise correction and subspace clustering are known to be challenging,robust LRR has shown effective even with large corruptions.
The robust LRR can be formulated by Type A, which can be
efficiently solved by iALM. In this experiment, we use the same
parameter settings and evaluation metrics suggested in [21].

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR
iALMeconSVD
Matlab 2010a
Matlab 2014aiALMLTSVDiALMBLWS

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

99,999399,997
1,599,981–
3,599,952

879,415
994,479
5,607,211
35,662,548

99,999399,998
1,599,984
3,599,966

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

7.34e-05
2.11e-07
1.81e-07–

1.42e-07

3.39e-02
7.77e-07
4.59e-07
7.76e-05

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

3.217
2751.910
21668.280–

1.982
147.055
1068.863
3360.379

1.722
11.07569.547203.051

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

0.140
119.415947.908
–

0.086
6.393
46.481
146.192

0.039
0.274
1.726
4.295

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

–

23

43
41
41
48

10000 × 100
2000 × 2000
4000 × 4000
6000 × 6000

-

1.62×
18.71×20.27×-1.87×248.48×
311.56×
-

23
23
23

Mu et al.

11
iALMFSVTiALMRSVD

iALMFRSVT-RP(Proposed)

Active Sub.

L1 filter

99,999399,997
1,599,981
3,599,952

132,666
399,991
1,599,968
3,599,938

231,604
402,932
1,663,440
3,752,782

9.33e-05
2.11e-07
1.81e-07
1.43e-07

3.00e-02
5.70e-07
3.56e-07
2.89e-07

5.16e-02
3.61e-07
2.80e-07
2.08e-07

0.782
5.123
30.27488.132

1.349
11.16379.471
246.738

3.179
94.300683.159
1914.431

0.034
0.226
1.348
3.942

0.047
0.465
3.311
10.281

-

23
23
23
23

29
24
24
24

-

4.11×
537.17×715.74×-

2.38 ×246.52 ×
272.66 ×-

1.01 ×29.18 ×31.72 ×-kSk0
999,999
99,99999,9993,998,401
399,997
399,993
15,995,318
1,599,981
1,599,979
35,990,553
3,599,952
3,599,955
Accuracy
2.66e-04
7.65e-01
7.34e-05
6.75e-05
1.03e-06
1.10e+00
2.11e-07
3.06e-07
4.83e-07
1.37e+00
1.81e-07
2.08e-07
3.83e-07
1.49e+00
1.42e-07
1.71e-07
Total Time
1.317
125.368
2.436
0.863
9.808
118.798348.132
7.928
53.468759.794
2527.14149.027164.752
2385.272
8238.109147.139
Avg. time for a single iteration
0.049
1.687
0.104
0.038
0.393
1.667
15.1890.345
2.185
10.664110.006
2.132
6.848
33.505358.899
6.397
Total no. iteration
24
73
23
23
24
70
23
23
24
70
23
23
24
70
23
23
Speed-up gain against the baseline
2.44×
0.03×
1.32×
3.73×
280.58×23.16×7.90×347.11×405.26×28.52×8.57×
441.97×-

TABLE 5: Quantitative comparison on RPCA. The accuracy is measured by NRMSE with ground-truth data, and the speed-up gain is computed against
the baseline method, iALMeconSVD on Matlab 2010a. For Mu et al. [29], we used the implementation and parameters provided by the authors. Since Mu etal. is based on the exact ALM algorithm [18], while other methods are based on the inexact ALM [18], we report the times for a single outer iteration of
Mu et al. as ‘Avg. time for a single iteration’. In addition, the highly stochastic property of the randomly projected nuclear norm makes the algorithm [29]sensitive to data and does not preserve the singular values and sparsity structure of the matrix. This result is consistent with another stochastic method,
LTSVD (see and compare kSk0 of LTSVD and Mu et al.). Also, since L1 filtering is incompatible to compare with other methods, the time for a single
iteration and the total number of iteration are not reported.

NRMSE

10 0

10

-2

10

-4

10

-6

Type
A
iALM-econSVD
iALM-FRSVT-WS [p=0]iALM-FRSVT-WS [p=2]iALM-FRSVT-WS [p=6]iALM-FRSVT-WS [p=10]
iALM-RSVD [p=0]iALM-RSVD [p=2]iALM-RSVD [p=6]iALM-RSVD [p=10]

0

5

10

15

20

25

30

35

B
C

40

(a) Over-sampling

NRMSE

argminiALM-econSVD
iALM-FRSVT-WS [η=0]iALM-FRSVT-WS [η=1]
iALM-FRSVT-WS [η=2]
iALM-RSVD [η=0]iALM-RSVD [η=1]
iALM-RSVD [η=2]

5

10

15

20

25

30

35

40

Iterations(b) Power iter.

Fig. 7: Varying parameter tests for (a) over-sampling and (b) power iterationin RPCA problem. X-axis: Number of iterations, Y-axis: NRMSE of low-rankkL−Lksolution L as kL GTk F .
GT

O=L+S

F

We apply FRSVT to the robust LRR for motion segmentation on Hopkins155 [40] and for face recognition on a partof Extended Yale Database B [16], respectively. All parametersfor FRSVT are the same as mentioned in Sec. 2. We use all
the 156 sequences in Hopkins155, and the first 10 classes inExtended Yale Database B, in which each class contains 64
face images captured under various illumination conditions –given 42 × 48 resized images, we construct a 2016 × 640 data
matrix by vectorizing each image. While Hopkins155 is only
weakly corrupted, Yale data is heavily corrupted by shadows,
specularity and noise. As shown in Table 7, our method speeds
up LRR without degrading of the accuracy.

αi kL(i) k∗ + λkEk1

O◦Γ=L+E

TABLE 6: Examples of NNM related objective functions. (A) Low-rank representation (LRR). (B) RPCA based on the non-convex truncated nuclear
norm, where k is the target rank. (C) Low-rank and sparse 3-order tensor
decomposition with alignment. Here, k·k2,1 is l2,1 norm,
P {αi } are balancing
parameters among the unfolding matrices L(i) , andαi = 1 is assumed.
We refer to the basic tensor algebra notations described in [44].

10 -6
0

3
P

L,E,Γ i=1

10 0

10 -4

Constraint
O = ZL + S

L,S

Iterations

10 -2

Objective function
argmin kLk∗ + λkSk2,1
L,S
P
argmin ni=k+1 σi (L) + λkSk1

Time per Motion
Time for 640 Faces

Computational Time (s)LRR+SVDLRR+SVD
Matlab 2010a
Matlab 2014a
1.230
1.016
419.440
75.216

LRR+Ours
0.452
44.590TABLE 7: Comparisons of the subspace segmentation algorithms on Hopkins155 [40] (Motion) and Extended Yale Database B [16] (Face). Motion
Segmentation Errors on the Hopkins155 of both LRR and LRR+Ours are
1.59%. The segmentation accuracies (%) on the Yale data are 79.06 forboth LRR and LRR+Ours. These results show the improvement of thecomputation time with retaining the same accuracy.

SVT

Update Z

Update E

Others

LRR

Ours

0

20

40

60

80

Computational Time (Seconds)

Fig. 8: Computation time comparison on the subspace clustering application [21]. Computation times are measured for face clustering on Extended
Yale Database B [16].

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

(a)

(b)

Fig. 9: Qualitative results for the weather artifact removal. (a) Sample input
images O. (b) Low-rank images Z.

Type B - Semi-Online Weather Artifact Removal We consider a weather artifact (e.g., snow or rain) removal problem
in a video sequence. Since the weather artifacts have nondeterministic appearance and a sparse yet random distribution
in the spatio-temporal domain, we model the artifacts as sparse
outlier and the scene as low-rank, while we want to retain
moving objects. We apply the low-rank and sparsity decomposition to n frames in a sliding window manner to leave moving
objects in the latent images.
Unfortunately, since finding a low-rank solution by the
nuclear norm relies on the high dimensionality of data [4],it could be degraded when only with a small number offrames as reported in Oh et al. [30]. With the assumption thatthe object motions are small within a few frames, we canencourage the low-rank solution to be rank-1 and decompose
sparse corruptions by Type B optimization with k = 1, where
O is constructed by stacking vectorized n images. It can beeffectively solved by alternating Partial SVT (PSVT) [30] andl1 minimization based on iALM. We replace PSVT by FRSVT
with the partial thresholding. We transfer the previous basis
estimation to the next n frame optimization, so it can be
regarded as a semi-online algorithm.
We set the trade-off parameter λ = √ 1
as suggested
max(m,n)by Oh et al. [30], and limit the maximum rank-l for computing
partial SVD in FRSVT as: l = 2 when n = 3, and l = 3 whenn = 5 and 10. We also apply the power iteration twice, i.e.,η = 2. Our algorithm produces n results simultaneously for nimages. With an n = 5 sliding window, the method based on
SVD takes 158.9ms per a single channel of a 384×288 image
on Matlab 2014a (318.3ms on Matlab 2010a), while our method
only takes 65.5ms. The qualitative results of our method can befound in Figs. 9 and 10, which show plausible snow removal.
Type C - Simultaneous Multi-Image Alignment and Rectification Simultaneous multi-image alignment and rectification
problem is Type C optimization suggested by Zhang et al. [44](called SRALT). The method combines the ideas of TILT [45]and RASL [44] on the third-order tensor structure, which
exploits the low-rank texture property and nuclear norm-based
misalignment error. We show the applicability of FRSVT to
NNM on the third-order tensor. We first show that SRALT withFRSVT produces consistent results with the original SRALT on
Windows dataset [36] in Fig. 12. Then, we present a new application of SRALT using gait data. We use the parameters α and0.5, α = [0.1, 0.1, 0.8]
λ in Type-C of Table 6 as: λ = √width·heightfor Windows, and α = [0.15, 0.15, 0.7] for the gait data.
Since gait is a biometric signal spanned on not only spatial
domain (2D), but also temporal domain, so it is natural to rep-12

resent the data by a third–order tensor [25]. In gait recognition,
the exact alignment of acquired data is frequently assumed,
but in real situations, the accuracy of the alignment varies
depending on pivot angles of cameras and locations of moving
people. Therefore, the SRALT framework is useful to aligngait data making sure that the assumption holds and to finda vertical angle as a preprocessing step. We observe that thehuman silhouettes have the low-rank texture property, and as
observed in Lu et al. [25], cyclic gait motions are spanned by a
few number of bases.
We conduct the experiment using Gait Challenge data
set [39], called USF Human ID version 1.7. There are 731 gait
samples, each with ten sampled gait cycles, and we use thefirst 300 gait samples. We resize the images into 32 × 22, so thesize of the input tensor is O ∈ R32×22×3000 (3000 unaligned
images with the size 22 × 32). The synthetic misalignment
generation parameters are described in the below. Randomly
drawn Euclidean (3 DoF) geometric transformations are applied to each gait image. The parameters are set to: 1) Scale
5π 2factor: fixed to 1, 2) Angle noise: N ( 10π
180 , ( 180 ) ), 3) Translation2
noise: N (0, 0.5 )+U (−1, 1), where N (µ, σ 2 ) denotes Gaussian
distribution with the mean µ and variance σ 2 , and U (a, b)
denotes uniform distribution defined on [a, b].
The set of geometric parameters Γ = {g1 , · · · , g3000 } (seeType C in Table 6) are defined with a p-group parameterizationas gj ∈ Rp , where 3-DoF parameterization (i.e., p = 3) for
Euclidean transformation is used in our experiments. Our
method is implemented on top of [44] by replacing their SVT
method with FRSVT. While the method of [44] converges
at the objective value 206.69 and takes about 6 hours 51
min. (24652s), our method takes only about 1 hour 45 min.
(6280s) and achieves an even smaller objective score 206.46.
The qualitative and detailed computation time comparisons are
summarized in Figs. 13 and 14, respectively.

5

D ISCUSSIONS

We have presented a fast approximate SVT method that exploits the property of iterative NNM procedures by range
propagation and adaptive rank prediction. The approximation
error bound shows that our method can produce reliableapproximation. The proposed method has been assessed using
the problems of affine constrained NNM, non-convex NNM
and NNM on tensor structure as well as the original NNM.
The empirical evaluations showed the consistent result withthe theoretical analysis, and our approach can reduce thecomputational time of applications, which involve low-rank
optimization, without degrading accuracy and convergence
behavior.
For convergence, while a general convergence analysis
against various host problems and algorithms is very interesting, it is going to be challenging. RPCA and the TypeA application are convex programs, thus if the algorithms
converge with satisfying the respective termination criteria(derived from the convergence criteria [18], [19]), then thealgorithms may produce a solution close to global optimal
within a small distance (i.e., up to -optimality). At this point,
we can only say that a host algorithm integrated with FRSVT
converges to -optimal, if and only if the problem is convex,
the host algorithm converges and FRSVT yields reasonably
accurate solutions.
The current major computational bottleneck of our method
is in the power iteration scheme. We are interested in further

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

(a)

(b)

(c)

13

(d)

(e)

Fig. 10: Comparisons according to varying sliding window lengths. (a) Sampled input image. (b) The exact PSVT based result with the sliding window
length 5. (c-e) The PSVT by FRSVT based results with the sliding window length {3, 5, 10} respectively. [Top] Low-rank images, L. [Bottom] Absolute ofsparse images, abs(S).
Computational Time (ms)Sliding Window n
35
Elapsed Time71.5 65.5

7
68.2

10
66.2

Fig. 11: Computation times of the FRSVT based semi-online weather artifact removal. The elapsed times are measured for a single color channel.
Since our algorithm produces n results simultaneously for n images at an
execution, we report the time for a single image for a single channel. This
results show that, when the number of inputs is small, our weather artifact
removal algorithm has linear complexity according to the length of a sliding
window, thus the elapsed times for a single image are similar.

(a) Input(b) SRALT [44](c) SRALT + FRSVT

Fig. 12: Qualitative comparisons for the simultaneously alignment and rectification on Windows data [36]. The final solutions produce similar results
while the computational speed is enhanced by FRSVT.

reducing the computational complexity by effectively relaxing
this computation block in the future.

ACKNOWLEDGMENTS
The first and fourth authors were supported by the NRF of
Korea grant funded by the Korea government, MSIP (No. 20100028680). The second author was partly supported by JSPS
KAKENHI Grant Number JP16H01732.

R EFERENCES
[1]
[2]
[3]
[4]
[5]
[6]

S. Boyd and L. Vandenberghe. Convex optimization. Cambridge
university press, 2004.
J.-F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding
algorithm for matrix completion. SIAM Journal on Optimization, 2010.
J.-F. Cai and S. Osher. Fast singular value thresholding without
singular value decomposition. Methods and Applications of Analysis,
2013.
E. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component
analysis? Journal of the ACM, 2011.
Y.-L. Chen and C.-T. Hsu. A generalized low-rank appearance model
for spatio-temporally correlated rain streaks. In IEEE International
Conference on Computer Vision (ICCV), 2013.
P. Drineas, R. Kannan, and M. W. Hahoney. Fast monte carloalgorithms for matrices II: Computing a low-rank approximation to amatrix. SIAM Journal of Computing, 2006.

[7]

A. Eriksson and A. van den Hengel. Efficient computation of robust
weighted low-rank matrix approximations using the l1 norm. IEEE
Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 2012.
[8] A. Ganesh, Z. Lin, J. Wright, L. Wu, M. Chen, and Y. Ma. Fast
algorithms for recovering a corrupted low-rank matrix. In IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive
Processing, 2009.
[9] G. H. Golub and C. F. V. Loan. Matrix computations 3rd Ed. Johns
Hopkins University Press, 1996.
[10] S. Gu, L. Zhang, W. Zuo, and X. Feng. Weighted nuclear norm
minimization with application to image denoising. In IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2014.
[11] E. T. Hale, W. Yin, and Y. Zhang. Fixed-point continuation forl1 -minimization: Methodology and convergence. SIAM Journal on
Optimization, 2008.
[12] N. Halko, P.-G. Martinsson, and J. A. Tropp. Finding structure withrandomness: Probabilistic algorithms for constructing approximate
matrix decompositions. SIAM Review, 2011.
[13] N. J. Higham. Computing the polar decomposition-with applications.
SIAM Journal on Scientific and Statistical Computing, 1986.
[14] Y. Hu, D. Zhang, J. Ye, X. Li, and X. He. Fast and accurate matrix completion via truncated nuclear norm regularization. IEEE Transactions
on Pattern Analysis and Machine Intelligence (TPAMI), 2013.
[15] S. Ji and J. Ye. An accelerated gradient method for trace norm
minimization. In International Conference on Machine Learning (ICML).
ACM, 2009.
[16] K.-C. Lee, J. Ho, and D. Kriegman. Acquiring linear subspaces forface recognition under variable lighting. IEEE Transactions on PatternAnalysis and Machine Intelligence (TPAMI), 2005.
[17] Z. Lin. Some software packages for partial svd computation. In
arXiv:1108.1548, 2011.
[18] Z. Lin, M. Chen, and Y. Ma. The augmented Lagrange multipliermethod for exact recovery of corrupted low-rank matrices. Technical
Report UILU-ENG-09-2215, UIUC, 2009.
[19] Z. Lin, R. Liu, and Z. Su. Linearized alternating direction method withadaptive penalty for low-rank representation. In Advances in Neural
Information Processing Systems (NIPS), 2011.
[20] Z. Lin and S. Wei. A block lanczos with warm start technique for accelerating nuclear norm minimization algorithms. In arXiv:1012.0365,
2010.
[21] G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma. Robust recovery ofsubspace structures by low-rank representation. IEEE Transactions on
Pattern Analysis and Machine Intelligence (TPAMI), 2013.
[22] G. Liu and S. Yan. Active subspace: Toward scalable low-ranklearning. Neural Computation, 2012.
[23] R. Liu, Z. Lin, Z. Su, and J. Gao. Linear time principal component
pursuit and its extensions using l1 filtering. Neurocomputing, 2014.
[24] C. Lu, J. Tang, S. Y. Yan, and Z. Lin. Generalized nonconvex nonsmooth low-rank minimization. In IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2014.
[25] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos. MPCA: Multilinear principal component analysis of tensor objects. IEEE Transactions
on Neural Networks, 2008.
[26] S. Ma, D. Goldfarb, and L. Chen. Fixed point and bregman iterative
methods for matrix rank minimization. Mathematical Programming,
2011.
[27] M. W. Mahoney. Randomized algorithms for matrices and data.
Foundations and Trends R in Machine Learning, 2011.
[28] K. Mohan and M. Fazel.
Iterative reweighted algorithms forIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, TO APPEAR

14

Fig. 13: Qualitative comparison of RASL [36] and SRALT [44] + Our FRSVT. [Top] Samples of input images. [Middle] Results obtained by RASL [36].
[Bottom] Results obtained by SRALT [44] + Our FRSVT. [Right side] Average images. (Top: for all input images, Middle: for RASL, Bottom: for SRALT+Ours).
Due to the angular bias shown in (a), RASL aligns samples to be consistent at the biased angle, while the SRALT based method in (c) not only aligns
images, but also corrects upright poses.
SVT1

SVT2

SVT3

[46] Y. Zheng, G. Liu, S. Sugimoto, S. Yan, and M. Okutomi. Practical lowrank matrix approximation under robust l1 -norm. In IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2012.

Others

SRALT

Ours

0

100

200

300

400

500

600

700

800

900

Computational Time (Seconds)

Fig. 14: Computation time comparison on registration and rectification
application [44]. The NNM for tensor consists of 3-way SVTs. Computationaltimes are measured for a single outer iteration (almost 100 inner iterations).
matrix rank minimization. Journal of Machine Learning Research,
13(Nov):3441–3473, 2012.
[29] Y. Mu, J. Dong, X. Yuan, and S. Yan. Accelerated low-rank visual
recovery by random projection. In IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2011.
[30] T.-H. Oh, H. Kim, Y.-W. Tai, J.-C. Bazin, and I. S. Kweon. Partial sum
minimization of singular values in RPCA for low-level vision. In IEEE
International Conference on Computer Vision (ICCV), 2013.
[31] T.-H. Oh, J.-Y. Lee, and I. S. Kweon. High dynamic range imaging by arank-1 constraint. In IEEE International Conference on Image Processing(ICIP), 2013.
[32] T.-H. Oh, J.-Y. Lee, Y.-W. Tai, and I. S. Kweon. Robust high dynamic
range imaging by rank minimization. IEEE Transactions on PatternAnalysis and Machine Intelligence (TPAMI), 2015.
[33] T.-H. Oh, Y. Matsushita, I. S. Kweon, and D. Wipf. A pseudo-bayesianalgorithm for robust PCA. In Advances in Neural Information Processing
Systems (NIPS), 2016.
[34] T.-H. Oh, Y. Matsushita, Y.-W. Tai, and I. S. Kweon. Fast randomized
singular value thresholding for nuclear norm minimization. In IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2015.
[35] T.-H. Oh, Y.-W. Tai, J.-C. Bazin, H. Kim, and I. S. Kweon. Partial
sum minimization of singular values in robust PCA: Algorithm
and applications. IEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI), 2015.
[36] Y. Peng, A. Ganesh, J. Wright, W. Xu, and Y. Ma. RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated
images. IEEE Transactions on Pattern Analysis and Machine Intelligence(TPAMI), 2012.
[37] G. Pierra. Decomposition through formalization in a product space.
Mathematical Programming, 1984.
[38] C. Sanderson. Armadillo: An open source c++ linear algebra libraryfor fast prototyping and computationally intensive experiments. Technical Report, NICTA, 2010.
[39] S. Sarkar, P. J. Phillips, Z. Liu, I. R. Vega, P. Grother, and K. W. Bowyer.
The humanid gait challenge problem: Data sets, performance, andanalysis. IEEE Transactions on Pattern Analysis and Machine Intelligence(TPAMI), 2005.
[40] R. Tron and R. Vidal. A benchmark for the comparison of 3-d motion
segmentation algorithms. In IEEE Conference on Computer Vision and
Pattern Recognition (CVPR), 2007.
[41] D. Wipf. Non-convex rank minimization via an empirical bayesian
approach. In Uncertainty in Artificial Intelligence (UAI), 2012.
[42] F. Woolfe, E. Liberty, V. Rokhlin, and M. Tygert. A fast randomized algorithm for the approximation of matrices. Applied and Computational
Harmonic Analysis, 25(3):335–366, 2008.
[43] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma. Robust
photometric stereo via low-rank matrix completion and recovery. In
Asian Conference on Computer Vision (ACCV), 2010.
[44] X. Zhang, D. Wang, Z. Zhou, and Y. Ma. Simultaneous rectification
and alignment via robust recovery of low-rank tensors. In Advancesin Neural Information Processing Systems (NIPS), 2013.
[45] Z. Zhang, A. Ganesh, X. Liang, and Y. Ma. TILT: Transform invariant
low-rank textures. International Journal of Computer Vision (IJCV), 2012.

Tae-Hyun Oh is a postdoctoral associate at
MIT/CSAIL. He received the B.E. degree (The highest ranking) in Computer Engineering from KwangWoon University in 2010, and the M.S. degree inElectrical Engineering from KAIST in 2012. He iscurrently working toward the Ph.D. degree at KAIST,South Korea. He was a visiting student in the VisualComputing group, Microsoft Research, Beijing andin the Cognitive group, Microsoft Research, Redmond in 2014 and 2016, respectively. He was arecipient of Microsoft Research Asia fellowship, agold prize of Samsung HumanTech thesis award, Qualcomm thesis award
and top research achievement awards from KAIST. His research interests
include robust computer vision and machine learning. He is a student
member of the IEEE.
Yasuyuki Matsushita received his B.S., M.S. and
Ph.D. degrees in EECS from the University of Tokyoin 1998, 2000, and 2003, respectively. From April
2003 to March 2015, he was with Visual Computing
group at Microsoft Research Asia. In April 2015, hejoined Osaka University as a professor. His researcharea includes computer vision, machine learningand optimization. He is on the editorial board ofIEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI), International Journal of Computer Vision (IJCV), IPSJ Journal of Computer Vision and Applications (CVA), The Visual Computer Journal, and Encyclopedia of Computer Vision. He served/is serving as a Program Co-Chair of
PSIVT 2010, 3DIMPVT 2011, ACCV 2012, ICCV 2017, and a General CoChair for ACCV 2014. He is a senior member of IEEE.
Yu-Wing Tai is a research director of youtu researchlab of Tencent SNG. He was a principle researchscientist of SenseTime Group Limited from September 2015 to December 2016. He was an associate professor working at KAIST from July 2009
to August 2015. He received his Ph.D. degree in
Computer Science from the National University ofSingapore (NUS) in 2009. From Sept 2007 to June
2008, he worked as a full-time student internshipin the Microsoft Research Asia (MSRA). He wasawarded the Microsoft Research Asia Fellowship in
2007, and the KAIST 40th Anniversary Academic Award for Excellent Professor in 2011 respectively. He was an area chair of ICCV 2011 and ICCV2015 respectively. He received a M.Phil and B.Eng (First Class Honors)
degree in Computer Science from the Hong Kong University of Science andTechnology (HKUST) in 2005 and 2003 respectively. He is a senior memberof the IEEE.
In So Kweon received the BS and MS degrees
in mechanical design and production engineeringfrom Seoul National University, in 1981 and 1983,
respectively, and the PhD degree in robotics from
Carnegie Mellon University, 1990. He is now a professor at KAIST, South Korea. He is a recipient of
the best student paper runner-up award at CVPR
09’. He was the program chair for ACCV 07’ and is
the general chair for ACCV 12’. He is also on theeditorial board of the IJCV.