Nonconvex Relaxation Approaches to Robust Matrix Recovery

1 Introduction

In many computer vision and machine learning problems, a data matrix can be represented as a low-rank component plus a sparse component. For example, video surveillance is the superposition of low-rank background and sparse foreground; corrupted images can be approximated by low-rank images plus sparse data noises; collaborative filtering problems can be solved by recovering a low-rank rating matrix and a sparse noise matrix from a partially observed data matrix. Therefore, it is of great interest to recover the low-rank component and the sparse component from the observed data matrix. However, solving matrix recovery problems is not trivial.

An intuitive idea is to formulate matrix recovery as a minimization problem with ℓ0-norm (i.e., the number of nonzero entries) loss and matrix rank penalty. Unfortunately, this optimization problem is NP-hard. An important alternative, well-known as robust principal component analysis (RPCA) [Cand`es et al., 2011; Wright et al., 2009], relaxes this problem into minimizing matrix ℓ1-norm loss and nuclear-norm penalty. This alternative is tractable because the resulting optimization problem is convex. When the rank of the underlying low-rank matrix is sufficiently low, RPCA can exactly recover the low-rank and sparse matrices with high probability [Cand`es et al., 2011].

RPCA has become an effective tool for computer vision applications. Beyond decomposing a data matrix into a lowrank component and a sparse component, RPCA can also be reformulated to handle a variety of computer vision problems, and it achieved state of the art results. For example, Liu et al. extended RPCA to solve subspace clustering problems; Peng et al. reformulated RPCA to align images; Zhang et al. introduced a variant of RPCA to camera calibration problems; Cheng et al. resorted to RPCA for solving image segmentation problems; Wang and Zhang employed RPCA to solve colorization problems.

However, RPCA has two major limitations. First, the rank of the underlying low-rank component is not sufficiently low in many applications, which violates the assumptions of RPCA. For example, natural images cannot be effectively approximated by matrices with very low rank, so RPCA may not be an effective tool for image inpainting. Under such a circumstance RPCA has relatively low performance as verified by our experiments.

