# A Generalized Robust Minimization Framework for  Low-Rank Matrix Recovery

Authors: Wen-Ze Shao; Qi Ge; Zong-Liang Gan; Hai-Song Deng; Hai-Bo Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper considers the problem of recovering low-rank matrices which are heavily corrupted by outliers or large errors. To improve the robustness of existing recovery methods, the problem is solved by formulating it as a generalized nonsmooth nonconvex minimization functional via exploiting the Schattenp-norm (0<p≤1) and Lq(0<q≤1) seminorm. Two numerical algorithms are provided based on the augmented Lagrange multiplier (ALM) and accelerated proximal gradient (APG) methods as well as efficient root-finder strategies. Experimental results demonstrate that the proposed generalized approach is more inclusive and effective compared with state-of-the-art methods, either convex or nonconvex.


---

## Body

## 1. Introduction

In many practical applications, such as removing shadows and specularities from face images, separating foregrounds and backgrounds from monitored videos, ranking, and collaborative filtering, the observed data matrixC can naturally be decomposed into a low-rank matrix A and a corrupted matrix B. That is, C=A+B, where B can be arbitrarily large and is usually assumed to be sparse and unknown. The problem is whether it is possible to recover the low-rank matrix A from the observed C. Recently, it has been shown that the answer is affirmative as long as the corrupted matrix B is sufficiently sparse and the rank of the low-rank matrix A is sufficiently low [1, 2]. That is to say, under certain conditions one can exactly recover from C the low-rank matrix A with high probability by solving the following convex optimization problem, that is, the idealization of robust principal component analysis (RPCA):
(1)min⁡A,B∥A∥*+λ∥B∥1,s.t.C=A+B,
where ∥A∥* denotes the nuclear norm of the low-rank matrix A, that is, the sum of the singular values of A, ∥B∥1 denotes the L1 norm of the matrix B when seen as a vector, and λ is a positive tuning parameter. Recently, there has been a lot of research focusing on solving the RPCA problem, for example, [3–6].In spite of the great theoretical success and wide practical applications of RPCA (1), its major limitation should be claimed due to the use of nuclear and L1 norms as regularizers. Specifically, compared with the intrinsic rank constraint, that is, rank (A)<r0, the nuclear norm regularizer will not only do more harm to the large singular values of A but also lead to weaker shrinkage of the disturbed small singular values. It is not hard to make a similar analysis with the case of the L1 norm regularizer. Then, the performance of RPCA (1) in dimensionality reduction and outlier separation will not be as good as expected in some scenarios. In Section 3, it has been empirically demonstrated that RPCA (1) is not robust to the case that either the matrix A is not sufficiently low-rank or the matrix B is more grossly corrupted.To improve the robustness of RPCA [3–6], this paper proposes a generalized nonsmooth nonconvex minimization framework for low-rank matrix recovery by exploiting the Schatten p-norm (0<p≤1) and Lq(0<q≤1) seminorm. And two numerical algorithms are deduced based on the augmented Lagrange multiplier (ALM) and accelerated proximal gradient (APG) methods as well as efficient root-finder strategies. Experimental results show that the proposed generalized approach is more inclusive and effective compared with state-of-the-art methods [3–6]. Notice that much recently a nonconvex relaxation approach for low-rank matrix recovery [7] is proposed exploiting a nonconvex penalty called minmax concave plus and also a nonconvex loss function. However, our approach is different from [7] and is better in the terms of recovery accuracy and robustness than [7] as well as the other two nonconvex methods [8, 9].The paper is organized as follows. Section2 provides the generalized nonsmooth nonconvex minimization framework, including the problem formulation and two numerical algorithms based on the ALM and APG methods. Section 3 verifies the recovery performance of the proposed method and compares it against state-of-the-art methods. Finally, the paper is concluded in Section 4.


## 2. Proposed Model and Algorithms

### 2.1. Problem Formulation

Taking into account both the recovery robustness and computational efficiency, the Schattenp-norm is used to better approximate the intrinsic rank constraint rank(A)<r0; similarly, the Lq seminorm is exploited to replace the L1 norm of a matrix when seen as a vector. It is now intuitive to generalize RPCA as follows:
(2)min⁡A,B∥A∥Spp+λ∥B∥Lqq,s.t.C=A+B,
where 0<p,q≤1, ∥A∥Sp, and ∥B∥Lq are defined, respectively, in the following. Assume A,B,C∈Rm×n; then the Lq seminorm of a matrix B when seen as a vector can be defined as
(3)∥B∥Lq=(∑k=1,l=1m,n|Bk,l|q)1/q,
where Bk,l is the k,lth element of B; and the Schatten p-norm of a matrix A can be defined as
(4)∥A∥Sp=(∑j=1min⁡{m,n}σjp)1/p,
where σj is the jth singular value of A and the singular value decomposition (SVD) of A is A=UΣVT. Clearly, as p=q=1, (2) reduces to convex RPCA in (1); as 0<p,q<1, (2) corresponds to a constrained nonsmooth and nonconvex minimization problem. Now, it comes to study the numerical iteration schemes of (2).In recent several years, communities of signal processing and computational mathematics show more and more interests in developing efficient algorithms for nonlinear nonsmooth optimization problems [10], such as iterative soft thresholding, split Bregman iteration, accelerated proximal gradient, augmented Lagrange multiplier, and so on, which have significantly simplified sparse optimization problems including RPCA (1). In a similar spirit to [1, 3, 4, 7], this paper exploits the accelerated proximal gradient (APG) and augmented Lagrange multiplier (ALM) methods to solve the generalized minimization problem (2), considering the fact that APG and ALM are the two most popular numerical algorithms currently. As for ALM, it has been shown [11] that under certain general conditions, ALM converges Q linearly to the optimal solution. As for APG, though little is known about the actual convergence of its produced sequences, the O(1/k2) rate of convergence of the objective function that they achieve is optimal [10]. However, we should note that the above mentioned convergence results are not applicable to the new problem (2) because of its nonsmooth and nonconvex properties. In spite of that, empirical studies in Section 3 demonstrate that the two deduced algorithms in the following can both solve (2) very well, with empirically fast convergence rate.


