Fast Randomized Singular Value Thresholding for Nuclear Norm Minimization
Tae-Hyun Oh∗
KAIST
thoh.kaist.ac.kr@gmail.com

Yasuyuki Matsushita
Microsoft Research Asia
/ Osaka University
yasumat@ist.osaka-u.ac.jp

Yu-Wing Tai
KAISTIn So Kweon
KAIST

yuwing@gmail.com

iskweon77@kaist.ac.kr

29, 28, 31, 36, 38]. As rank minimization using the l0 -norm
is an NP-hard problem, it is typically relaxed using the nuclear norm (i.e.,  · ∗ , the sum of all the singular values),
which is the convex envelope of the rank function.
As a simple example, a nuclear norm minimization(NNM) problem is expressed as:

Abstract
Rank minimization problem can be boiled down toeither Nuclear Norm Minimization (NNM) or Weighted
NNM (WNNM) problem. The problems related to NNM(or WNNM) can be solved iteratively by applying a closedform proximal operator, called Singular Value Thresholding (SVT) (or Weighted SVT), but they suffer from high computational cost to compute a Singular Value Decomposition (SVD) at each iteration. In this paper, we propose an
accurate and fast approximation method for SVT, called fast
randomized SVT (FRSVT), where we avoid direct computation of SVD. The key idea is to extract an approximate basis for the range of a matrix from its compressed matrix.
Given the basis, we compute the partial singular values ofthe original matrix from a small factored matrix. While thebasis approximation is the bottleneck, our method is alreadyseveralfold faster than thin SVD. By adopting a range propagation technique, we can further avoid one of the bottleneck at each iteration. Our theoretical analysis provides
a stepping stone between the approximation bound of SVD
and its effect to NNM via SVT. Along with the analysis, ourempirical results on both quantitative and qualitative studies show our approximation rarely harms the convergence
behavior of the host algorithms. We apply it and validate
the efﬁciency of our method on various vision problems,e.g. subspace clustering, weather artifact removal, simultaneous multi-image alignment and rectiﬁcation.

X∗ = argminX f (X) + τ X∗ ,(1)where X ∈ Rm×n , and τ > 0 is a regularization parameter.
The function f (X) can be deﬁned according to different
applications, e.g., f (X) = O − X1 in robust principal
component analysis (RPCA) [3], f (X) = 12 AX − B2F
in multivariate regression and multi-class learning [25], andf (X) = 12 πΨ (O) − πΨ (X)2F in matrix completion [3].
Here, O is measured data,  · 1 and  · F denote thel1 and the Frobenius norms, and πΨ (·) is an orthogonal
projection operator setting [πΨ (X)]i,j = [X]i,j for (i, j) ∈Ψ and 0 otherwise. Also, weighted nuclear norm minimization (WNNM), which can be non-convex according toweights [9, 13, 27], is used to better approximate rank(·).
Unless the loss function f (X) is a proximity term, mostprevious works use ﬁrst-order optimization algorithms,
e.g., dual method [7], accelerated proximal gradient [14],augmented Lagrange multiplier [17], alternating direction
method [17, 18], and iterative reweighted least squares [23].
All these approaches, NNM, WNNM and their regularizedversions, iteratively solve the simple NNM subproblem deﬁned as the following nuclear norm plus the proximity term:
Problem (Nuclear norm minimization). For τ ≥ 0 and
A ∈ Rm×n ,1. Introduction

1X∗ = argminX τ X∗ + X − A2F ,2