Second, RPCA does not have good performance in robust dimensionality reduction. As was mentioned in [Wright et al., 2009; Cand`es et al., 2011], RPCA seeks to overcome the brittleness of the conventional PCA when the data are grossly corrupted. Unfortunately, the performance of RPCA in dimensionality reduction has neither theoretical guarantee nor experiment verification. Unlike truncated singular value decomposition (SVD) which discards only the small singular values, RPCA shrinks all the singular values equally. As is discussed in the following paragraphs, such behavior of RPCA hurts the "good" singular values and thereby leads to biased results. The empirical results in Figure 1(b) 1(c) 1(d) verify this point.

Since RPCA is an important tool in computer vision and machine learning, fundamental improvements in RPCA methodology can significantly contribute to many real world applications. This work seeks to better recover the lowrank and sparse matrices via a more robust and less biased formulation. Our method is motivated by the nonconvex sparsity-inducing penalties in [Fan and Li, 2001; Zhang, 2010a]. Figure 1 demonstrates that our method fulfills matrix recovery and dimensionality reduction much more effectively than RPCA.

Figure 1: Results on matrix dimenionality deduction. 1(a) a color image with 50% entries added with Gaussian noises i.i.d. from N(0, 0.12 ); 1(b) 1(c) 1(d) are the low-rank components recovered from 1(a) by RPCA, each of which is of rank 200, 100, and 50, respectively. 1(e) 1(f) 1(g) are recovered by our method (NRMR).

Fan and Li pointed out that the ℓ1-norm penalty overpenalizes large entries of vectors. Moreover, they proposed three criteria for good penalty functions: unbiasedness, sparsity and continuity at the origin. The ℓ1-norm satisfies both sparsity and continuity, but it is biased. Based on these three properties, Fan and Li proposed a new penalty function called the smoothly clipped absolute deviation penalty (SCAD). Recently, Zhang proved that a so-called minmax concave plus (MCP) penalty also possesses the three properties and achieves better performance than SCAD. Both SCAD and MCP are nonconvex and nearly unbiased. Extensive experiments in [Breheny and Huang, 2011; Fan and Li, 2001; Zhang, 2010a; Shi et al., 2011; Zhang et al., 2012; Zhang and Tu, 2012] have demonstrated the superiority of SCAD and MCP over the ℓ1-norm penalty.

It is a well known fact that the matrix rank is the number of nonzero singular values of a matrix, and the nuclear norm is the sum of all the singular values of a matrix. We can thereby relate the matrix rank to the ℓ0-norm of a vector; in the same way we relate the matrix nuclear norm to the ℓ1-norm of a vector. On one hand, this relation implies that the nuclear norm over-penalizes large singular values in the same way that the ℓ1-norm over-penalizes large entries; in other words, the nuclear norm results in a biased estimator as well as the ℓ1-norm does. On the other hand, this relation encourages us to devise a new nonconvex loss function or penalty for dealing with bias in the rank-minimization problems [Mazumder et al., 2010].

In particular, we consider the natural extension of MCP functions on matrices. Moreover, we also propose and study a matrix γ-norm as a nonconvex relaxation of the matrix rank. Roughly speaking, the matrix γ-norm is the minimax concave (MCP) function of the matrix singular values. The γ-norm is characterized by a positive factor (say γ), and it is tighter to the matrix rank than the nuclear norm is. Moreover, its limit at γ →∞ is the nuclear norm.

We apply the matrix MCP function as a nonconvex relaxation of the ℓ0-norm loss and the matrix γ-norm as a nonconvex relaxation of the matrix rank penalty to handle robust matrix recovery. Our work offers the following contributions.

• We introduce the nonconvex sparsity-inducing penalty MCP function to the rank-minimization problem.
• We develop the Nonconvex relaxation based Robust Matrix Recovery (NRMR) model for the matrix recovery problem. NRMR alleviates over-penalization on "good" variables and thus better recovers the low-rank and sparse matrices.
• We devise a majorization-minimization augmented Lagrange multiplier algorithm (MM-ALM) to solve the NRMR model. We also propose a speedup strategy which solves NRMR as efficiently as the ALM algorithm solves RPCA.

The remainder of this paper is organized as follows. In Section 3 we formulate the matrix recovery problem. In Section 4 we apply MCP function to rank-minimization problems and obtain a called γ-norm with which we formulate the NRMR model. In Section 5 we devise an algorithm calledMM-ALM to solve the NRMR model. In Section 6 we empirically evaluate the performance of NRMR mainly in comparison with RPCA.

2 Notations

The notations in this paper are defined as follows. For a matrix A = [Aij] ∈ Rm×n , let ‖A‖0 be the ℓ0-norm (i.e., the number of nonzero entries of A), ‖A‖1 = ∑i,j |Aij| be the ℓ1-norm, ‖A‖F = (∑i,j A2ij )1/2 = (∑ri=1 σ2 i (A))1/2 be the Frobenius norm, and ‖A‖∗ = ∑ri=1 σi(A) be the nuclear norm, where σi(A) be the i-th largest singular value of A and r = min{m, n}. Finally, let I denote the identity matrix with appropriate size, 1n = [1, 1, · · · , 1]T ∈ Rn denote the all-one vector.

3 ProblemFormulation

Given an m×n matrix D, we are concerned with the problem of recovering a low-rank matrix L and a sparse matrix S such that L+S = D. It can be formulated as the following optimization problem:

Since the problem in (1) is NP-hard, it is intractable in practical applications. Based on the fact that ‖L‖∗ and ‖S‖1 are the best convex approximations of rank(L) and ‖S‖0, respectively, Cand`es et al. proposed the RPCA model, which is defined by

An efficient algorithm of solving RPCA is an augmented Lagrange multiplier algorithm (ALM), which was studied in [Lin et al., 2009]. Like other algorithms for a penalized nuclear norm minimization problem, the ALM algorithm is based on the singular value shrinkage operator [Cai et al., 2010].

4 Methodology

In Section 4.1 we introduce a nonconvex function that called the matrix MCP norm. In Section 4.2 we apply MCP to the rank minimization problem and obtained a nonconvex low-rank-inducing penalty called the γ-norm. With the matrix MCP norm and the γ-norm, in Section 4.3 we formulate the Nonconvex relaxation based Robust Matrix Recovery (NRMR) model. In Section 4.4 we discuss how to tune the parameters of NRMR.

4.1 The MatrixMCP Norm

Given a vector β = (β1, . . . , βp)T ∈ Rp , λ > 0, and γ > 1, the MCP function [Zhang, 2010a] is defined by

where

Here (z)+ = max{z, 0}. The matrix MCP norm is naturally defined by

Here and later, we denote ψγ(t) = ψ1,γ(t) and Mγ(A) = M1,γ(A) for notational simplicity. The function Mγ(·) is continuous w.r.t. γ. When γ → ∞, ψγ(t) → |t| and Mγ(·) becomes the ℓ1-norm; when γ → 1, it gives rise to a hard threshold operator corresponding to the ℓ0 norm. Thus, Mγ(·) bridges the ℓ1 and ℓ0 norm. The matrix MCP norm enjoys several properties as follows.

Proposition 1. The matrix MCP norm defined in (3) satisfies the following properties:

1. Mγ(A) ≥ 0, with equality iff A= 0;
2. Mγ(A) is concave w.r.t. |A|, where |A| = [|Aij|];
3. Mγ(A) is increasing in γ, Mγ(A) ≤ ‖A‖1, and limγ→∞ Mγ (A) = ‖A‖1.

It is worth pointing out that like the matrix rank and the ℓ0-norm, the matrix MCP norm and the following matrix γnorm are not real norms, because they are nonconvex and do not satisfy the triangle inequality of a norm.

4.2 The Matrix γ-Norm

Let σ(A) = (σ1(A), . . . , σr(A))T be a function from Rm×n to Rr+ where r = min{m, n}. Interestingly, we have that rank(A) = ‖σ(A)‖0, ‖A‖∗ = ‖σ(A)‖1, and ‖A‖F = ‖σ(A)‖2. This can help us better understand the connection between the vector ℓ1-norm and the matrix nuclear norm.

Fan and Li (2001) and Zhang (2010) proved that the ℓ1norm over-penalizes large components, leading to a biased estimator. Moreover, they showed that a nonconvex penalty such as SCAD and MCP yields an asymptotically unbiased estimator. This motivates us to apply the MCP function to the rank minimization problem to better approximate the matrix rank. We refer to this nonconvex function as the matrix γnorm and denote it by ‖A‖γ.

Particularly, the γ-norm of an m×n matrix A is defined by

Clearly, ‖A‖γ is nonconvex w.r.t. A. Moreover, it is easy to verify that it possesses the following properties.

1. ‖A‖γ ≥ 0, with equality iff A= 0;
2. ‖A‖γ is an increasing function of γ, ‖A‖γ ≤ ‖A‖∗, and limγ→+∞ = ‖A‖∗;
3. ‖A‖γ is unitarily invariant; that is, ‖UAV‖γ=‖A‖γ whenever U (m×m) and V (n×n) are orthonormal.

4.3 Formulation of NRMR

We use the matrix MCP norm as a surrogate of the ℓ0-norm and the γ-norm as a surrogate of the matrix rank, simultaneously. Accordingly, we have

where ‖L‖γ1 and Mγ2 (S) are characterized by parameters γ1 and γ2, respectively. Since our model is based on the idea of nonconvex relaxation, we call it Nonconvex relaxation based Robust Matrix Recovery (NRMR). From the properties of the matrix MCP norm and the γ-norm, it is clear that NRMR is a tighter approximation to Problem 1 than RPCA.

4.4 The Tuning Parameters

Finally, compared with RPCA which has only one tuning parameter λ, NRMR appears to be more complicated. It is therefore useful to discuss how to tune the parameters λ, γ1 and γ2. The parameter λ is the most significant among the three parameters and should be carefully chosen. Similar to RPCA, λ should be selected from the neighborhood of 1/√max{m, n}. Smaller λ gives rise to lower rank of L∗ . As for γ1 and γ2, experiments show that the results are insensitive to them in a large range, but γ1 and γ2 should be neither too large nor too small. When γ → +∞, NRMR simply becomes RPCA; when γ → 1, the local convergence results in poor solutions. Empirically speaking, γ1 and γ2 should be set to be small real numbers strictly greater than 1. We fix γ1 = γ2 = 4 in our experiments.

5 Algorithms

This section is organized as follows: in Sections 5.1 we discuss the generalized singular value shrinkage operator for the γ-norm; in Section 5.2 we present the MM-ALM algorithm for solving NRMR; in Section 5.3 we give a more efficient algorithm which is based on MM-ALM.

5.1 Generalized Singular Value Shrinkage Operator

The Singular Value Threshold (SVT) algorithm is proposed in [Cai et al., 2010] for solving the nuclear norm penalized problems. For τ ≥ 0, the singular value shrinkage operator Sτ is defined by

where X = UXΣXVTX is the singular value decomposition (SVD) of X. Similarly, we define the generalized singular value shrinkage operator Sτ,Λ and generalized shrinkage operator Dτ,W by

Here Λ and W are arbitrary elementwise nonnegative matrices. The optimality of Sτ,Λ and Dτ,W is shown in Theorem 3.

Theorem 3. For each τ ≥ 0, γ > 1, X,Y,Xold ∈ Rm×n , let Λ = (I - ΣXold /γ )+, W = (1m1Tn - |Xold |/γ )+, then Sτ,Λ and Dτ,W defined in (6) obey

where

is the locally linear approximation (LLA) of Mγ(A) given Aold .

5.2 The MM-ALMAlgorithm

The majorization-minimization (MM) method in [Hunter and Li, 2005] is a family of algorithms such as the locally quadratic approximation algorithm (LQA) in [Fan and Li, 2001], the locally linear approximation algorithm (LLA) in [Zou and Li, 2008; Zhang, 2010b], etc. Many optimization problems with nonconvex penalties can be solved by MM, e.g. those in [Fan and Li, 2001; Zou and Li, 2008; Zhang, 2010a; 2010b]. MM proceeds by iteratively solving a

Algorithm 1 NRMR via MM-ALM Algorithm

1: Input: Observation D ∈ Rm×n , parameters λ, γ1, γ2. 2: Initialize Y(0) , L∗0, S∗0; μ > 0; k = 0; 3: repeat 4: L(0) k+1 ←L∗k; S(0) k+1 ← S∗k; 5: Λ←Diag(1n - σ(L∗k)/γ1 )+; 6: W← (1m1Tn - |S∗k|/γ2 )+; 7: repeat 8: L(j+1) k+1 ←S1/μ,Λ(D-S(j) k+1 + Y(j) /μ); 9: S(j+1) k+1 ←Dλ/μ,W(D-L(j+1) k+1 + Y(j) /μ); 10: Y(j+1) ← Y(j) + μ(D - L(j+1) k+1 - S(j+1) k+1 ); 11: j ← j + 1; 12: until converged 13: L∗k+1 ←L(j+1) k+1 ; S∗k+1 ←S(j+1) k+1 ; k ← k + 1; 14: until converged 15: return L∗k+1, S∗k+1.

convex problem which locally approximates the original nonconvex problem in each step.

Our MM-ALM algorithm is a kind of the MM method. MM-ALM consists of an outer loop and an inner loop. In each iteration, the outer loop replaces the nonconvex problem by its LLA to form a weighted RPCA, while the inner loop is an ALM algorithm solving the resulting weighted RPCA problem. MM-ALM repeats these two loops until converged. MM-ALM is described in Algorithm 1.

LLA based outer loop

The outer loop is based on LLA. The basic idea is shown in the following proposition which indicates that the LLA majorizes the original problem.

Proposition 4. Let p(x) be a concave function on (0,∞). Then

Since the objective functionf(L, S) in Problem 4 is concave w.r.t. (σ(L), |S|), in each step we approximate ‖L‖γ1 + λMγ2 (S) by its LLA at (σ(Lold ), |Sold |) and obtain the following problem:

which majorizes the NRMR problem in (4). The problem above is indeed a weighted RPCA. Thus we can resort to the ALM algorithm in [Lin et al., 2009], which solves RPCA, to solve Problem 8.

The ALM algorithm based inner loop

We employ the ALM algorithm to solve Problem 8. The augmented Lagrangian function is

The ALM algorithm solves Problem 8 by alternately minimizing Lμ w.r.t. L and S, and maximizing w.r.t. Y; each step is guaranteed by Theorem 3. This constitutes the inner loop of MM-ALM.

Local Convergence

According to Proposition 4 and the properties of the ALM algorithm, we have that the objective function values of Problem 4 in each iteration of Algorithm 1 obey

Thus the objective function value is monotonically nonincreasing. If we denote the optimal low rank matrix by L∗ and sparse matrix by S∗ , then (L∗ , S∗ ) is a fixed point of Algorithm 1, that is, if we denote (Lk+1, Sk+1) = N(Lk, Sk), then we have (L∗ , S∗ ) = N(L∗ , S∗ ).

Though the convergence is local, local optimal solutions are still effective. The empirical results in [Breheny and Huang, 2011; Fan and Li, 2001; Shi et al., 2011; Zhang, 2010a; Zhang et al., 2012; Gong et al., 2012; Xiang et al., 2012] all show that the local optimal solutions to the nonconvex problems usually outperform global optimal solution to the LASSO [Tibshirani, 1996] problem. In Section 6 we will empirically demonstrate that the local optimal solutions to NRMR are superior over the global optimal solutions to RPCA.

5.3 Speedup Strategy for MM-ALM

As we see, each step of the MM procure results in a weighted RPCA; the MM-ALM algorithm repeatedly calls the ALM algorithm to solve the weighted RPCA problem in the inner loop. Thus the computation costs can be many times more than the ALM algorithm solving RPCA. To alleviate the computations, we propose an efficient strategy.

A feasible strategy is the one-step LLA studied in [Zou and Li, 2008]. The one-step LLA runs the outer loop only once instead of waiting to converge. In comparison, Algorithm 1 is actually the multi-stage LLA in [Zhang, 2010b]. Here we use the solutions to RPCA for initialization. Our off-line experiments reveal that the full MM-ALM leads to only marginal improvement over one-step LLA by sacrificing much more time. Experiments show that one-step LLA at most doubles the time of the ALM algorithm for solving RPCA.

6 Experiments

In this section, we conduct experiments mainly in comparison with RPCA, which is the state-of-the-art low-rank matrix recovery method. Though there are several nonconvex approaches such as log-det heuristic [Fazel et al., 2003], Matrix ALPS [Kyrillidis and Cevher, 2012], and SpaRCS [Waters et al., 2011], none of them were claimed to outperform RPCA in terms of accuracy, so we do not compare with these nonconvex methods.

In our experiments, RPCA is solved by the ALM algorithm in [Lin et al., 2009], and our NRMR is solved by the one-step LLA algorithm; we do not use the full MM-ALM algorithm because it is time consuming. We fix γ1 = γ2 = 4 for NRMR in all the following experiments. We use the relative square error (RSE) to evaluate matrix recovery accuracy. The RSE

is defined as follows:

where L∗ is the solution to RPCA or NRMR, and L is the ground truth. We use the number of singular value decompositions (SVD) to evaluate time efficiency, because the running time of all the aforementioned algorithms are dominated by the SVD in each iteration.

In Section 6.1 we conduct experiments on a set of synthetic data. In Section 6.2 we conduct experiments on natural image for the sake of presenting and comparing the results intuitively.

6.1 Synthetic Data

We first compare between RPCA and NRMR on a set of synthetic data. The data are generated in the following way. We generate L as a product of two m × r matrices whose entries are sampled i.i.d. from Gaussian distribution N (0, 1/m). The sparse component S is constructed by setting a proportion of entries to be ±1 and the rest to be zeros. Then we use RPCA and NRMR to recover the low-rank and sparse components from the data matrix D = L + S. We tune the parameters λ of both RPCA and NRMR and report the results where RSE is minimized. We report the RSE, number of SVD, and rank(L∗ ) in Table 1.

(b) noisy

(c) RPCA

(d) NRMR

(e) RPCA

(f) NRMR

Figure 3: Denoising results on images with Gaussian noises. 3(a) original images; 3(b) noisy images; 3(c) and 3(d) are obtained by tuning λ such that the RSEs are minimized; 3(e) and 3(f) are of the same rank (rank(L∗ ) = 100).

6.2 Natural Images

We use 50 images from the Berkeley Segmentation Dataset [Martin et al., 2001]. Each image consists of three components: Lred,Lgreen,Lblue ∈ Rm×n , each entry of which takes value in [0, 1). Our off-line experiments show that direct concatenation of the three components to a matrix leads to better results than segmenting the image into patches followed by concatenating the patches. So we directly stack the three components into an m×3n matrix L:

Then we add noises to L to generate contaminated data matrix D. We add two kinds of noises: i.i.d. Gaussian noises N(0, 0.22 ) to 50% pixels of the images, or salt and pepper noises to 20% pixels.

For both RPCA and NRMR, we tune the parameter λ that the recovered low-rank matrix has a rank of rank(L∗ ) = 100. We report the resulting RSE and average number of SVD in Figure 2. We also give in Figure 3 some visual comparisons among the recovered low-rank matrices, where λ is tuned such that RSE is minimized or that rank(L∗ ) = 100.

We can see from Table 1 and Figure 2 that the local optimal of NRMR always outperforms the global optimal of RPCA

(in terms of RSE), without adding much computation costs. The recovered images in Figure 3 all demonstrate that the results of our NRMR better approximate the original images.

7 Concluding Remarks

In this paper we have proposed a novel approach to the matrix recovery problem. In particular, we have defined robust matrix recovery as an optimization problem with a nonconvex loss function and a nonconvex penalty function. We have also devised a majorization-minimization augmented Lagrange multiplier algorithm for solving this nonconvex problem. Experiment results demonstrate that our nonconvex approach significantly outperforms RPCA-a classical convex approach-in terms of accuracy without much extra computation cost.

Our work has shed light on some important properties of the nonconvex low-rank-inducing penalty. The matrix γnorm studied in this paper is a tighter approximation to the matrix rank than the nuclear norm is, and it broadens our vision on matrix low-rank learning - from a convex relaxation to a more general nonconvex relaxation. Furthermore, it might be interesting to further explore the γ-norm as well as other potential nonconvex low-rank-inducing penalties to better solve the low-rank matrix learning problem.

Acknowledgments

This work has been supported in part by the Natural Science Foundations of China (No. 61070239) and the Scholarship Award for Excellent Doctoral Student granted by Ministry of Education.