### 2.2. ALM-Based Algorithm

This subsection exploits ALM to solve problem (2), which is a nonconvex extension of [3]. First of all, define functions f(A,B) and h(A,B) as
(5)f(A,B)=∥A∥Spp+λ∥B∥Lqq,h(A,B)=C-A-B.
According to ALM, the Lagrange function for (2) is given as
(6)L(A,B,X,μ)=f(A,B)+〈X,h(A,B)〉+μ2∥h(A,B)∥L22,
where X is a matrix of Lagrange multipliers and μ≥0 is the augmented Lagrange penalty parameter. It is seen that A,B can be solved iteratively by alternating minimization of L(A,B,X,μ). In the meanwhile, a continuation strategy is applied to μ≥0 in order to improve both the accuracy and efficiency of low-rank matrix recovery. Specifically, the iteration process is described in Algorithm 1.Algorithm 1:
Set1<ρ<2,λ,μ-,ε and K. Initialize k=0,B0,X0,μ0
while Residual error = ∥C-Ak+1-Bk+1∥L2max⁡{∥C∥∞/λ,∥C∥2}>ε
Ak+1=arg⁡minA⁡L(A,Bk,Xk,μk);
Bk+1=arg⁡minB⁡L(Ak+1,B,Xk,μk);
Xk+1=Xk+μk(C-Ak+1-Bk+1);
μk+1=min⁡{ρμk,μ-};
k=k+1;
if  k>K  break; end
endIn the following,Ak+1,Bk+1 are solved, respectively, by
(7)min⁡Aμk2∥(C-Bk+Xkμk)-A∥L22+∥A∥Spp,(8)min⁡Bμk2∥(C-Ak+1+Xkμk)-B∥L22+λ∥B∥Lqq.
As for (7), suppose that C-Bk+Xk/μk=ΦkΩk(Ψk)T; then (7) can be rewritten as
(9)min⁡Σ≥0μk2∥ΦkΩk(Ψk)T-UΣVT∥L22+trace(Σp).
According to von Neumann’s theorem on singular values [12], we have
(10)∥ΦkΩk(Ψk)T-UΣVT∥L22=trace((Ωk)TΩk)+trace(ΣTΣ)-2trace(ΦkΩk(Ψk)TUΣVT)≥trace((Ωk)TΩk)+trace(ΣTΣ)-2trace((Ψk)TΣ)=∥Ωk-Σ∥L22.
Hence, the minimization problem (7) can be approximated by
(11)min⁡Σ≥0μk2∥Ωk-Σ∥L22+trace(Σp).
More concretely, that is,
(12)min⁡σj≥0μk2(σj-ωjk)2+(σj)p.
Let B~k=C-Ak+1+Xk/μk. It is not hard either to find that the minimization problem (8) can be reduced to
(13)min⁡Bi,jμk2λ(Bi,j-B~i,jk)2+(Bi,j)q.
From (12) and (13), the kernel of Algorithm 1 is the following root-finder problem:
(14)min⁡ςβ2(ς-ξ)2+(ς)α.Here, we borrow the numerical idea in [13] to solve (14). For α=1,1/2,2/3, analytical solutions are used as calculated by Algorithms 2 and 3 in [13]. For all other α values, the numerical root-finder method Newton-Raphson is exploited to solve (14).


### 2.3. APG-Based Algorithm

