advertisement

RoseHulman Undergraduate Mathematics Journal Conditions for Robust Principal Component Analysis Michael Hornsteina Volume 12, No. 2, Fall 2011 Sponsored by Rose-Hulman Institute of Technology Department of Mathematics Terre Haute, IN 47803 Email: mathjournal@rose-hulman.edu http://www.rose-hulman.edu/mathjournal a Stanford University Rose-Hulman Undergraduate Mathematics Journal Volume 12, No. 2, Fall 2011 Conditions for Robust Principal Component Analysis Michael Hornstein Abstract. Principal Component Analysis (PCA) is the problem of finding a lowrank approximation to a matrix. It is a central problem in statistics, but it is sensitive to sparse errors with large magnitudes. Robust PCA addresses this problem by decomposing a matrix into the sum of a low-rank matrix and a sparse matrix, thereby separating out the sparse errors. This paper provides a background in robust PCA and investigates the conditions under which an optimization problem, Principal Component Pursuit (PCP), solves the robust PCA problem. Before introducing robust PCA, we discuss a related problem, sparse signal recovery (SSR), the problem of finding the sparsest solution to an underdetermined system of linear equations. The concepts used to solve SSR are analogous to the concepts used to solve robust PCA, so presenting the SSR problem gives insight into robust PCA. After analyzing robust PCA, we present the results of numerical experiments that test whether PCP can solve the robust PCA problem even if previously proven sufficient conditions are violated. Acknowledgements: I would like to thank my advisor Prof. Deanna Needell for introducing me to robust PCA and helping me through the process of writing the paper. Her feedback greatly improved the manuscript. I would like to thank the reviewer for helpful suggestions which improved both the proofs and the readability of the paper. This paper was written for the Vertical Integration of Research and Education program (VIGRE) in the Stanford Statistics Department. Page 138 1 RHIT Undergrad. Math. J., Vol. 12, No. 2 Introduction Principal Component Analysis (PCA) is a central problem in statistics which finds a low-rank approximation to a matrix [7, 9, 10]. If each column of a matrix M represents a variable, and each entry within a column represents a measurement of the variable, then PCA finds a smaller set of variables that retain much of the information present in the original variables. PCA provides insight into data by identifying the coordinate directions along which the data vary the most. By approximating a dense matrix with a low-rank matrix, PCA also generates a compressed representation of data. PCA is defined formally as the optimization problem of finding the best rank-r approximation to an n × n matrix. The approximation error is measured using the induced `2 norm. For an n × n matrix M , and a fixed rank r less than or equal to n, PCA can be expressed as the following optimization problem: minimize subject to kM − Lk2 rank(L) ≤ r. PCA is sensitive to errors added to the matrix M . In particular, PCA is highly sensitive to sparse errors with large magnitudes [2]. Corrupting a single entry of the matrix can markedly change the solution to the PCA optimization problem. This phenomenon can magnify the effects of data corruption, which often occurs in practice through round-off error, malicious tampering, or other sources. Researchers have searched for methods to make PCA less sensitive to errors. To make PCA less sensitive to errors, the PCA optimization problem can be replaced with a different optimization problem that separates out the errors. Robust PCA is the problem of performing low-rank approximation in a way that is insensitive to sparse errors with large magnitudes. More concretely, robust PCA is the problem of separating a matrix into the sum of a low-rank matrix and a sparse matrix of errors. Given a low-rank matrix L0 and a sparse matrix of errors S0 , the goal of robust PCA is to recover L0 and S0 from their sum M = L0 + S0 . Performing this separation prevents the errors from obscuring the low-rank component. Like PCA, which has applications in multiple fields, robust PCA also has numerous applications. Applications of robust PCA are surveyed in [2]. Robust PCA can be used to analyze videos by separating the background of a video from moving objects. If each column of a matrix encodes a video frame, then a low-rank approximation represents the background and the sparse “errors” represent moving objects. Another example made famous recently by NetFlix is collaborative filtering, in which entry (i, j) of a matrix represents the preference of user i for category j. A low-rank approximation to the matrix can be used to predict users’ preferences for categories that they have not explicitly ranked. Robust PCA can be used to separate out errors that arise in the matrix due to data tampering or people with atypical preferences. To solve the robust PCA problem, Candès et al. introduce an optimization problem called Principal Component Pursuit (PCP) [2]. This paper investigates the conditions RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 139 under which PCP successfully recovers a low-rank matrix L0 and a sparse matrix S0 from the sum M = L0 + S0 . To understand robust PCA and its solution via PCP, it is helpful to first consider a related problem, sparse signal recovery (SSR). SSR is the problem of identifying the sparsest vector in the set of solutions to an underdetermined system of linear equations. SSR can be solved by replacing an intuitively clear but intractable optimization problem with a tractable convex optimization problem. The concepts used to solve SSR are analogous to the concepts used to solve robust PCA. We will discuss SSR before robust PCA to establish the background for robust PCA. The paper is organized as follows. Section 2 introduces the sparse signal recovery problem. Section 3 discusses robust PCA by analogy with sparse signal recovery. Section 4 investigates the conditions under which Principal Component Pursuit solves the robust PCA problem and provides the results of numerical experiments that test whether the conditions presented in [2] are necessary. 2 2.1 Sparse Signal Recovery Overview The sparse signal recovery problem is the problem of identifying the sparsest solution (i.e. the one with fewest nonzero entries) to an underdetermined system of equations Ax = y, where x and y are vectors and A is a wide matrix. Because A has fewer rows than columns, the system of equations has multiple solutions. Recovering the sparsest solution is only possible if there is a unique sparsest vector in the solution set to Ax = y. We can also understand the SSR problem by thinking of A as a linear transformation. In an SSR problem instance, the matrix A maps the vector x0 to the vector y. We observe y, but not x0 . Because A has a nonzero null space, A maps multiple vectors to y. The SSR problem is to recover x0 from among the set of vectors that A maps to y. If A maps multiple equally sparse vectors to y, then we have no way to distinguish x0 from those other vectors, so it is impossible to recover x0 . Hence we need a condition on A to ensure that x0 is the unique sparsest vector in the set of solutions to Ax = y. We call a vector S-sparse if it has S or fewer nonzero entries: kxk0 = #{i : xi 6= 0} ≤ S, where the `0 quasi-norm kxk0 denotes the number of nonzero entries of x. For convenience, we refer to this quasi-norm as a norm. The condition on A is that we require A to be one-to-one on S-sparse vectors. If A is one-to-one on S-sparse vectors, and x0 is an S-sparse vector such that Ax0 = y, then A does not map any other S-sparse vectors to y. Hence x0 is the unique sparsest solution to Ax = y, so the SSR problem has a unique solution. Later we discuss a condition that can be imposed on A to ensure that A is one-to-one on S-sparse vectors. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 140 2.2 Solving SSR by Formulating an Optimization Problem SSR can be solved by formulating a tractable optimization problem. To explain the optimization problem, we first discuss a simpler optimization problem using the `0 norm. The optimization problem that directly corresponds to the SSR problem is minkxk0 subject to Ax = y. This optimization problem has no efficient solution, and is NP-hard in general [11]. We adopt the strategy of replacing the `0 objective function with a related objective function that yields the same solution (under appropriate conditions), but that can be minimized efficiently. Convex functions can be minimized efficiently using convex optimization software [1], so we replace the `0 norm with a convex function. The convex norm that we choose comes from the family of norms called `p norms, defined as follows: kxkp := n X ! p1 |xi |p . i=1 Examples of the `p norm are the `2 norm, v u n uX |xi |2 , kxk2 = t i=1 which is Euclidean length, and the `1 norm, kxk1 = n X |xi |, i=1 which is the sum of the absolute values of the components of x. To choose an `p norm for the SSR objective function, we examine the unit balls of the `p norms. The unit ball for a norm is the set of points of norm 1. The unit ball for the `0 norm is the set of points along the coordinate axes. We are looking for a convex norm that approximates the `0 norm. The unit balls for the `p norms are shown in Figure 1, with p decreasing from the outermost ball to the innermost. As p decreases, the `p unit balls shrink down to hug the coordinate axes more closely, thus better approximating the `0 unit ball. As shown in Figure 1, for p ≥ 1, the `p norm is convex, but for p < 1, the `p quasi-norm is not convex, so p = 1 is the smallest value of p for which the `p norm is convex. Replacing the `0 norm with the `1 norm, the resulting optimization problem is minkxk1 subject to Ax = y. An advantage of the `1 norm is that the minimization problem can be solved using linear programming [1]. Under appropriate conditions, the solution to this optimization problem RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 141 Figure 1: This figure shows `p balls for several values of p. For p ≥ 1, the `p norm is convex, and for 0 ≤ p < 1, the `p quasi-norm is not convex. The `0 unit ball is the coordinate axes. Page 142 RHIT Undergrad. Math. J., Vol. 12, No. 2 is the same as the solution to the `0 optimization problem [3, 4, 5, 6]. If A is one-to-one on S-sparse vectors, then the SSR problem is well-defined (that is, it has at most one solution). If A satisfies a stronger condition, called the restricted isometry property, then in addition to being well-defined, `1 minimization has the same solution as `0 minimization, so SSR can be solved efficiently [3, 5]. 2.3 Conditions Under Which the SSR Problem Is Well-Defined We have seen that if A is one-to-one on S-sparse vectors, then there is at most one S-sparse solution to Ax = y, so the SSR problem is well-defined. Although this condition clarifies when the SSR problem is well-defined, it is not intuitively obvious which matrices are oneto-one on S-sparse vectors. We now discuss a necessary and sufficient condition for a matrix A to be one-to-one on S-sparse vectors. The condition helps us to identify matrices that are one-to-one on S-sparse vectors and gives some geometric intuition regarding what it means for a matrix to be one-to-one on S-sparse vectors. Proposition 2.1. The matrix A is one-to-one on S-sparse vectors if and only if every subset of 2S columns of A is linearly independent. Proof. Assume that every subset of 2S columns of A is linearly independent. We show that A is one-to-one on S-sparse vectors. Let x1 and x2 be S-sparse vectors, and assume that Ax1 = Ax2 . Then A(x1 − x2 ) = 0. Because x1 and x2 are S-sparse vectors, x1 − x2 is a 2S-sparse vector, so A(x1 − x2 ) is a linear combination of at most 2S columns of A. Because every subset of 2S columns of A is linearly independent, it must be that x1 − x2 = 0, so x1 = x2 . Thus A is one-to-one on S-sparse vectors. Now assume that A is one-to-one on S-sparse vectors. We show that every subset of 2S columns of A is linearly independent. Let B be a subset of 2S columns of A. Because A is one-to-one on S-sparse vectors, and B is a subset of the columns of A, B must be one-to-one on S-sparse vectors. Let x be a nonzero 2S-dimensional vector. The vector x can be written as x = x1 − x2 , where x1 and x2 are distinct S-sparse vectors. It follows that Bx = Bx1 − Bx2 . If Bx = 0, then Bx1 = Bx2 , which contradicts the assumption that B is one-to-one on S-sparse vectors. Therefore for any nonzero 2S-dimensional vector x, Bx is nonzero, so the columns of B are linearly independent. We can also characterize matrices A that are one-to-one on S-sparse vectors in terms of their nullspaces. Corollary 2.2. The matrix A is one-to-one on S-sparse vectors if and only if the nullspace of A does not contain a nonzero 2S-sparse vector. Proof. A nonzero 2S-sparse vector x in the nullspace of A corresponds to a linearly dependent subset of at most 2S columns of A. This can be seen by interpreting the components of x as the coefficients of a linear combination of columns of A. By Proposition 2.1, it follows that A is one-to-one on S-sparse vectors if and only if there is no 2S-sparse vector in the nullspace of A. RHIT Undergrad. Math. J., Vol. 12, No. 2 2.4 Page 143 Conditions Under Which SSR Can Be Solved Efficiently The condition that A is one-to-one on S-sparse vectors guarantees that Ax = y has at most one S-sparse solution but does not guarantee that an S-sparse solution can be found efficiently. If A satisfies a stronger condition, called the restricted isometry property (RIP), then the sparsest solution also achieves the minimum `1 norm, so SSR can be solved efficiently by `1 minimization. We present the restricted isometry property [5], which is defined in terms of the S-restricted isometry property and the isometry constant δS . An “isometry” is a mapping that preserves length, so the restricted isometry property is connected with orthonormal matrices, which preserve the length of vectors that they act on. The S-restricted isometry property requires that all submatrices of A consisting of S columns approximately act as isometries. That is, all submatrices of S columns of A approximately preserve the length of vectors to which they are applied. The S-restricted isometry constant, δS , is a property of A that measures the extent to which submatrices of S columns of A deviate from being isometries. If the value of δS is sufficiently small, then subsets of S columns of A are approximately orthonormal, so A satisfies the S-restricted isometry property. We now formally define the δS -restricted isometry constant. If T is a set of column indices, we define AT to be the submatrix of A consisting of the columns specified by T . The δS -restricted isometry constant [5] is the smallest quantity that satisfies the inequality (1 − δS )kck22 ≤ kAT ck22 ≤ (1 + δS )kck22 , for all submatrices AT of A with S columns and all vectors c in the domain of AT . Dividing through by kck2 , we get the inequality 1 − δS ≤ kAT ck22 ≤ 1 + δS . kck22 This form of the inequality reveals that if the S-restricted isometry constant of A is δS , then each submatrix consisting of S√columns of A rescales √ the magnitude of each vector in its domain by a factor of at least 1 − δS and at most 1 + δS . We now prove that the 2S-restricted isometry property implies that A is one-to-one on S-sparse vectors. Proposition 2.3. If δ2S < 1, then A is one-to-one on S-sparse vectors. Proof. Let AT be an arbitrary submatrix of A with 2S columns. If δ2S < 1, then 1 − δ2S > 0, so kAT ck22 0 < 1 − δ2S ≤ kck22 for all nonzero vectors c in the domain of AT . By the above inequality, it follows that the vector AT c is nonzero for every nonzero vector c in the domain of AT . Thus the columns of AT are linearly independent. Because AT is an arbitrary submatrix of A with 2S columns, it follows that every subset of 2S columns of A is linearly independent. By Proposition 2.1, A is one-to-one on S-sparse vectors. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 144 Next we prove two propositions that provide geometric insight into RIP. In Proposition 2.4, we show that if δS is small, then the columns of A have norm approximately equal to one. In Proposition 2.6, we show that if δS is small, then pairs of distinct columns of A are approximately orthogonal. Proposition 2.4. Let A be a matrix with S-restricted isometry constant δS for some S ≥ 1. Then the norm of the ith column of A satisfies p p 1 − δS ≤ kai k2 ≤ 1 + δS . Proof. The √ standard basis vector ei is S-sparse, so by the S-restricted isometry property, kAei k2 ≤ 1 + δS kei k2 . Therefore p p kai k2 = kAei k2 ≤ 1 + δS kek2 = 1 + δS . √ Likewise, kAei k2 ≥ 1 − δS kei k2 , so p p kai k2 = kAei k2 ≥ 1 − δS kek2 = 1 − δS . The following lemma will help us prove Proposition 2.6. Lemma 2.5. The quantity kAxk22 can be expressed as X X kAxk22 = x2i kai k22 + xi xj hai , aj i. i i6=j Proof. Using the relationship between norms and inner products, we can express kAxk22 as * n + n X X X kAxk22 = hAx, Axi = xi xj hai , aj i. x i ai , x i ai = i=1 i=1 i,j The last step follows from the bilinearity of inner products. We group the summation into terms wtih i = j and terms with i 6= j: X xi xj hai , aj i = i,j n X x2i hai , ai i + i=1 X xi xj hai , aj i. i6=j Because hai , ai i = kai k22 , the right-hand side can be rewritten as kAxk22 = n X i=1 x2i kai k22 + X i6=j xi xj hai , aj i. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 145 We show that if δS is small, and A has real entries, then distinct columns of A are approximately orthogonal. Lemma 2.1 of [4] implies an analogous result for matrices with complex entries. Proposition 2.6. Let A be a matrix with real entries and S-restricted isometry constant δS . Then any two distinct columns of A satisfy |hai , aj i| ≤ 2δs . Proof. Let A be a matrix with real entries, and let δS be the S-restricted isometry constant of A. Let x be a vector with xi = xj = √12 and all other components equal to zero. Note that kxk2 = 1. By Lemma 2.5, X kAxk22 = x2i kai k22 + x2j kaj k22 + xk xl hak , al i = = 1 √ 2 2 kai k22 + k6=l 2 1 √ 2 kaj k22 + 1 √ 2 2 hai , aj i + 1 √ 2 2 haj , ai i 1 kai k22 + kaj k22 + hai , aj i. 2 By the S-restricted isometry property and the fact that kxk22 = 1, it follows that, 1−δS ≤ kAxk22 ≤ 1 + δS , so 1 − δS ≤ 1 kai k22 + kaj k22 + hai , aj i ≤ 1 + δS . 2 First we prove an upper bound on hai , aj i. By the above inequality, it follows that hai , aj i ≤ 1 + δS − 1 kai k22 + kaj k22 . 2 By Proposition 2.4, kai k22 ≥ 1 − δS and kaj k22 ≥ 1 − δS , so 1 hai , aj i ≤ 1 + δS − (1 − δS + 1 − δS ) 2 = 2δS . The lower bound on hai , aj i follows by an analogous argument: 1 kai k22 + kaj k22 2 1 ≥ 1 − δS − (1 + δS + 1 + δS ) 2 = −2δS . hai , aj i ≥ 1 − δS − We have shown that −2δS ≤ hai , aj i ≤ 2δS , so |hai , aj i| ≤ 2δS . RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 146 The key consequence of RIP is that if a matrix satisfies RIP, then `1 minimization recovers the sparsest vector in the solution set to Ax = y. The following theorem states this result. Theorem 2.7 (from [3, 4, 5]). If δ3S + 3δ4S < 2, then for any S-sparse vector x∗ that satisfies Ax∗ = y, x∗ is the unique solution to arg minkxk1 x subject to Ax = y. This theorem provides the conditions under which the `1 minimization problem, which can be solved efficiently, has the same solution as the `0 minimization problem. If a matrix satisfies the conditions of Theorem 2.7, and there is an S-sparse solution to the system of equations Ax = y, then SSR can be solved efficiently by `1 minimization. 2.5 The Structure of SSR Having looked at the details of SSR, we now summarize the structure of the SSR problem. SSR starts with an underdetermined system Ax = y. The true solution, x0 , satisfies an additional constraint: it is S-sparse. For a general underdetermined system, there may be multiple S-sparse vectors in the set of solutions, but we impose an additional condition on A, ensuring that there is at most one S-sparse solution: we require A to be one-to-one on S-sparse vectors. If x0 is an S-sparse vector such that Ax0 = y, and A is one-to-one on S-sparse vectors, then x0 is the unique solution to the optimization problem minkxk0 subject to Ax = y. In general, this optimization problem cannot be solved efficiently; it is NP hard [11]. If A satisfies RIP—a stronger condition than being one-to-one on S-sparse vectors—then the `1 minimum is the same as the `0 minimum, so SSR can be solved efficiently with the following optimization problem: minkxk1 subject to Ax = y. The table in Figure 2 summarizes the structure of SSR. When we discuss robust PCA, we will see that it has a parallel structure. 3 3.1 Robust Principal Component Analysis The Robust PCA Problem The matrix M is formed as the sum of a low-rank matrix L0 and a sparse matrix S0 : M = L0 + S0 . We observe M , but we do not observe L0 and S0 . The goal of robust PCA is to recover the matrices L0 and S0 . We can view robust PCA as solving the underdetermined system RHIT Undergrad. Math. J., Vol. 12, No. 2 1. 2. 3. Underdetermined system Extra constraint Conditions on problem instance to ensure problem is well-defined (that is, has at most one solution) 4. Actual optimization problem 5. A sufficient condition for an efficient solution 6. Tractable optimization problem Page 147 Ax = y x0 sparse A is one-to-one on S-sparse vectors minkxk0 subject to Ax = y RIP minkxk1 subject to Ax = y Figure 2: Structure of the sparse signal recovery problem. of equations M = L + S, where M is a given matrix, and L and S are the unknowns. We must impose additional conditions to ensure that L0 and S0 are the unique solution. We solve the system of equations subject to the additional constraints that L0 is low-rank and S0 is sparse. As stated, the problem is not yet well defined, because there may be multiple ways to decompose M into the sum of a low-rank matrix and a sparse matrix. As we will see later, there are conditions on L0 and S0 such that for a range of values of the rank of L0 and the sparsity of S0 , we can recover L0 and S0 using a tractable optimization problem. We now present an example that demonstrates the need for more precise conditions in defining the robust PCA problem. The example illustrates a trade-off between rank and sparsity when decomposing M : by increasing the rank of the low-rank component, we can decrease the number of nonzero entries in the sparse component. Without additional conditions on the low-rank and sparse components, it is unclear which decomposition is the “correct” one. Consider the two matrices L1 and S1 , in which L1 is rank one and S1 has four nonzero entries: 1 1 1 1 1 0 0 0 1 1 1 1 S1 = 0 2 0 0 . L1 = 1 1 1 1 0 0 3 0 1 1 1 1 0 0 0 4 Let M = L1 + S1 . If we subtract the 1 from row one, column one of S1 and add it to the corresponding entry of L1 , we get the matrices L2 and S2 : 2 1 1 1 0 0 0 0 1 1 1 1 S2 = 0 2 0 0 . L2 = 1 1 1 1 0 0 3 0 1 1 1 1 0 0 0 4 The matrix L2 has rank two, and S2 has three nonzero entries. Because L2 +S2 = L1 +S1 , we Page 148 RHIT Undergrad. Math. J., Vol. 12, No. 2 have found distinct decompositions of M that exhibit a trade-off between rank and sparsity. Accordingly, we need conditions on the low-rank matrix L0 and the sparse matrix S0 that distinguish the decomposition M = L0 + S0 from other decompositions with different ranksparsity trade-offs. We also need conditions to guarantee that there are no decompositions of M with equivalent rank and sparsity. Motivated by the need to ensure a unique solution, we will focus on the conditions under which a particular optimization problem, Principal Component Pursuit, recovers the matrices L0 and S0 from the sum M = L0 + S0 . Under the conditions of Theorem 3.4, Principal Component Pursuit recovers L0 and S0 for a range of values of the rank of L0 and the number of nonzero entries of S0 [2]. 3.2 Solving Robust PCA by Formulating an Optimization Problem To solve the robust PCA problem, Candès et al. introduced an optimization problem called Principal Component Pursuit (PCP) [2]. Under appropriate conditions on L0 and S0 , the PCP optimization problem recovers the original decomposition (L0 , S0 ). To motivate Principal Component Pursuit, we introduce a simpler optimization problem. An “ideal” optimization problem for separating a matrix into low-rank and sparse components is min rank(L) + λkSk0 subject to L + S = M, because this optimization problem searches for a low-rank matrix and a sparse matrix that add up to M . The weighting parameter λ determines which term influences the objective function more, achieving a balance between minimizing the rank of L and maximizing the sparsity of S. This optimization problem cannot in general be solved efficiently [11]. As with the sparse signal recovery problem, we adopt the strategy of replacing an objective function with a different objective function that is convex and can therefore be minimized efficiently. In the SSR problem, we replaced the `0 norm with the `1 norm. Likewise, for robust PCA we replace the `0 norm of S with the `1 norm. Note that we are not using the induced `1 norm; we are using the entrywise `1 norm, which is the sum of the absolute values of the components of the matrix. In addition to replacing kSk0 with a convex function, we replace rank(L) with a convex function. Let the singular value decomposition of L be L = U ΣV ∗ , and let σi be the ith singular value of L. The rank of a matrix is equal to the number of nonzero singular values. Instead of minimizing the number of nonzero singular values of L, we minimize the sum of the singular values, which is a convex function. In fact, the sum of the singular values is the entrywise `1 norm of Σ. The sum of the singular values is called the “nuclear norm,” denoted r X kLk∗ := σi i=1 where r is the rank of L. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 149 Replacing the rank of L with the nuclear norm, and the `0 norm of S with the `1 norm, the optimization problem becomes minkLk∗ +λkSk1 subject to M = L + S for some fixed constant λ [2]. This optimization problem is Principal Component Pursuit. 3.3 Conditions Under which Robust PCA Is Well-Defined Candès et al. impose conditions on L0 and S0 under which PCP recovers the original decomposition [2]. To understand the conditions, it is helpful to examine the motivations for them presented in [2]. Each condition is motivated by an obstacle that hinders recovery of (L0 , S0 ). Looking at the obstacles helps us understand why the conditions work and whether there are matrices that can be recovered even if they do not satisfy the conditions. 3.4 Conditions on the Low-Rank Component L0 The first obstacle occurs if the low-rank matrix L0 is sparse. Because we are trying to separate a matrix M into a low-rank component and a sparse component, the separation is not well-defined if the low-rank component is also sparse. We show that if the low-rank component is also sparse, then there are multiple ways to decompose the matrix M into the sum of a low-rank matrix and a sparse matrix. Proposition 3.1. Suppose M is an n × n matrix with M = L0 + S0 , where L0 is a rank-r matrix and S0 has S nonzero components. If there is a k-sparse vector in the column space of L0 , then we can express M as M = L1 + S1 , where L1 has rank r − 1 and S1 has at most kn additional nonzero components. Proof. Let v1 be a k-sparse vector in the column space of L0 . Extend v1 to an orthonormal basis {v1 , · · · , vr } of the column space of L0 . Any vector in the column space of L0 can be written as a vector in span(v1 ) plus a vector in span(v2 , · · · , vr ). Let li be the ith column of L0 , and let si be the ith column of S0 . We write the column li as a vector in span(v1 ) plus a vector in span(v2 , · · · , vr ). That is, let li = xi + wi , where xi ∈ span(v1 ) and wi ∈ span(v2 , · · · vr ). Given the above notation, we define the matrices L1 and S1 as follows. The ith column of L1 is wi , and the ith column of S1 is xi + si . The ith column of L1 and the ith column of S1 sum to form the ith column of M , so that L1 + S1 = M : wi + (xi + si ) = (wi + xi ) + si = li + si = mi . Because each wi is in the span of v2 , · · · , vr , it follows that L1 has rank r − 1. Because xi is in the span of v1 , and v1 is k-sparse, it follows that xi is k-sparse. The original S0 is S-sparse, and S1 consists of k-sparse vectors added to the columns of S0 , so S1 has at most S + kn nonzero components. We have thus provided a way to decompose M into a matrix of rank r − 1 plus a matrix with at most S + kn nonzero components. Page 150 RHIT Undergrad. Math. J., Vol. 12, No. 2 To rule out sparse matrices L0 , we impose a condition on L0 that ensures that there is no sparse vector in the column space of L0 . The conditions we impose are called the incoherence conditions [2]. The incoherence conditions are expressed in terms of the singular value decomposition of L0 . Stating the incoherence conditions requires the following notation. The singular value decomposition of L0 is L0 = U ΣV ∗ . Let r be the rank of L0 . Let Ur be the matrix consisting of the first r columns of U , and let Vr be the matrix consisting of the first r columns of V . The incoherence conditions impose three constraints on the left and right singular vectors of L0 : 1. The projections of the standard basis vectors onto the column space of L0 are small. 2. The projections of the standard basis vectors onto the row space of L0 are small. 3. The projections of the rows of Ur onto the rows of Vr are small. We will show that the incoherence conditions rule out matrices L0 with sparse vectors in their column spaces. Proposition 3.3 states that if the projections of the standard basis vectors onto a unit vector x are small, then x is not sparse. Thus Proposition 3.3 implies that if the first incoherence condition holds, then the columns of L0 are not sparse. The following lemma helps in the proof of Proposition 3.3. Lemma 3.2. Let x be a vector with `2 norm equal to 1 and at most S nonzero components. Then x has a component whose absolute value is at least √1S . Proof. Let x be a vector with at most S nonzero components, and let I be the set of indices corresponding to the nonzero components of x. If the magnitude of each nonzero component √ of x is strictly smaller than 1/ S, then v uX sX 2 s 2 u 1 1 2 √ kxk2 = |xi | < t ≤ S √ = 1. S S i∈I i∈I The first√inequality holds because the nonzero components of x have magnitude strictly less than 1/ S. The second inequality holds because x has at most S nonzero components, so there are at most S terms in the summation. We have shown that if x has at most S nonzero √ components, and the magnitudes of the nonzero components are strictly smaller than 1/ S, then kxk2 is strictly smaller than 1. It follows that if kxk2 = 1 and x has at most √ S nonzero components, then x has at least one component whose magnitude is at least 1/ S. We now show that if the projections of the standard basis vectors onto a unit vector x are sufficiently small, then x is not sparse. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 151 Proposition 3.3. If the magnitude of the projection of each standard basis vector onto x is smaller than √1S and kxk2 = 1, then x has at least S + 1 nonzero components. Proof. Because ei is a unit vector, |hei , xi| = |xi | is the magnitude of the projection of ei onto x. If |xi | < √1S for each component of x, then because kxk2 = 1, it follows from Lemma 3.2 that x has at least S + 1 nonzero components. The first incoherence condition guarantees that the projections of the standard basis vectors onto the column space of L0 are small. In particular, the projections of the standard basis vectors onto the columns of L0 are small. Therefore, by Proposition 3.3, the columns of L0 are not sparse, so the incoherence conditions ensure that L0 is not sparse. We now present the formal definition of the incoherence conditions. Recall that L0 is an m × n matrix of rank r with SVD L0 = U ΣV ∗ . The matrices Ur and Vr are the matrices consisting of the first r columns of U and V , respectively. The incoherence parameter µ is a property of the matrix L0 . It is the smallest value that satisfies all three inequalities (from [2]): µr , m µr maxkVr∗ ei k22 ≤ , i r n µr . kUr Vr∗ k∞ ≤ mn maxkUr∗ ei k22 ≤ i (1) (2) (3) The third inequality uses the entrywise `∞ norm, written k·k∞ , which is the largest absolute value of any entry in the matrix. Satisfying the incoherence conditions means having a small value of µ. For the first incoherence condition (the first inequality), a small value of µ means that the projections of the standard basis vectors onto the column space of L0 are small. For the second incoherence condition, a small value of µ means that the projections of the standard basis vectors onto the row space of L0 are small. For the third incoherence condition, a small value of µ means that the projections of the rows of Ur onto the rows of Vr are small. 3.5 Conditions on the Sparse Component S0 The second obstacle that undermines our ability to recover the original matrices (L0 , S0 ) is if S0 is low-rank. If the sparse component is also low-rank, then there is ambiguity regarding whether columns of the sparse component should actually be included within the low-rank component. Suppose that L0 has rank r and S0 has rank k. Then the matrix L1 = L0 + S0 has rank at most r + k. Therefore the decomposition M = L1 + S1 , where S1 is the zero matrix, is another decomposition of M into a low-rank component and a sparse component. More generally, if a subset of columns of S0 spans a k-dimensional subspace, then those columns can be subtracted from S0 and added to the corresponding columns of L0 to produce a new RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 152 decomposition (L1 , S1 ). The matrix L1 has rank at most r + k, and S1 has fewer nonzero components than S0 (because a subset of columns of S0 has been replaced with zeros). This argument demonstrates that if S0 is low-rank, then we can produce a different decomposition of M into the sum of a low-rank matrix and a sparse matrix. Because we do not know the rank of L0 in advance, and we do not know the number of nonzero components of S0 in advance, we have no way to distinguish the original decomposition (L0 , S0 ) from the other decomposition (L1 , S1 ). We now discuss an assumption that rules out matrices S0 that are low-rank. In [2], the authors adopt the assumption that the locations of the nonzero entries in S0 have been chosen randomly. Let S be the number of nonzero entries in the matrix S0 . We assume that each set of S entries in the matrix has an equal probability of being chosen as the subset of nonzero entries. Because the locations of the nonzero entries are randomly distributed, the probability that S0 is low-rank is small. In practice, the locations of the errors are not necessarily chosen randomly, but this assumption makes the analysis easier. Candès et al. show that if this assumption is satisfied (together with the incoherence conditions), then the PCP optimization problem recovers the original decomposition [2]. 3.6 The Main Theorem The main theorem of [2] states that if the above conditions on L0 and S0 are satisfied, then PCP recovers the original decomposition with high probability. (Recall that the only source of randomness is the assumption that the locations of the nonzero entries of S0 are randomly chosen.) The theorem is as follows: Theorem 3.4 (from [2]). Let L0 and S0 be n×n matrices that satisfy conditions (1) – (3). Let the locations of the nonzero entries of S0 be uniformly distributed. Then PCP succeeds with high probability if the following two inequalities hold, rank(L0 ) ≤ ρr µ−1 kS0 k0 ≤ ρs n2 , n , log 2 (n) (4) (5) where ρr and ρs are constants. The first inequality states that the rank of L0 is smaller than some quantity that depends on the incoherence parameter. Hence the first inequality ensures that the rank is low and the incoherence conditions are satisfied. The second inequality states that the number of nonzero components of S0 is smaller than some fraction of the number of entries. These conditions are sufficient for PCP to successfully recover the original decomposition, but they may not be necessary. There may be matrices (L0 , S0 ) that are successfully recovered by PCP even though they do not satisfy the above conditions. We investigate the conditions in more detail in the next section. RHIT Undergrad. Math. J., Vol. 12, No. 2 3.7 Page 153 The Structure of Robust PCA We summarize the structure of the robust PCA problem and highlight its parallels with sparse signal recovery. The goal of robust PCA is to separate a matrix M into a low-rank matrix L0 and a sparse matrix of errors S0 . This goal corresponds to solving an underdetermined system of equations M = L + S, with the additional constraints that L is “low-rank” and S is “sparse.” Because there can be multiple decompositions of a matrix into low-rank and sparse components, these constraints must be made more precise for the robust PCA problem to be well-defined. Instead of precisely defining the robust PCA problem, we focus on the conditions under which Principal Component Pursuit recovers matrices L0 and S0 from their sum M = L0 + S0 . To separate the matrix M into low-rank and sparse components, we could formulate the following optimization problem min rank(L) + λkSk0 subject to M = L + S. This objective function is “ideal” because it is conceptually simple—it is a direct formalization of our intuition that L0 is low-rank and S0 is sparse. However, it has no efficient solution, so Candès et al. replace it with the tractable convex optimization problem Principal Component Pursuit (PCP) [2]: minkLk∗ + λkSk1 subject to M = L + S. Under conditions on L0 and S0 imposed by Theorem 3.4, which ensure that L0 is not sparse and S0 is not low-rank, PCP recovers L0 and S0 from their sum. 4 Principal Component Pursuit and the Incoherence Conditions We now investigate the conditions under which PCP successfully solves the robust PCA problem. The incoherence conditions discussed in the previous section are sufficient to guarantee that the PCP algorithm recovers the original low-rank matrix. The motivation behind the incoherence conditions is to ensure that the low-rank matrix L0 is not sparse. While the motivation is intuitively clear, it is difficult to verify that the incoherence conditions are satisfied without observing L0 . Furthermore, the incoherence conditions may not be necessary. We present experiments that demonstrate that recovery is possible even if the conditions of Theorem 3.4 are violated. Ideally, we would like to identify conditions on (L0 , S0 ) that are 1. necessary and sufficient, 2. easily checked, 3. satisfied by a broad class of matrices. Page 154 RHIT Undergrad. Math. J., Vol. 12, No. 2 Robust PCA 1. Underdetermined system L+S =M 2. Extra constraints L0 low-rank, S0 sparse 3. Conditions on problem in- Existence of a unique destance to ensure problem is composition of M into M = well defined (that is, has at L+S such that rank(L) ≤ r most one solution) and kSk0 ≤ s, with r and s fixed constants 4. “Ideal” optimization prob- min rank(L) + λkSk0 sublem ject to L + S = M 5. Tractable optimization minkLk∗ + λkSk1 subject to problem L+S =M 6. Conditions under which Incoherence Conditions the tractable optimiza- (Theorem 3.4) tion problem recovers the correct solution SSR Ax = y x0 sparse A is one-to-one on S-sparse vectors minkxk0 subject to Ax = y minkxk1 subject to Ax = y Restricted Isometry Property Figure 3: This table presents the structure of robust PCA and SSR, illustrating the parallels between the two problems. In Section 4.1, we consider whether randomly generated matrices satisfy the incoherence conditions. In Section 4.2, we examine the relationship between the incoherence conditions and the success of PCP. 4.1 Incoherence of Randomly Generated Matrices To investigate whether randomly generated matrices satisfy the incoherence conditions, use the following procedure. Suppose X is an n × r matrix with each entry generated independently from a normal distribution of mean 0 and variance √1n . Suppose Y is also an n × r matrix with each entry generated independently from a normal distribution with mean 0 and variance √1n . Then XY T is an n × n matrix of rank r. This observation results in a method to generate rank-r matrices. To obtain an n × n matrix of rank r, we generate matrices X and Y as described above, then take the product XY T . We test whether matrices generated by the above procedure satisfy the incoherence conditions. All experiments were performed in Matlab. For every dimension between 10 and 100 and every rank between 1 and the dimension, we generated a matrix according to the above procedure. For each matrix, we determined the incoherence parameter—the smallest value of µ that satisfies the incoherence inequalities (1) – (3). The plot in Figure 4 shows the incoherence parameter as a function of rank and dimension. The only meaningful values in the plot are the values below the diagonal, because the rank cannot exceed the dimension of the matrix. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 155 100 90 120 80 100 70 Rank 60 80 50 60 40 30 40 20 20 10 0 0 10 20 30 40 50 60 Dimension 70 80 90 100 0 Figure 4: Incoherence parameter as a function of rank and dimension. The incoherence parameter of a matrix is the smallest value of µ that satisfies inequalities (1) – (3). For each dimension between 10 and 100 and each rank between 1 and the dimension, we generated a matrix and calculated the incoherence parameter of the matrix. Each matrix was generated as a product XY T , where X and Y are n × r matrices with entries generated from a normal distribution of mean 0 and variance √1n . As shown in Figure 4, there is a region of high incoherence parameter corresponding Student Version of MATLAB to a horizontal band at the bottom of the plot. For matrices with rank between 1 and 5, the incoherence parameter tends to be markedly greater than for matrices with higher rank. The plot suggests that the incoherence parameter is higher for low-rank matrices but decays rapidly as the rank increases. This is expected, because in each of the three inequalities, the incoherence parameter is inversely proportional to the rank: m maxkUr∗ ei k22 , r i n µ ≥ maxkVr∗ ei k22 , r i mn µ≥ kUr Vr∗ k2∞ . r µ≥ The decay in the incoherence parameter as the rank increases suggests that as the rank increases, the factors maxi kUr∗ ei k2 , maxi kVr∗ ei k2 , and kUr Vr∗ k2∞ do not increase as rapidly as the rank. A question that arises from the plot is which of the three inequalities requires µ to have the highest value, thus serving as the binding constraint. We performed an experiment to test which constraint is binding for matrices with normally distributed entries. Following is RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 156 100 15 90 80 70 10 Rank 60 50 40 5 30 20 10 0 0 10 20 30 40 50 60 Dimension 70 80 90 100 0 Figure 5: “Left incoherence parameter.” Minimum value of µ that satisfies the first incoherence inequality, as a function of rank and dimension. We generated one matrix for each value of dimension and rank. The matrices were generated using the procedure described in Figure 4. a description of the experiment, with the results displayed in Figure 5. For each dimension between 10 and 100, and each rank between 1 and the dimension, we generated a matrix using of MATLAB the same procedure as for Figure 4. For Student each Version matrix, we calculated the smallest value of µ that satisfies the first incoherence inequality. We refer to that value as the “left incoherence parameter,” because the first incoherence inequality involves the left singular vectors of the matrix. The plot in Figure 5 shows the left incoherence parameter as a function of rank and dimension for the matrices that we generated. This plot has the same form as the plot in Figure 4. That is, there is a horizontal band at the bottom of the plot corresponding to a region of high left incoherence parameter, and the left incoherence parameter decays as the rank increases. But the magnitudes of the left incoherence parameters in Figure 5 are smaller than the magnitudes of the incoherence parameters in Figure 4. This suggests that for “generic” matrices, the first incoherence inequality is not binding. A similar plot for the second incoherence condition suggests that the second incoherence condition is not binding either. Hence the third inequality is binding. The third incoherence inequality requires that the projections of the rows of Ur onto the rows of Vr are small. To investigate the incoherence conditions in more detail, we vary the rank and dimension separately. For every rank between 1 and 200, we generated a 200 × 200 matrix using the procedure described in Figure 4. Figure 6 shows the incoherence parameter as a function of the rank. After a sharp drop-off, the plot fluctuates around a level value. A future question RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 157 80 70 Incoherence Parameter 60 50 40 30 20 10 0 20 40 60 80 100 Rank 120 140 160 180 200 Figure 6: Incoherence parameter as a function of rank, with the dimension fixed at 200×200. The matrices were generated as products XY T , where X and Y are 200 × r matrices with entries generated from a normal distribution with mean 0 and variance √1n . to investigate is why the incoherence parameter does not seem to depend on the rank after the initial drop-off. We also vary the dimension while fixing the rank. For every dimension between 10 and 200, we generated a rank 10 matrixStudent using the above procedure. Figure 7 shows the Version of MATLAB incoherence parameter as a function of the dimension. The incoherence parameter seems to increase as a function of dimension, as expected from the incoherence inequalities. Understanding the incoherence of randomly generated matrices gives insight into whether “generic” matrices satisfy the incoherence conditions. 4.2 Incoherence Conditions and the Success of PCP We investigate how the performance of PCP depends on the incoherence parameter of the low-rank matrix. For each dimension between 10 and 100 and each rank between 1 and the dimension, we generate a matrix using the procedure described in Figure 4. For each matrix we calculate the amount by which the matrix violates the first assumption of the main PCP theorem, Theorem 3.4. The first assumption states that rank(L0 ) ≤ ρr µ−1 n . log 2 (n) The amount by which the inequality is violated is equal to the left-hand side minus the right-hand side. We set the constant ρr equal to 1. Figure 8 shows the PCP assumption RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 158 60 Incoherence Parameter 50 40 30 20 10 0 0 20 40 60 80 100 120 Dimension 140 160 180 200 Figure 7: Incoherence parameter as a function of dimension, with the rank fixed at 10. The matrices were generated as products XY T , where X and Y are n × 10 matrices with entries generated from a normal distribution with mean 0 and variance √1n . violations as a function of rank and dimension. All of the matrices exhibit violations; none satisfy the first assumption of the main theorem. As the rank increases, the amount by which the inequality is violated increases. The violations do not seem to depend on the dimension. We performed tests to reveal whetherStudent successful occur despite violations of the Version ofrecoveries MATLAB sufficient conditions of Theorem 3.4. We generated 30×30 matrices L0 and S0 , and we tested whether the PCP objective function correctly recovers L0 and S0 from their sum M = L0 +S0 . The low-rank matrices L0 were generated using the procedure described in Figure 4. For the sparse matrices S0 , the locations of the nonzero entires were chosen uniformly at random, the signs of the entries were chosen randomly (positive with probability 1/2), and the magnitudes of the nonzero entries were equal to 10. To solve the PCP optimization problems we used CVX, a software package for convex optimization [8]. Figure 9 shows the proportion of successful recoveries as a function of rank(L0 )/30 and of the proportion of nonzero entries in S0 . For each value of the rank of L0 and the proportion of nonzero entries in S0 , we conducted ten trials of PCP. A trial is considered successful if the relative error of the recovered low-rank matrix is sufficiently small, measured with the Frobenius norm. That is, a trial is successful if the following condition holds kL0 − Lrec kF < 0.001, kL0 kF where Lrec is the low-rank matrix recovered by PCP. Figure 9 shows that despite the violations of the PCP conditions, there are successful recoveries. As the rank of L0 increases RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 159 Rank 100 90 90 80 80 70 70 60 60 50 50 40 40 30 30 20 20 10 10 0 0 10 20 30 40 50 60 Dimension 70 80 90 100 0 Figure 8: PCP assumption violations as a function of rank and dimension. For each dimension between 10 and 100 and each rank between 1 and the dimension, we generated a low-rank matrix L0 . We measured the violation of inequality (4) in Theorem 3.4. That is, the plot displays the value of rank(L0 )µ − ρr logn2 n . and the proportion of nonzero entries in S0 increases, the fraction of successful recoveries decreases. Student Version of MATLAB 5 Conclusion The above experiments suggest that successful recoveries are possible for matrices that violate the incoherence conditions. Future work could investigate the role of each of the three incoherence inequalities separately. In particular, what is the role of the third incoherence inequality, and why is it binding for the matrices generated in the above experiments? Because the conditions of Theorem 3.4 are sufficient but not necessary, future work could search for weaker conditions that are nonetheless sufficient, as well as conditions that are easier to check. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 160 1 0.4 0.9 Proportion of nonzero entries in S 0.8 0.3 0.7 0.6 0.5 0.2 0.4 0.3 0.1 0.2 0.1 0 0.05 0.1 0.15 0.2 0.25 rank/n 0.3 0.35 0.4 0 Figure 9: Proportion of successful recoveries as a function of the rank of the low-rank matrix L0 and the fraction of nonzero entries in the sparse matrix of errors S0 , for 30 × 30 matrices. Ten trials were performed for each value of the rank of L0 and sparsity of S0 . The lowrank matrices L0 were generated using the procedure described in Figure 4. For the sparse matrices S0 , the locations of the nonzero entires were chosen uniformly at random, the signs of the entries were chosen randomly (positive with probability 1/2), and the magnitudes of the nonzero entries were equal to 10. We used CVX to solve the PCP optimization problems. Student Version of MATLAB Trials were judged by the relative error in recovering the low-rank matrix, measured by the −Lrec kF < 0.001, where Lrec is the Frobenius norm. A trial was considered successful if kL0kL 0 kF matrix recovered by PCP. RHIT Undergrad. Math. J., Vol. 12, No. 2 Page 161 References [1] S.P. Boyd and L. Vandenberghe. Convex optimization. Cambridge Univ Pr, 2004. [2] E. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component analysis? Arxiv preprint arXiv:0912.3599, 2009. [3] E. Candès, M. Rudelson, T. Tao, and R. Vershynin. Error correction via linear programming. ANN IEEE SYMP FOUND, 2005. [4] E.J. Candès. The restricted isometry property and its implications for compressed sensing. Comptes Rendus Mathematique, 346(9-10):589–592, 2008. [5] E.J. Candès and T. Tao. Decoding by linear programming. IEEE Trans. Inform. Theory, 51(12):4203–4215, 2005. [6] D.L. Donoho and P.B. Stark. Uncertainty principles and signal recovery. SIAM J. Appl. Math, 49(3):906–931, 1989. [7] C. Eckart and G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3):211–218, 1936. [8] M. Grant and S. Boyd. CVX: Matlab software for disciplined convex programming, version 1.21. http://cvxr.com/cvx, July 2010. [9] H. Hotelling. Analysis of a complex of statistical variables into principal components. Journal of educational psychology, 24(6):417, 1933. [10] I.T. Joliffe. Principal component analysis, 1986. [11] S. Muthukrishnan. Data streams: Algorithms and applications. Now Publishers, MA, 2005.