Low-rank matrix recovery problems arise in many engineering and applied science problems, and rank minimization techniques have attracted tremendous interests. Rank
minimization is a crucial regularizer to derive a low-ranksolution, and is required in many mathematical models incomputer vision and machine learning [4, 9, 13, 20, 22, 26,(2)where the optimal X∗ can be obtained by the singular value
thresholding operator deﬁned as:
Deﬁnition 1 (Singular value thresholding1 [1]). The
problem (2) has a closed form solution by the singular value∗ This work was done while the ﬁrst author was an intern at Microsoft
Research Asia.

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

1

thresholding (SVT) operator Sτ (·) as
,X∗ = Sτ (A) = UA Sτ (ΣA )VA(3)

where Sτ (x) = sgn(x) · max(|x| − τ, 0) is the soft shrinkis the SVD of A.
age operator [10], and UA ΣA VA
The major computational bottleneck for NNM and
WNNM problems is the necessity of solving Eq. (2) multiple times, where SVD computation occupies the largest
computation cost (e.g., O(mn min(m, n)) for a SVD [8].
In this paper, we propose a fast SVT technique to accelerate general NNM and WNNM methods. Our method ismotivated by the previous study of a randomized SVD proposed by Halko et al. [11], and we extend the original general method in several respects for better solving the NNM
and WNNM problems that we focus in this paper. As a result, we propose an algorithm that we call fast randomized
SVT (FRSVT). We present the connection between FRSVT
and low-rank approximation with both theoretical and empirical analyses, and show the effectiveness of the proposed
method via simulations and a few applications. Speciﬁcally,
this paper makes the following contributions:• We develop a successive truncated low-rank decomposition that can be generally applied to NNM and WNNM
problems. Our method achieves high speed and rarely
harms the convergence behavior of the host algorithms
without loss of accuracy.
• By exploiting the proximity of the range space over iterations, our method propagates the estimated basis of theprevious iteration to the next iteration for acceleration. We
call the proposed technique range propagation (RP).
• We provide theoretical analysis between the low-rank approximation and our FRSVT method. 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.

2. Related Works
Candès et al. [3] showed that, under some mild condition, the solution of NNM is equivalent to the solutionof rank minimization in conjunction with sparse outlier
model [3]. Inspired by the success of the convex surrogate for the rank minimization, low-rank observation is exploited by various computer vision applications, such as
rain removal [4], de-noising [9], inpainting [13], motion
segmentation by subspace clustering [20], structure frommotion [22], background subtraction [26], tag transduction[26], high dynamic range imaging [29], batch image alignment [31], photometric stereo [36], image rectiﬁcation [38],and nuclear norm regularized learning models [25].

Due to the high computational complexity of SVD in
the SVT operator, a fast SVT is always necessary for both
small- and large-scale problems for real-timeness and scalability, and there have been continued efforts to improve
the speed of low-rank algorithms. Liu et al. [21] efﬁciently
solve the NNM on a small matrix by factorizing a matrix
into two small matrices. Although the work achieved significant speed-up, global optimum cannot be guaranteed, because the factorization introduces non-convexity like otherfactorization methods [39, 6].
With retaining the advantage of convexity, Liu et al. [22]exactly solve RPCA on a small sub-sampled matrix andpropagate the seed solution to other parts via 1 ﬁlteringin a linear time. Since the work only focuses on RPCA, afast SVT method is still needed as a tool for NNM to be
applied to large-scale problems. Cai et al. [2] avoid explicit SVD computation using the dual of SVT. Since their
method uses Newton iterations with an inverse matrix, theinput matrix should be preprocessed by the complete orthogonal decomposition [8] (COD) to ensure a non-singularsquare matrix. The standard COD consists of twice of QR
decompositions with column/row pivoting, which requires
O(mn min(m, n)), so the reduction of computation complexity is still limited.
In other thread of works, it has been shown that the exact
SVD computation is unnecessary in the inner loop of NNM.
Mu et al. [26] proposed a compressed optimization by random projection. Ma et al. [25] solved NNM related problems with a linear-time approximate SVD [5]. However,these methods are occasionally unstable, because the input matrix itself is approximated by sampling or projection,
where the original information is impaired and the randomness leads to unstable incorrect results. In contrast to thesemethods, our method only approximates subspace bases toguarantee that most spectrum information is retained.

3. Fast Randomized Singular Value Thresholding (FRSVT)The basic idea of our method shares the idea of Liu etal. [22] and Ma et al. [25], in that the solution of Eq. (2)can be found by applying SVT to a small matrix instead of
the original large matrix. Instead of sampling columns orrows of a matrix as in [22, 25], our method extracts a small
core matrix by ﬁnding orthonormal bases with the unitaryinvariant property. Speciﬁcally, since the NNM deﬁned in
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),where Sτ (·) is the SVT operator.

(4)

Algorithm 1 Fast Randomized Singular Value Thresholding (FRSVT) algorithm
Input : A ∈ Rm×n , τ > 0, l = k + p > 0 and q ≥ 0. For
range propagation, the orthonormal column matrix Q̃ of the previous
iteration.
if not Range propagation thenSample Gaussian random matrix Ω ∈ Rn×l
Y = AΩ
Q = QR CP(Y)else
Sample Gaussian random matrix Ω ∈ Rn×p
Y = AΩ
QY = PartialOrthogonalization(Q̃, Y)
Q = [Q̃, QY ]end ifrepeatQ = QR(AA Q)
until η times[H, C] = QR(A Q)[W, P] = PolarDecomposition(C)[V, D] = EigenDecomposition(P)Sτ (A) = (QV) Sτ (D) (HWV)
Output : Sτ (A), Q̃ = QVProof. For the derivation, refer to the supplementary material. This also holds for most of the WSVT cases.
In general, SVT requires SVD computation, and its complexity is O(mn2 )2 . Based on Proposition 1, we can avoid
expensive computation by instead computing SVT on asmaller matrix B ∈ Rn×n , when Q ∈ Rm×n is available.
Given Q ∈ Rm×k that best approximates A by a rankk matrix, the complexity of SVD of B ∈ Rk×n becomesO(nk 2 ). Therefore, when k  n, the computation speed
can be signiﬁcantly 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, apartial (or truncated) SVD is frequently used to reduce thecomplexity in many prior arts. Our method similarly ﬁnds
a rank-k approximation (k < n) of the original matrix A as
A ≈ Âk = QB. It saves the computation of the ﬁrst step
as well as the second step, because the size of matrices Q
and B are reduced to m × k and k × n, respectively. By exploiting the observation that the major orthonormal k bases
evolve slowly over iterations, our method efﬁciently initializes matrix Q at each iteration by bypassing expensive random range estimation (Range propagation in Sec. 3.1). In
addition, by avoiding direct SVD computation, we can further reduce the computation as described in Sec. 3.2. Also,
we describe a target rank prediction technique for further
improving speed in Sec. 3.3.
2 Without loss of generality, we assume m ≥ n in this paper. In thecase m < n, the proposed method is analogous.

3.1. Finding Approximate Range
Inspired by Halko et al. [11], we ﬁrst estimate the orthonormal bases Q = [q1 , · · · , ql ] (where l ≥ k) such thatspan(Q) ⊆ Range(A) from a matrix compressed by random projection. Intuition of their randomized range ﬁnding
algorithm is as follows. By multiplying a random vector ωj ,a random 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 rank-k projection of A, ofwhich range is the target to capture, and E represents small
perturbation, then sample vector yj can be obtained by
yj = Ak ωj + Eωj .

(5)

Even though unwanted E may be included in {yj }, since
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, {yj } could
not span the entire Range(Ak ). 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 ], wherel = k + p, Q can be obtained by orthonormalizing Y.
When r = rank(A) < l, r bases are enough to span theentire Range(A). Indeed, when rank(A) < l, the range
ﬁnding algorithm is fairly accurate and close to the exact
method as we will see the theoretical analysis in Sec. 4.
Thus, we can reduce the dimension of Q for almost freeby estimating the rank and dominant bases by the QR decomposition with Column Pivoting (QR-CP) to Y. While
Halko et al. also proposed a complex algorithm to adaptively and approximately determine the number of bases
by different randomization for sampling, orthogonalization,and re-orthogonalization, our method is based on an exact method with the same complexity using QR-CP, andwe adaptively predict the sampling rate in a simpler manner (we will see in Sec. 3.3). Fortunately, in LAPACK, anefﬁcient QR-CP routine using level-3 BLAS (dgeqp3) isavailable. 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 andorgqr routines, respectively. After obtaining Q ∈ Rm×s ,where s = min(l, r), we can compute a small matrix B by
B = Q A ∈ Rs×n .
Range Propagation (RP) for Fast Range Finding
Based on the observation that Range(A(i) ) at the i-th iteration of NNM related problems is similar to the one atthe (i − 1)-th iteration, our method uses the singular vectors at the (i − 1)-th step as an initial approximation ofrange bases Q(i) at the i-th step. To capture the change ofthe range space more, additional p sample vectors {y} areNormalized magnitude