This subsection exploits AGP as well as the continuation technique to solve problem (2), which is a nonconvex extension of [4]. First of all, a relaxed minimization problem is produced from (2); that is,
(15)min⁡A,BF(A,B)=min⁡A,Bνf(A,B)+g(A,B),
where g(A,B)=∥h(A,B)∥L22/2 and ν≥0 is a relaxation parameter. Obviously, F(A,B) is different from the Lagrange function of (2) in Section 2.2. However, instead of directly minimizing F(A,B), a sequence of separable quadratic approximations to F(A,B) is minimized, denoted as Q(A,B,Λ,Γ),
(16)Q(A,B,Λ,Γ,ν)=g(Λ,Γ)+〈∇g(Λ,Γ),(A,B)-(Λ,Γ)〉+Lg2∥(A,B)-(Λ,Γ)∥L22+νf(A,B),
where (Λ,Γ) are specifically chosen points and Lg=2. Then, A,B can be solved iteratively by alternating minimization of Q(A,B,Λ,Γ,ν) with reasonable choices of Λ,Γ. That is,
(17)min⁡A,BQ(A,B,Λ,Γ,ν)=min⁡A,Bνf(A,B)+Lg2∥(A,B)-(GΛ,GΓ)∥L22,
where
(18)GΛ=Λ-12∇g(Λ,Γ),GΓ=Γ-12∇g(Λ,Γ).
To assure both the accuracy and efficiency of minimizing Q(A,B,Λ,Γ,ν), two key strategies are taken into account deliberately. For one thing, (Λ,Γ) are determined by iterative smoothed computation as suggested in [14]. For another, the continuation technique is also applied to ν≥0, just similar to Algorithm 1. Moreover, the stopping criterion is identical to the one proposed in [15] and utilized in [4]. The iteration scheme is presented in Algorithm  2 specifically.Algorithm 2:
Set0<η<1,λ,ν-,τ,J. Initialize j=1,A0/-1,B0/-1,Λ0,Γ0,ν0.
while Residual error = ∥(Λj-Aj+1)+(Bj+1-Γj)∥L22max⁡(1,∥Λj∥L22+∥Γj∥L22)>τ
Λj=Aj+tj-1-1tj(Aj-Aj-1),Γj=Bj+tj-1-1tj(Bj-Bj-1);
GΛj=Λj-12(Λj+Γj-C),GΓj=Γj-12(Λj+Γj-C),
Aj+1=arg⁡minA⁡Q(A,Bj,Λj,Γj,νj);
Bj+1=arg⁡minB⁡Q(Aj+1,B,Λj,Γj,νj);
νj+1=max⁡{ηνj,ν-};
j=j+1;
ifj>J  break; end
endIn the following,Aj+1,Bj+1 can be solved, respectively, by
(19)min⁡ALg2∥A-GΛj∥L22+νj∥A∥Spp,(20)min⁡BLg2∥B-GΓj∥L22+νjλ∥B∥Lqq.
Similar to (12) and (13) in Section 2.2, both (19) and (20) are instances of rooter-find problems (14) and hence can be solved efficiently by borrowing the numerical idea in [13].


## 2.1. Problem Formulation

Taking into account both the recovery robustness and computational efficiency, the Schattenp-norm is used to better approximate the intrinsic rank constraint rank(A)<r0; similarly, the Lq seminorm is exploited to replace the L1 norm of a matrix when seen as a vector. It is now intuitive to generalize RPCA as follows:
(2)min⁡A,B∥A∥Spp+λ∥B∥Lqq,s.t.C=A+B,
where 0<p,q≤1, ∥A∥Sp, and ∥B∥Lq are defined, respectively, in the following. Assume A,B,C∈Rm×n; then the Lq seminorm of a matrix B when seen as a vector can be defined as
(3)∥B∥Lq=(∑k=1,l=1m,n|Bk,l|q)1/q,
where Bk,l is the k,lth element of B; and the Schatten p-norm of a matrix A can be defined as
(4)∥A∥Sp=(∑j=1min⁡{m,n}σjp)1/p,
where σj is the jth singular value of A and the singular value decomposition (SVD) of A is A=UΣVT. Clearly, as p=q=1, (2) reduces to convex RPCA in (1); as 0<p,q<1, (2) corresponds to a constrained nonsmooth and nonconvex minimization problem. Now, it comes to study the numerical iteration schemes of (2).In recent several years, communities of signal processing and computational mathematics show more and more interests in developing efficient algorithms for nonlinear nonsmooth optimization problems [10], such as iterative soft thresholding, split Bregman iteration, accelerated proximal gradient, augmented Lagrange multiplier, and so on, which have significantly simplified sparse optimization problems including RPCA (1). In a similar spirit to [1, 3, 4, 7], this paper exploits the accelerated proximal gradient (APG) and augmented Lagrange multiplier (ALM) methods to solve the generalized minimization problem (2), considering the fact that APG and ALM are the two most popular numerical algorithms currently. As for ALM, it has been shown [11] that under certain general conditions, ALM converges Q linearly to the optimal solution. As for APG, though little is known about the actual convergence of its produced sequences, the O(1/k2) rate of convergence of the objective function that they achieve is optimal [10]. However, we should note that the above mentioned convergence results are not applicable to the new problem (2) because of its nonsmooth and nonconvex properties. In spite of that, empirical studies in Section 3 demonstrate that the two deduced algorithms in the following can both solve (2) very well, with empirically fast convergence rate.


## 2.2. ALM-Based Algorithm

This subsection exploits ALM to solve problem (2), which is a nonconvex extension of [3]. First of all, define functions f(A,B) and h(A,B) as
(5)f(A,B)=∥A∥Spp+λ∥B∥Lqq,h(A,B)=C-A-B.
According to ALM, the Lagrange function for (2) is given as
(6)L(A,B,X,μ)=f(A,B)+〈X,h(A,B)〉+μ2∥h(A,B)∥L22,
where X is a matrix of Lagrange multipliers and μ≥0 is the augmented Lagrange penalty parameter. It is seen that A,B can be solved iteratively by alternating minimization of L(A,B,X,μ). In the meanwhile, a continuation strategy is applied to μ≥0 in order to improve both the accuracy and efficiency of low-rank matrix recovery. Specifically, the iteration process is described in Algorithm 1.Algorithm 1:
Set1<ρ<2,λ,μ-,ε and K. Initialize k=0,B0,X0,μ0
while Residual error = ∥C-Ak+1-Bk+1∥L2max⁡{∥C∥∞/λ,∥C∥2}>ε
Ak+1=arg⁡minA⁡L(A,Bk,Xk,μk);
Bk+1=arg⁡minB⁡L(Ak+1,B,Xk,μk);
Xk+1=Xk+μk(C-Ak+1-Bk+1);
μk+1=min⁡{ρμk,μ-};
k=k+1;
if  k>K  break; end
endIn the following,Ak+1,Bk+1 are solved, respectively, by
(7)min⁡Aμk2∥(C-Bk+Xkμk)-A∥L22+∥A∥Spp,(8)min⁡Bμk2∥(C-Ak+1+Xkμk)-B∥L22+λ∥B∥Lqq.
As for (7), suppose that C-Bk+Xk/μk=ΦkΩk(Ψk)T; then (7) can be rewritten as
(9)min⁡Σ≥0μk2∥ΦkΩk(Ψk)T-UΣVT∥L22+trace(Σp).
According to von Neumann’s theorem on singular values [12], we have
(10)∥ΦkΩk(Ψk)T-UΣVT∥L22=trace((Ωk)TΩk)+trace(ΣTΣ)-2trace(ΦkΩk(Ψk)TUΣVT)≥trace((Ωk)TΩk)+trace(ΣTΣ)-2trace((Ψk)TΣ)=∥Ωk-Σ∥L22.
Hence, the minimization problem (7) can be approximated by
(11)min⁡Σ≥0μk2∥Ωk-Σ∥L22+trace(Σp).
More concretely, that is,
(12)min⁡σj≥0μk2(σj-ωjk)2+(σj)p.
Let B~k=C-Ak+1+Xk/μk. It is not hard either to find that the minimization problem (8) can be reduced to
(13)min⁡Bi,jμk2λ(Bi,j-B~i,jk)2+(Bi,j)q.
From (12) and (13), the kernel of Algorithm 1 is the following root-finder problem:
(14)min⁡ςβ2(ς-ξ)2+(ς)α.Here, we borrow the numerical idea in [13] to solve (14). For α=1,1/2,2/3, analytical solutions are used as calculated by Algorithms 2 and 3 in [13]. For all other α values, the numerical root-finder method Newton-Raphson is exploited to solve (14).


## 2.3. APG-Based Algorithm

This subsection exploits AGP as well as the continuation technique to solve problem (2), which is a nonconvex extension of [4]. First of all, a relaxed minimization problem is produced from (2); that is,
(15)min⁡A,BF(A,B)=min⁡A,Bνf(A,B)+g(A,B),
where g(A,B)=∥h(A,B)∥L22/2 and ν≥0 is a relaxation parameter. Obviously, F(A,B) is different from the Lagrange function of (2) in Section 2.2. However, instead of directly minimizing F(A,B), a sequence of separable quadratic approximations to F(A,B) is minimized, denoted as Q(A,B,Λ,Γ),
(16)Q(A,B,Λ,Γ,ν)=g(Λ,Γ)+〈∇g(Λ,Γ),(A,B)-(Λ,Γ)〉+Lg2∥(A,B)-(Λ,Γ)∥L22+νf(A,B),
where (Λ,Γ) are specifically chosen points and Lg=2. Then, A,B can be solved iteratively by alternating minimization of Q(A,B,Λ,Γ,ν) with reasonable choices of Λ,Γ. That is,
(17)min⁡A,BQ(A,B,Λ,Γ,ν)=min⁡A,Bνf(A,B)+Lg2∥(A,B)-(GΛ,GΓ)∥L22,
where
(18)GΛ=Λ-12∇g(Λ,Γ),GΓ=Γ-12∇g(Λ,Γ).
To assure both the accuracy and efficiency of minimizing Q(A,B,Λ,Γ,ν), two key strategies are taken into account deliberately. For one thing, (Λ,Γ) are determined by iterative smoothed computation as suggested in [14]. For another, the continuation technique is also applied to ν≥0, just similar to Algorithm 1. Moreover, the stopping criterion is identical to the one proposed in [15] and utilized in [4]. The iteration scheme is presented in Algorithm  2 specifically.Algorithm 2:
Set0<η<1,λ,ν-,τ,J. Initialize j=1,A0/-1,B0/-1,Λ0,Γ0,ν0.
while Residual error = ∥(Λj-Aj+1)+(Bj+1-Γj)∥L22max⁡(1,∥Λj∥L22+∥Γj∥L22)>τ
Λj=Aj+tj-1-1tj(Aj-Aj-1),Γj=Bj+tj-1-1tj(Bj-Bj-1);
GΛj=Λj-12(Λj+Γj-C),GΓj=Γj-12(Λj+Γj-C),
Aj+1=arg⁡minA⁡Q(A,Bj,Λj,Γj,νj);
Bj+1=arg⁡minB⁡Q(Aj+1,B,Λj,Γj,νj);
νj+1=max⁡{ηνj,ν-};
j=j+1;
ifj>J  break; end
endIn the following,Aj+1,Bj+1 can be solved, respectively, by
(19)min⁡ALg2∥A-GΛj∥L22+νj∥A∥Spp,(20)min⁡BLg2∥B-GΓj∥L22+νjλ∥B∥Lqq.
Similar to (12) and (13) in Section 2.2, both (19) and (20) are instances of rooter-find problems (14) and hence can be solved efficiently by borrowing the numerical idea in [13].