1

Image
Gaussian
Video

0.5

0
0
100
200
Singular value index

Figure 1. Illustration of singular value decaying. [Left] Gaussian
random, video and image samples. [Right] Decaying graphs ofsingular values. The Red, Green, Blue lines represent the graphsof a Gaussian random matrix, video, and image, respectively. The
spectrum of visual data decays signiﬁcantly fast.

newly sampled. We append {y} to the previous singular
vector matrix Q̃(i−1) as Q(i) = [Q̃(i−1) , y1 , · · · , yp ], andapply partial orthogonalization only for newly added {y}by the modiﬁed Gram-Schmidt procedure [8]. The numberof bases can subsequently be reduced by checking the rankwith QR-CP in the ﬁrst step of the power iteration.
Power Iteration Among the overall process in our algorithm, since the only approximation step is to estimate
the orthonormal column matrix Q only, the accuracy of
our algorithm depends only on this step. In Eq. (5), if
the magnitude of the action (i.e., spectrum) of Ak is not
dominant against E, the directions of sample vectors arebiased and may be included in Range(Ak )⊥ . This introduces accuracy degradation to the rest of the process. To
resolve this issue, Halko et al. [11] proposed the power iteration scheme, which makes the spectrum difference between Ak and E larger by estimating Q on (AA )η A. It
improves the chance of better capturing the range of Ak
from Y = (AA )η AΩ, while the singular vectors remainunchanged. Halko et al. also showed that η = 2 or 4 power
iterations are sufﬁcient for usual data of interest, and highly
accurate range ﬁnding can be achieved. As shown in Fig. 1,
decaying singular values of visual data is much faster thanGaussian random matrix. Our empirical tests also show thatη = 2 is enough and it is used in all our experiments.

3.2. 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 semi-deﬁnite matrices, SVD can be more efﬁciently
computed by Eigen decomposition (ED), which is generally
faster than SVD at least twice in our empirical tests. To apply ED to a general matrix, we form a positive semi-deﬁnite
matrix by the following decomposition:
Deﬁnition 2 (Polar decomposition [12]). Let X ∈ Cm×n ,m ≥ n. There exists a matrix W ∈ Cm×n and a unique
Hermitian positive semi-deﬁnite matrix P ∈ Cn×n suchthat
X = WP,W∗ W = I,

where I is the identity matrix. If rank(X) = n, then P is
positive deﬁnite 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. [12], which has a quadratic convergence
behavior. In our experiment, it converges at a small number of iterations (typically 7) with various different data,
which is consistent with the result of [2, 12]. Due to therequirement of the inverse operator in Newton iterations,
it is only applicable to non-singular square matrix. SinceB ∈ Rn×s is a full column rank matrix, the non-singularsquare matrix can be simply obtained from B = HC by
QR decomposition, where we call C ∈ Rs×s a core matrix
that is always non-singular and square. Contrary to Sec. 3.1,
no column pivoting is required.
We sequentially apply the polar decomposition and ED
on the core matrix to be 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 singular value matrix of B. Finally, Sτ (A) can beapproximated 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 either side ofrandom matrix multiplication). Overall algorithm is summarized in Algorithm 1.

3.3. Adaptive Rank Prediction (AP) Heuristic
For SVT, only singular vectors corresponding to the singular values that are greater than a certain threshold areneeded, and full SVD is unnecessary. Since the rank ofA(i) is unknown before SVD, predicting its rank can avoid
unnecessary computation. We observe that, in many NNM
related problems, the rank of A(i) is monotonically increasing or decreasing over iterations, and the rank is stabilized
as the number of iteration increases. As we shall see in thetheorem of error bound in Sec. 4, over-sampling is always
useful to reduce the expected error bound of FRSVT. Thus,
optimistically predicting rank allows to achieve both computational efﬁciency and stability.
The speed advantage of our method will be lost with ahigher sampling rate. In such a case, we resort to the truncated SVT by upper bounding the target rank. As shown innatural image statistics of Fig. 1, the rank of A(i) is generally stabilized at low-rank in many computer vision application. Usually, the ﬁnal accuracy is not harmed, as seen

in the successes of the truncated SVD in the NNM related
problems [17, 18, 19, 21].
Based on these observations, we deﬁne over-sampling
rate p as:
pi+1 =a,ρn,

if ri < li ,otherwise,(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), and · denotes the ceil operation. The prediction rule for sampling rate is deﬁned asli+1 = min(ri + pi+1 , b), where b = γn 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 singularvalues of A(i) that are larger than the threshold, i.e., theestimated rank of Sτ (A(i) ). Initially, we set l0 = 0.1b. In
the case of the range propagation, the number of columns inQ̃(i−1) becomes ri , and pi = li − ri .
When ri < li , Eq. (7) slightly over-samples, otherwise
it optimistically predicts the rank of the next iteration by
the a larger over-sampling rate. By virtue of low-ranknessof visual data shown in Fig. 1, the optimistic rule leads toaccurate estimate.

3.4. Computational Complexity
The computation of the sample matrix Y by random projection requires O(mnl), but we can reduce it to O(mn) by
the range propagation technique. The random matrix can be
multiplied to either left or right hand side of A. We multiply the random matrix to the longer side (left side, whenn < m) of A, so that the complexity of power iteration canbe further reduced. Under m ≥ n, the power scheme consists of O(mns) matrix multiplication and O(ns2 ) QR. The
polar decomposition and ED takes O(s3 ). Finally, the solution matrix composition takes O(mns). Most expensive
computations with O(mns) are simple matrix multiplications, so it can be easily parallelized.

4. Approximation BoundTo analyze the behavior of the proposed algorithm, it ismeaningful to see the approximation error bound. Except
for the randomized step in sampling Y, the other steps of
our algorithm is based on exact methods. Thus, the accuracy is only affected by the randomized range estimation, which computes a rank-k approximation of a matrix.
Since there is the only one approximation step, the proposed
method shares the same theorems with Halko et al. [11] as:
Theorem 1 (Average Frobenius3 error bounds by randomization [11]). Let Âk be a rank-k approximation ma3 For easy derivation, we provide the squared versions (the square of
the Frobenius norm) for Theorems 1 and 2. The original error bound (theFrobenius norm) can be obtained by taking the square root.

trix A ∈ Rm×n . For a target rank k ≥ 2 and an oversampling parameter p ≥ 2, where k + p ≤ min(m, n),draw an standard Gaussian matrix Ω ∈ Rn×(k+p) . Then,the theoretical minimum error and the upper bound satisfy
σj2 (A) ≤ EA−Âk 2F ≤ poly(v)·σj2 (A),j>kj>kwhere E denotes expectation with respect to the randommatrix, poly(v) is a function for v = {k, p} (without thepower iteration) or {k, p, η} (with the power iteration).
Halko et al. show that the upper bound of the error isclose to theoretical minimum error with high probability
in conjunction with some improving techniques (e.g., oversampling, power iteration), and the bound is pessimistic.
From Theorem 1, the following theorem is our main result for the bound of the approximate SVT.
Theorem 2 (Average error bound of the approximate
SVT). Let Sτ (·) be the SVT operator. Then, the average
error satisﬁes the following inequality.
ESτ (A)−Sτ (Âk )2F ≤ poly(v)·σj2 (A) −G(A),where G(A) =

j>kj>k min(σj (A), τ )

2

≥ 0.