## 3. Experimental Results

### 3.1. Experimental Settings

In this section, simulation experiments are designed and conducted to show the validity of our proposed approach. It firstly needs to produce available data usingC=A*+B*, in which A* and B* are, respectively, the true low-rank and sparse matrices that we wish to recover. Without loss of generality, A* is generated as a product of two m×r matrices whose entries are sampled i.i.d. from Gaussian distribution N(0,1/m), and the sparse matrix B* is constructed by setting a proportion of entries to be ±1 and the rest to be zeros. More specifically, if r and spr represent, respectively, the matrix rank and sparsity ratio, then the MATLAB v7.0 scripts for generating A* and B* can be given as(i)
A = 1/m*randn(m,r)*1/m*randn(r,m);(ii)
B = zeros(m,m);(iii)
p = randperm(m*m);(iv)
L = round(spr*m*m);(v)
B(p(1:L)) = sign(randn(L,1));In the following experiments, setm to 500, r to 50, 100, 150, and 200, and spr to 5%, 10%, 15%, and 20%. To be noted, the matrix recovery problem (2) roughly changes from easy to hard as r or spr changes from small to large. To assess the accuracy of low-matrix recovery, the relative squared error (RSE) is used, defined as
(21)RSE=∥A~-A*∥L2∥A*∥L2,
where A~ is the recovered low-rank matrix. And the number of SVD’s is used to evaluate the computational efficiency, since the running time of Algorithms 1 and 2 as well as [1, 3, 4, 7] is dominated by the SVD in each iteration. The experiments in this paper are conducted on a Lenovo computer equipped with an Intel Pentium (R) Core i5-3470 CPU (3.20 GHz) and 8 GB of RAM.


### 3.2. Comparison between Algorithm1 and RPCA [3]