Proof. The proof can be found in the appendix.
In Theorems 1 and 2, in the case without power iteration,poly(·) is deﬁned as poly(k, p) = (1 + k/(p − 1)). The
polynomial bound with power iteration with the Frobenius
norm representation has no simple form, so instead we canobserve the behavior of the error bound with power iteration
by referring to the spectral bound of Halko et al.
Consequently, our Theorem 2 asserts that the bound of
the approximate SVT is tighter than the error by rank-k approximation in Theorem 1 on average. Thus, we can seethat FRSVT can be safely used in practice. From Theorem
2, we see the following properties:• σk+1 (A) → and Sτ (A) − Sτ (Âk )2F may approach
to close to zero when k ≥ rank(A), which means it is almost exact.
• Fortunately, the bound becomes rapidly tighten as p or η
increase in both practice and theory, similar to Halko et al..
• The bound is independent of the size of the matrix.

5. Experimental Results
In this section, we ﬁrst evaluate the efﬁciency of the proposed method in comparison with other methods using simulation data. The evaluation is conducted by examining theperformance of a single SVT computation and also by assessing the performance of RPCA computation, which isarguably the most relevant NNM problem in computer vision today. We then show NNM applications in computervision using real-world data: subspace clustering, semionline weather artifact removal, and simultaneous multiimage alignment and rectiﬁcation. All the experiments are
conducted on a PC with Intel i7-3.4GHz and 16GB RAM.
The same shared parameters were used among algorithms.

5.1. Evaluation using Simulation Data
We quantitatively evaluate our method in comparison toother methods with synthetic matrices sampled from a standard Gaussian distribution. We evaluate the computationaltimes on Matlab 2010a and 2014a 64bits. Since the recent Matlab has been intensively optimized on Intel CPU(mainly due to improvement of Intel MKL), the computation efﬁciency has been noticeably improved. Therefore,
it is worth reporting performance difference on these two
versions of Matlab, because most related works have beenassessed with Matlab older than 2011a [2, 19, 21, 22, 26].
For a fair comparison, we turn off multi-threading functions
including maxNumCompThreads(1) in Matlab. We describe the implementation details of each method in the
supplementary material. Also, more complete testings andcomparisons can be found in the supplementary material.
Single SVT test Figure 2 compares speed and accuracy of SVT computation using SVT methods, such as
Matlab built-in SVD4 (baseline), Lanczos [16], FSVT [2],LTSVD [5] and our FRSVT (with / without RP). Except forthe baseline SVD, the others produce the truncated SVD ofthe input matrix, and it is used for SVT computation. For arank-k approximation, we compute rank-(k+p) approximations by setting the over-sampling rate p to 2 for Lanczos,
5 for LTSVD, 5 for FRSVT. While LTSVD in Fig. 2-(e)is faster than ours, the approximation error is signiﬁcantly
higher than other truncated SVD, and it results in slower
convergence in RPCA as we will see in the following.
Robust PCA test To see the convergence behavior of the
SVT methods for RPCA [3], in Table 1, we compare ourmethod to various SVT methods, such as SVD, LTSVD,BLWS [19], FSVT, RSVD [11], using an inexact augmented Lagrange multiplier method [17] (iALM, or called
alternating directional multiplier method)5 . We apply the
proposed adaptive rank prediction in Sec. 3.3 to LTSVD,RSVD and our FRSVT. LTSVD shows convergence degradation to achieve comparable accuracy with others due torough 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.
4 We apply either the economic size or full SVD according to the matrix
shape and report faster one, but we denote them as econSVD only.
5 We report the result of Mu et al. [26] in the supplementary material as
we could not reproduce their result.

Algorithms
iALMeconSVD10iALMeconSVD14
iALMLTSVD [5]iALMBLWS [19]iALMFSVT [2]iALMRSVD [11]iALMFRSVT (ours)
iALMFRSVT-RP (ours)

#NM
23
23
41
24
23
23
23
23

- #NM: The number of NNM,
- TT: Total elapsed time (sec),

SIT
947.9
46.5
1.7
2.2
110.0
3.6
1.5
1.3

TT
21668
1069
70
53
2527
81
33
30

SPG
–
20×
312×
405×
9×
268×
665×
716×

ERR
1.8e-71.8e-7
4.6e-7
4.8e-71.8e-71.8e-71.8e-71.8e-7- SIT: Elapsed time (sec) of a single iteration,- SPG: Speed-up gain against the baseline,- ERR: AGT − ÂF /AGT F .

Table 1. Quantitative Comparisons on RPCA. 4000 × 4000 matrices are used. The results of other sizes and convergence graphs areshown in the supplementary material. In iALM procedure, #NM
is the number of iteration to solve the NNM subproblem, whichcorresponds to the total number of iterations.
Type
A
B

Objective function
argmin L∗ + λS2,1
L,S

argmin ni=k+1 σi (L) + λS1

Cargmin

L,S

3

L,E,Γ i=1

αi L(i) ∗ + λE1Constraint
O = ZL + S
O=L+SO◦Γ=L+E

Table 2. Examples of NNM related objective functions. (a) Lowrank representation. (b) RPCA based on the non-convex truncated nuclear norm, where k is the target rank to be encouraged.
(c) Low-rank and sparse 3-order tensor decomposition with alignment. Here,  · 2,1 is l2,1 norm, {αi } are the balance parameters
among the unfolding matrices L(i) and αi = 1 is assumed. We
refer the basic tensor algebra notations denoted in [37].

5.2. More Applications
We show more general applications of our FRSVT
method to various types of low-rank optimization problems
summarized in Table 2, such as afﬁne constrained NNM(Type A), non-convex truncated NNM (Type B) and NNM
on tensor structure (Type C). We show them with typicalcomputer vision applications in the following.
Type A - Robust Subspace Clustering by Low-Rank
Representation (LRR) Many visual data is often characterized by a mixture of multiple subspaces. One of therecent 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, LRR has been shown robust
against large corruptions. The robust LRR can be formulated by Type A, which can be efﬁciently solved by iALM.
In this experiment, we use the same parameter settings andevaluation metrics suggested in [20].
We apply FRSVT to the robust LRR for motion segmentation on Hopkins155 [35] and for face recognition on a partof Extended Yale Database B [15] respectively. We use all
156 sequences in Hopkins155, and only the ﬁrst 10 classes
in Yale, 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

60
40

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

−2

10

−3

10

20
0
0

−4

500

1000

1500

Target rank(a) Computation time
Full rank 2000 × 2000
Matlab 2010a

2000

10

0

500

1000

1500

2000

80

−5

60
40

10

LTSVD
LANCZOSFRSVTFRSVT−RP

−10

10

20
0
0

−15

500

1000

1500

2000

10

0

Target rank

Target rank(b) Accuracy
Full rank 2000 × 2000

30

0

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP(c) Computation time
Low rank 2000 × 2000
Matlab 2010a

500

1000

1500

2000

Target rank(d) Accuracy
Low rank 2000 × 2000

25
20

0

10

SVD
LTSVD
LANCZOSFRSVTFRSVT−RP

Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP−1

10

Accuracy

80

100

Accuracy−1Computation Time (sec)

Theoretical
LTSVD
LANCZOSFRSVTFRSVT−RP10

Accuracy

Computation Time (sec)

120

10

100

Computation Time (sec)

0

120

15
10

−2

10

−3

10

5
0
0

−4

500

1000

1500

Target rank(e) Computation time
Full rank 2000 × 2000
Matlab 2014a

2000

10

0

500

1000

1500

2000

Target rank(f) Accuracy
Full rank 2000 × 2000

Figure 2. SVT comparisons across SVD methods. [Left] X-axis: Target rank, Y-axis: Elapsed time (Sec). [Right] X-axis: Target rank,
Y-axis: A∗ − Âk F /A∗ F , where A∗ = Sτ (DGT ), Âk = Sτ (Dk ), DGT is the input data, and Dk is approximated by each method.
For the low-rank matrix test, we generate input data matrices of which rank correspond to the target rank by multiplying two Gaussian
random matrices with m × r and r × n, while in other tests full rank matrices
are used. In (b,d,f), the theoretical minimum error boundby the rank truncation in Sτ (Dk ) is provided for guidance, and it is deﬁned as j>k σi (A∗ ), where σi (A∗ ) is computed by the built-in
SVD of Matlab with the assumption that the SVD is exact. For the low-rank case in (d), the theoretical minimum error is 0, so we omit thetheoretical bound. For other results with other sizes, refer to the supplementary material.

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.590

Table 3. Comparisons of the subspace segmentation algorithms on
Hopkins155 [35] (Motion) and Extended Yale Database B [15](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 for both LRR and LRR+Ours. These
results show the improvement of the computation time with retaining the same accuracy. Details on the computational time areshown in the supplementary material.

vectorizing each image. While Hopkins155 is only weakly
corrupted, Yale data is heavily corrupted by shadows, specularity and noise. As shown in Table 3, our method speeds
up LRR without degrading of the accuracy.
Type B - Semi-Online Weather Artifact Removal We
consider a weather artifacts (e.g. snow or rain) removal
problem in a video sequence. Since the weather artifacts
have non-deterministic appearance and sparse yet random
distribution in the spatio-temporal domain, we model theartifacts 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 windowmanner to leave moving objects in the latent images.
Unfortunately, since ﬁnding low-rank solution by the nuclear norm is based on the blessing of high dimensionality [3], it could be degenerated with small n frames as reported in Oh et al. [30]. With the assumption that the object motions are small during few frames, we can encourage
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]and l1 minimization based on iALM. We replace PSVT by

(a)

(b)

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

FRSVT with the partial thresholding. We transfer the previous basis estimation to the next n frame optimization, so
it can be regarded as semi-online algorithm.
Our algorithm produces n results simultaneously for nimages. With an n = 5 sliding window, the method basedon 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.52ms. The qualitative results ofour method can be found in Fig. 3, which shows plausible
snow removal effects.
Type C - Simultaneous Multi-Image Alignment and Rectiﬁcation Simultaneous multi-image alignment and rectiﬁcation problem is Type C optimization suggested byZhang et al. [37] (called SRALT). The method combines
the ideas of TILT [38] and RASL [37] on the 3rd -order tensor structure, which exploits low-rank texture property andnuclear norm-based misalignment error respectively. Withshowing the applicability of FRSVT to NNM on a 3rd -order

SVT1

SVT2

SVT3

Others