In the literature, although several different numerical algorithms solving RPCA have been reported [3–6], the ALM method [3] is shown possessing the best performance in both accuracy and efficiency. Hence, this subsection compares Algorithm 1 with its convex and reduced version, that is, RPCA [3]. As implementing Algorithm 1, the parameters ρ,λ,ε,K are uniformly set as ρ=1.5,λ=1/500,ε=1e-7,K=100, and the parameter μ- is set as μ-=μ0·1e-7, where μ0 is set as 1.25/σ~ and σ~ is the largest singular value of C. Besides, B0,X0are both set as zero matrices. As for value choices of p and q, we set them as 0.85 based on empirical studies, despite the fact that it may produce more accurate recovery with choices adaptive to different r and spr.Experimental results of Algorithm1 and [3] are provided in Tables 1, 2, 3, and 4 corresponding to different settings. As the sparsity ratio spr equals 5%, it is obviously observed that Algorithm 1 performs perfectly in recovering the true rank of A* and is better than [3] in the term of RSE. It is also noticed that, as the sparsity ratio spr becomes larger, the recovery accuracy of both Algorithm 1 and [3] reduces, too. But it is still the case that Algorithm 1 behaves better than [3] no matter in the term of RSE or true rank recovery.Table 1
Low-rank recovery with Algorithm1 and [3] (spr = 5%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
6.21
e
-
6
4.68e − 6
20
30
48
50
100
8.31
e
-
6
5.99e − 6
30
30
100
100
150
2.00
e
-
3
3.49e − 5
30
100
270
150
200
7.03
e
-
3
3.85e − 5
30
100
306
200Table 2
Low-rank recovery with Algorithm1 and [3] (spr = 10%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
7.55e − 6
9.22
e
-
6
30
30
50
50
100
3.32
e
-
4
2.81e − 5
30
100
192
100
150
4.70
e
-
3
4.56e − 5
30
100
309
150
200
9.30
e
-
3
1.30e − 3
30
100
325
361Table 3
Low-rank recovery with Algorithm1 and [3] (spr = 15%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
1.75
e
-
5
1.52e − 5
30
30
50
50
100
2.50
e
-
3
3.87e − 5
30
100
287
100
150
7.90
e
-
3
6.01e − 5
30
100
328
151
200
1.23
e
-
2
3.60e − 3
30
100
334
398Table 4
Low-rank recovery with Algorithm1 and [3] (spr = 20%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
6.04
e
-
5
2.12e − 5
30
30
176
50
100
5.70
e
-
3
4.35e − 5
30
100
319
100
150
1.11
e
-
2
2.10e − 3
30
100
334
333
200
1.56
e
-
2
6.60e − 3
30
100
337
403One more point should be claimed which is that slightly lower RSE’s can be achieved by Algorithm1 as setting K=200. However, since the improvement in recovery accuracy is very limited, we just choose K=100 for computational efficiency.


### 3.3. Comparison between Algorithm2 and RPCA [4]

As running Algorithm2, the parameters η,λ,τ,J are uniformly set as η=0.9<1,λ=1/500,τ=1e-7, and J=200, and the parameter ν- is set as ν-=ν0·1e-9, where ν0 is set as the largest singular value of C, that is, σ~. In addition, A0/-1,B0/-1,Λ0, and Γ0 are all set as zero matrices. As for p and q, similar to the above manner, they are both set as 0.9 based on intensive empirical studies.Experimental results of Algorithm2 and [4] are provided in Table 5, 6, 7, and 8 corresponding to different settings. It is also remarkable that Algorithm 2 recovers the true rank of A* in almost all the scenarios, which is much superior to Algorithm 1. Its second advantage over Algorithm 1 is that it achieves slightly more robust recovery as B* is much grossly corrupted and A* is not sufficiently low-rank; for example, rank(A*) is 200. In spite of that, it is observed in other majority of cases that Algorithm 1 outperforms Algorithm 2 by O(1e-1) in the term of RSE. Therefore, it can be concluded that both algorithms possess their own advantages and disadvantages, and on the whole, Algorithm 1 shows better performance in terms of both recovery accuracy and efficiency.Table 5
Low-rank recovery with Algorithm2 and [4] (spr = 5%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.10
e
-
3
7.12e − 4
100
100
50
50
100
1.10
e
-
3
4.99e − 4
100
120
100
100
150
9.06
e
-
4
2.59e − 4
100
120
150
150
200
2.90
e
-
3
1.61e − 4
120
140
256
200Table 6
Low-rank recovery with Algorithm2 and [4] (spr = 10%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.40
e
-
3
6.35e − 4
100
120
50
50
100
1.40
e
-
3
4.29e − 4
100
120
100
100
150
2.50
e
-
3
2.47e − 4
120
120
241
150
200
7.80
e
-
3
2.91e − 4
120
140
290
200Table 7
Low-rank recovery with Algorithm2 and [4] (spr = 15%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.70
e
-
3
5.73e − 4
100
120
50
50
100
2.10
e
-
3
3.36e − 4
100
120
114
100
150
6.50
e
-
3
2.71e − 4
100
140
283
150
200
1.21
e
-
2
7.00e − 4
120
100
140
210Table 8
Low-rank recovery with Algorithm2 and [4] (spr = 20%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
2.00
e
-
3
5.28e − 4
100
120
70
50
100
4.10
e
-
3
3.31e − 4
120
120
259
100
150
1.11
e
-
2
4.36e − 4
120
140
297
150
200
1.62
e
-
2
2.10e − 3
120
140
305
282


### 3.4. Comparison between Proposed Approach and [7]

In the literature several nonconvex approaches for low-rank matrix recovery have also been proposed, for example, [7–9]. However, only [7] announces that it outperforms ALM-based RPCA [3] in the term of recovery accuracy.Table9 presents the RSE, number of SVD, and recovered rank achieved by [7] with sparsity ratios equal to 5% and 20%. Making comparison among Tables 1, 4, 5, 8, and 9, we can claim that both Algorithms 1 and 2 outperform [7] in terms of RSE and true rank recovery when A* is not sufficiently low-rank or B* is much grossly corrupted. In the meanwhile, we should also note that our method is computationally less efficient than [7] because of slightly more SVD’s used in each iteration, which is one of the future works to be studied.Table 9
Low-rank recovery with nonconvex approach [7].
rank(A*) r
RSE
# SVD
rank(A~)
5%
20%
5%
20%
5%
20%
50
8.29
e
-
7
2.91
e
-
6
40
60
50
50
100
9.97
e
-
7
1.67
e
-
1
68
70
100
195
150
7.39
e
-
5
3.94
e
-
1
72
73
419
259
200
1.19
e
-
1
5.36
e
-
1
73
73
471
266


### 3.5. Empirical Analysis on the Convergence of Algorithms1 and 2

As mentioned earlier, existed convergence results on ALM and APG in [10, 11] are not applicable to problem (2) due to the usage of nonconvex Lq seminorm and Schatten p-norm, which makes it difficult to conduct theoretical convergence analysis of the proposed algorithms, either. In spite of that, the empirical analysis can be made by plotting the residual error curve against iteration number for each algorithm. Specifically, the residual error curves are provided as the sparsity ratio equals 20% for both Algorithms 1 and 2, as shown, respectively, in Figures 1 and 2. It is obvious that the two deduced algorithms are of empirically fast convergence in each recovery scenario. Actually, this observation is also valid to other easier recovery cases with lower sparsity ratios. In addition, the number of iterations can be deduced from the residual error curves for each recovery problem.Residual error curves against iteration number for ALM-based Algorithm1 (spr = 20%); (a) rank(A*)=50, (b) rank(A*)=100, (c) rank(A*)=150, and (d) rank (A*)=200.
(a)
(b)
(c)
(d)Residual error curves against iteration number for APG-based Algorithm2 (spr = 20%); (a) rank(A*)=50, (b) rank(A*)=100, (c) rank(A*)=150, and (d) rank(A*)=200.
(a)
(b)
(c)
(d)


## 3.1. Experimental Settings

In this section, simulation experiments are designed and conducted to show the validity of our proposed approach. It firstly needs to produce available data usingC=A*+B*, in which A* and B* are, respectively, the true low-rank and sparse matrices that we wish to recover. Without loss of generality, A* is generated as a product of two m×r matrices whose entries are sampled i.i.d. from Gaussian distribution N(0,1/m), and the sparse matrix B* is constructed by setting a proportion of entries to be ±1 and the rest to be zeros. More specifically, if r and spr represent, respectively, the matrix rank and sparsity ratio, then the MATLAB v7.0 scripts for generating A* and B* can be given as(i)
A = 1/m*randn(m,r)*1/m*randn(r,m);(ii)
B = zeros(m,m);(iii)
p = randperm(m*m);(iv)
L = round(spr*m*m);(v)
B(p(1:L)) = sign(randn(L,1));In the following experiments, setm to 500, r to 50, 100, 150, and 200, and spr to 5%, 10%, 15%, and 20%. To be noted, the matrix recovery problem (2) roughly changes from easy to hard as r or spr changes from small to large. To assess the accuracy of low-matrix recovery, the relative squared error (RSE) is used, defined as
(21)RSE=∥A~-A*∥L2∥A*∥L2,
where A~ is the recovered low-rank matrix. And the number of SVD’s is used to evaluate the computational efficiency, since the running time of Algorithms 1 and 2 as well as [1, 3, 4, 7] is dominated by the SVD in each iteration. The experiments in this paper are conducted on a Lenovo computer equipped with an Intel Pentium (R) Core i5-3470 CPU (3.20 GHz) and 8 GB of RAM.


## 3.2. Comparison between Algorithm1 and RPCA [3]

In the literature, although several different numerical algorithms solving RPCA have been reported [3–6], the ALM method [3] is shown possessing the best performance in both accuracy and efficiency. Hence, this subsection compares Algorithm 1 with its convex and reduced version, that is, RPCA [3]. As implementing Algorithm 1, the parameters ρ,λ,ε,K are uniformly set as ρ=1.5,λ=1/500,ε=1e-7,K=100, and the parameter μ- is set as μ-=μ0·1e-7, where μ0 is set as 1.25/σ~ and σ~ is the largest singular value of C. Besides, B0,X0are both set as zero matrices. As for value choices of p and q, we set them as 0.85 based on empirical studies, despite the fact that it may produce more accurate recovery with choices adaptive to different r and spr.Experimental results of Algorithm1 and [3] are provided in Tables 1, 2, 3, and 4 corresponding to different settings. As the sparsity ratio spr equals 5%, it is obviously observed that Algorithm 1 performs perfectly in recovering the true rank of A* and is better than [3] in the term of RSE. It is also noticed that, as the sparsity ratio spr becomes larger, the recovery accuracy of both Algorithm 1 and [3] reduces, too. But it is still the case that Algorithm 1 behaves better than [3] no matter in the term of RSE or true rank recovery.Table 1
Low-rank recovery with Algorithm1 and [3] (spr = 5%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
6.21
e
-
6
4.68e − 6
20
30
48
50
100
8.31
e
-
6
5.99e − 6
30
30
100
100
150
2.00
e
-
3
3.49e − 5
30
100
270
150
200
7.03
e
-
3
3.85e − 5
30
100
306
200Table 2
Low-rank recovery with Algorithm1 and [3] (spr = 10%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
7.55e − 6
9.22
e
-
6
30
30
50
50
100
3.32
e
-
4
2.81e − 5
30
100
192
100
150
4.70
e
-
3
4.56e − 5
30
100
309
150
200
9.30
e
-
3
1.30e − 3
30
100
325
361Table 3
Low-rank recovery with Algorithm1 and [3] (spr = 15%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
1.75
e
-
5
1.52e − 5
30
30
50
50
100
2.50
e
-
3
3.87e − 5
30
100
287
100
150
7.90
e
-
3
6.01e − 5
30
100
328
151
200
1.23
e
-
2
3.60e − 3
30
100
334
398Table 4
Low-rank recovery with Algorithm1 and [3] (spr = 20%).
rank(A*) r
RSE
# SVD
Rank (A~)
[3]
Algorithm1
[3]
Algorithm1
[3]
Algorithm1
50
6.04
e
-
5
2.12e − 5
30
30
176
50
100
5.70
e
-
3
4.35e − 5
30
100
319
100
150
1.11
e
-
2
2.10e − 3
30
100
334
333
200
1.56
e
-
2
6.60e − 3
30
100
337
403One more point should be claimed which is that slightly lower RSE’s can be achieved by Algorithm1 as setting K=200. However, since the improvement in recovery accuracy is very limited, we just choose K=100 for computational efficiency.


## 3.3. Comparison between Algorithm2 and RPCA [4]

As running Algorithm2, the parameters η,λ,τ,J are uniformly set as η=0.9<1,λ=1/500,τ=1e-7, and J=200, and the parameter ν- is set as ν-=ν0·1e-9, where ν0 is set as the largest singular value of C, that is, σ~. In addition, A0/-1,B0/-1,Λ0, and Γ0 are all set as zero matrices. As for p and q, similar to the above manner, they are both set as 0.9 based on intensive empirical studies.Experimental results of Algorithm2 and [4] are provided in Table 5, 6, 7, and 8 corresponding to different settings. It is also remarkable that Algorithm 2 recovers the true rank of A* in almost all the scenarios, which is much superior to Algorithm 1. Its second advantage over Algorithm 1 is that it achieves slightly more robust recovery as B* is much grossly corrupted and A* is not sufficiently low-rank; for example, rank(A*) is 200. In spite of that, it is observed in other majority of cases that Algorithm 1 outperforms Algorithm 2 by O(1e-1) in the term of RSE. Therefore, it can be concluded that both algorithms possess their own advantages and disadvantages, and on the whole, Algorithm 1 shows better performance in terms of both recovery accuracy and efficiency.Table 5
Low-rank recovery with Algorithm2 and [4] (spr = 5%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.10
e
-
3
7.12e − 4
100
100
50
50
100
1.10
e
-
3
4.99e − 4
100
120
100
100
150
9.06
e
-
4
2.59e − 4
100
120
150
150
200
2.90
e
-
3
1.61e − 4
120
140
256
200Table 6
Low-rank recovery with Algorithm2 and [4] (spr = 10%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.40
e
-
3
6.35e − 4
100
120
50
50
100
1.40
e
-
3
4.29e − 4
100
120
100
100
150
2.50
e
-
3
2.47e − 4
120
120
241
150
200
7.80
e
-
3
2.91e − 4
120
140
290
200Table 7
Low-rank recovery with Algorithm2 and [4] (spr = 15%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
1.70
e
-
3
5.73e − 4
100
120
50
50
100
2.10
e
-
3
3.36e − 4
100
120
114
100
150
6.50
e
-
3
2.71e − 4
100
140
283
150
200
1.21
e
-
2
7.00e − 4
120
100
140
210Table 8
Low-rank recovery with Algorithm2 and [4] (spr = 20%).
rank(A*) r
RSE
# SVD
Rank (A~)
[4]
Algorithm2
[4]
Algorithm2
[4]
Algorithm2
50
2.00
e
-
3
5.28e − 4
100
120
70
50
100
4.10
e
-
3
3.31e − 4
120
120
259
100
150
1.11
e
-
2
4.36e − 4
120
140
297
150
200
1.62
e
-
2
2.10e − 3
120
140
305
282


## 3.4. Comparison between Proposed Approach and [7]

In the literature several nonconvex approaches for low-rank matrix recovery have also been proposed, for example, [7–9]. However, only [7] announces that it outperforms ALM-based RPCA [3] in the term of recovery accuracy.Table9 presents the RSE, number of SVD, and recovered rank achieved by [7] with sparsity ratios equal to 5% and 20%. Making comparison among Tables 1, 4, 5, 8, and 9, we can claim that both Algorithms 1 and 2 outperform [7] in terms of RSE and true rank recovery when A* is not sufficiently low-rank or B* is much grossly corrupted. In the meanwhile, we should also note that our method is computationally less efficient than [7] because of slightly more SVD’s used in each iteration, which is one of the future works to be studied.Table 9
Low-rank recovery with nonconvex approach [7].
rank(A*) r
RSE
# SVD
rank(A~)
5%
20%
5%
20%
5%
20%
50
8.29
e
-
7
2.91
e
-
6
40
60
50
50
100
9.97
e
-
7
1.67
e
-
1
68
70
100
195
150
7.39
e
-
5
3.94
e
-
1
72
73
419
259
200
1.19
e
-
1
5.36
e
-
1
73
73
471
266


## 3.5. Empirical Analysis on the Convergence of Algorithms1 and 2

As mentioned earlier, existed convergence results on ALM and APG in [10, 11] are not applicable to problem (2) due to the usage of nonconvex Lq seminorm and Schatten p-norm, which makes it difficult to conduct theoretical convergence analysis of the proposed algorithms, either. In spite of that, the empirical analysis can be made by plotting the residual error curve against iteration number for each algorithm. Specifically, the residual error curves are provided as the sparsity ratio equals 20% for both Algorithms 1 and 2, as shown, respectively, in Figures 1 and 2. It is obvious that the two deduced algorithms are of empirically fast convergence in each recovery scenario. Actually, this observation is also valid to other easier recovery cases with lower sparsity ratios. In addition, the number of iterations can be deduced from the residual error curves for each recovery problem.Residual error curves against iteration number for ALM-based Algorithm1 (spr = 20%); (a) rank(A*)=50, (b) rank(A*)=100, (c) rank(A*)=150, and (d) rank (A*)=200.
(a)
(b)
(c)
(d)Residual error curves against iteration number for APG-based Algorithm2 (spr = 20%); (a) rank(A*)=50, (b) rank(A*)=100, (c) rank(A*)=150, and (d) rank(A*)=200.
(a)
(b)
(c)
(d)


## 4. Conclusions and Discussions

In this paper, a generalized robust minimization framework is proposed for low-rank matrix recovery by exploiting the Schattenp-norm (0<p≤1) and Lq(0<q≤1) seminorm. And two numerical algorithms are deduced based on the ALM and APG methods as well as efficient root-finder techniques. Experimental results demonstrate that the proposed algorithms possess their own advantages and disadvantages and both perform more effectively than state-of-the-art methods, either convex or nonconvex, in terms of both RSE and true rank recovery.Note that this paper does not consider the influence of additive noise on the proposed algorithms, which actually corresponds to the problem of noisy RPCA [16, 17]. As claimed in [17], noisy RPCA is intrinsically different from the RPCA problem, that is, the focus of this paper. Indeed, the proposed algorithms in this paper are not quite robust to additive noise, just the same as many existing approaches to RPCA, for example, [1–4, 6–9]. To some degree, this observation coincides with the investigations in [18, 19]; that is, the Lq seminorm as a sparsity-enforcing penalty is vulnerable against the influence of additive noise on the data, as it resembles the L0 seminorm when q approaches 0, in spite of the fact that q in Algorithms 1 and 2 is chosen, respectively, as 0.85 and 0.9. Our future research topic is to extend the proposed algorithms to the noisy RPCA problem with applications to the field of image and vision computing.



---

*Source: 656074-2014-05-06.xml*