500

600

SRALT

Ours

0

(a)

(b)

(c)

Figure 4. Qualitative comparison with RASL [31] and SRALT [37]
+ Our FRSVT. (a) Samples of input images. (b) Results obtainedby RASL [31]. (c) Results obtained by SRALT [37] + Our FRSVT.
(d) The average images. (Top: for all input images, Middle: for
RASL, Bottom: for SRALT+Ours). Due to the angular bias shownin (a), RASL aligns the samples to be consistent at the biased angle, while the SRALT based method in (c) not only align the images, but also correct upright poses.

tensor, we present a new application of SRALT to gait data.
Since gait is a biometric signal spanned on not only spatial domain (2D), but also temporal domain, so it is natural
to represent the data by 3rd -order tensor [24]. In gait recognition, the exact alignment of acquired data is frequently assumed, but in real situation, the accuracy of the gait recognition could be varied according to the pivot angle of the camera and locations of moving persons of both training and test
data. Therefore, SRALT framework is useful to align many
gait data and to ﬁnd vertical angles as a preprocessing step.
We observe that the human silhouettes have the low-rank
texture property, and as observed in Lu et al. [24], cyclic
gait motions are spanned by a few number of bases.
We conduct the experiments on the Gait Challenge data
sets6 [34], whose size is O ∈ R32×22×3000 (3000 unaligned
images with the size 22 × 32). The set of the geometrictransformations Γ = {g1 , · · · , g3000 }, where gj ∈ Rp is a
p-group parameterization, and is parameterized by a 3-DoF
Euclidean transformation in our experiments. We randomly
apply translation and rotation to introduce misalignments.
Our method is implemented on the top of [37] by replacing
their SVT method with FRSVT. While the method of [37]is converged 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 reaches the even smaller objectivescore 206.46. The qualitative and the detailed computation
time comparisons are reported in Figs. 4 and 5, respectively.

6. Conclusions
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 our method can produce reliable approximation. The proposed method has been assessed using
the problems of afﬁne constrained NNM, non-convex NNM
6 The dataset is described in the supplementary material.

100

200

(d)

300

400

700

800

900

Computational Time (Seconds)

Figure 5. Computational time comparison on the registration andrectiﬁcation application [37]. The NNM for tensor consists of 3way SVTs. Computational times are measured for a single outer
iteration (almost 100 inner iterations).

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 low-rank applications without losing
accuracy and hurting convergence behavior. The major bottleneck of our method is in the power scheme, and we areinterested in further reducing the computational complexity
by effectively relaxing this computation block in the future.

Appendices – Proof of Theorem 2
The following propositions are useful for deriving the approximation bound.
Proposition 2 (Dual pseudo-contraction). Let Pτ (·) is a dual
operator of Sτ (·), such that X = Sτ (X) + Pτ (X). Then,
Sτ (A) − Sτ (B)2F ≤ A − B2F − Pτ (A) − Pτ (B)2F .
Proof. Since the derivation is straightforward based on thepseudo-contraction [32], we omit the details.
Proposition 3. Let Âk be a rank-k approximation of A.

2min Pτ (A) − Pτ (Âk )2F =j>k min(σj (A), τ ) .
Proof. This can be directly reached from Proposition 2.
Proof of Theorem 2.

By Propositions 2 and 3,

Sτ (A) − Sτ (Âk )2F
≤ A − Âk 2F − Pτ (A) − Pτ (Âk )2F

≤ A − Âk 2F −min(σj (A), τ )2 .
j>kTake the expectation to the both sides,
ESτ (A) − Sτ (Âk )2F

≤ EA − Âk 2F −min(σj (A), τ )2j>k≤ poly(v) ·σj2 (A) −min(σj (A), τ )2 ,
j>k

j>k

(by Theorem 1)σj2 (A) − G(A)= poly(v) ·
j>k



min(σj (A), τ )2
G(A) =j>k

Obviously, G(A) ≥ 0 by the deﬁnition.

*Acknowledgments We would like thank Zhouchen Lin for helpful comments, Xiaoqin Zhang and Yadong Mu for sharing their source
codes. The ﬁrst author was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No.
2010-0028680).

References[1] J.-F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on
Optimization, 2010. 1[2] J.-F. Cai and S. Osher. Fast singular value thresholding without singular value decomposition. Methods and Applicationsof Analysis, 2013. 2, 4, 6
[3] E. Candès, X. Li, Y. Ma, and J. Wright. Robust principal
component analysis? Journal of the ACM, 2011. 1, 2, 6, 7
[4] 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).
IEEE, 2013. 1, 2
[5] P. Drineas, R. Kannan, and M. W. Hahoney. Fast monte carloalgorithms for matrices II: Computing a low-rank approximation to a matrix. SIAM Journal of Computing, 2006. 2,
6[6] A. Eriksson and A. van den Hengel. Efﬁcient computation ofrobust weighted low-rank matrix approximations using the l1
norm. IEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI), 2012. 2[7] 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. 1[8] G. H. Golub and C. F. V. Loan. Matrix computations 3rd Ed.
Johns Hopkins University Press, 1996. 2, 4
[9] 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). IEEE, 2014. 1, 2
[10] E. T. Hale, W. Yin, and Y. Zhang. Fixed-point continuation
for l1 -minimization: Methodology and convergence. SIAM
Journal on Optimization, 2008. 2, 1
[11] N. Halko, P.-G. Martinsson, and J. A. Tropp. Finding
structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review,
2011. 2, 3, 4, 5, 6
[12] N. J. Higham. Computing the polar decomposition-with applications. SIAM Journal on Scientiﬁc and Statistical Computing, 1986. 4[13] 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. 1, 2
[14] S. Ji and J. Ye. An accelerated gradient method for trace
norm minimization. In ICML. ACM, 2009. 1[15] K.-C. Lee, J. Ho, and D. Kriegman. Acquiring linear subspaces for face recognition under variable lighting. IEEE
Transactions on Pattern Analysis and Machine Intelligence(TPAMI), 2005. 6, 7, 2

[16] Z. Lin. Some software packages for partial svd computation.
In arXiv:1108.1548, 2011. 6, 2
[17] Z. Lin, M. Chen, and Y. Ma. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices. Technical Report UILU-ENG-09-2215, UIUC, 2009.
1, 5, 6, 3
[18] Z. Lin, R. Liu, and Z. Su. Linearized alternating direction
method with adaptive penalty for low-rank representation. In
Advances in Neural Information Processing Systems (NIPS),
2011. 1, 5
[19] Z. Lin and S. Wei. A block lanczos with warm start techniquefor accelerating nuclear norm minimization algorithms. In
arXiv:1012.0365, 2010. 5, 6, 2
[20] G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma. Robust
recovery of subspace structures by low-rank representation.
IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 2013. 1, 2, 6
[21] G. Liu and S. Yan. Active subspace: Toward scalable lowrank learning. Neural Computation, 2012. 2, 5, 6
[22] R. Liu, Z. Lin, Z. Su, and J. Gao. Linear time principal component pursuit and its extensions using l1 ﬁltering. Neurocomputing, 2014. 1, 2, 6
[23] 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). IEEE,
2014. 1[24] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos. MPCA:
Multilinear principal component analysis of tensor objects.
IEEE Transactions on Neural Networks,, 2008. 8[25] S. Ma, D. Goldfarb, and L. Chen. Fixed point and bregman
iterative methods for matrix rank minimization. Mathematical Programming, 2011. 1, 2
[26] Y. Mu, J. Dong, X. Yuan, and S. Yan. Accelerated low-rankvisual recovery by random projection. In IEEE Conference
on Computer Vision and Pattern Recognition (CVPR). IEEE,
2011. 1, 2, 6, 3
[27] T.-H. Oh, H. Kim, Y.-W. Tai, J.-C. Bazin, and I. S. Kweon.
Partial sum minimization of singular values in rpca for lowlevel vision. In IEEE International Conference on Computer
Vision (ICCV). IEEE, 2013. 1, 6
[28] T.-H. Oh, J.-Y. Lee, and I. S. Kweon. High dynamic range
imaging by a rank-1 constraint. In IEEE International Conference on Image Processing (ICIP). IEEE, 2013. 1[29] T.-H. Oh, J.-Y. Lee, Y.-W. Tai, and I. S. Kweon. Robust
high dynamic range imaging by rank minimization. IEEE
Transactions on Pattern Analysis and Machine Intelligence(TPAMI), 2014. 1, 2[30] T.-H. Oh, Y.-W. Tai, J.-C. Bazin, H. Kim, and I. S.
Kweon. Partial sum minimization of singular values inrobust pca: Algorithm and applications. arXiv preprint
arXiv:1503.01444, 2015. 7[31] 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 PatternAnalysis and Machine Intelligence (TPAMI), 2012. 1, 2, 8, 5[32] G. Pierra. Decomposition through formalization in a productspace. Mathematical Programming, 1984. 8

[33] C. Sanderson. Armadillo: An open source c++ linear algebra library for fast prototyping and computationally intensive
experiments. Technical Report, NICTA, 2010. 2[34] S. Sarkar, P. J. Phillips, Z. Liu, I. R. Vega, P. Grother, andK. W. Bowyer. The humanid gait challenge problem: Data
sets, performance, and analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 2005. 8[35] 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). IEEE,
2007. 6, 7
[36] 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. 1, 2[37] X. Zhang, D. Wang, Z. Zhou, and Y. Ma. Simultaneous rectiﬁcation and alignment via robust recovery of low-rank tensors. In Advances in Neural Information Processing Systems(NIPS), 2013. 6, 7, 8, 5[38] Z. Zhang, A. Ganesh, X. Liang, and Y. Ma. TILT: Transform invariant low-rank textures. International Journal ofComputer Vision (IJCV), 2012. 1, 2, 7
[39] Y. Zheng, G. Liu, S. Sugimoto, S. Yan, and M. Okutomi.
Practical low-rank matrix approximation under robust l1 norm. In IEEE Conference on Computer Vision and Pattern
Recognition (CVPR). IEEE, 2012. 2