Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
53,007 characters · 10 sections · 0 citation commands
Robust Principal Component Analysis with Non-Sparse Errors
A low-rank component in high dimensional data sets is often the object of interest. In asset return analysis, for example, a low-rank matrix represents systematic risks (Ross, 1976). In psychology, the main personality traits form a low-rank matrix (e.g., Cattell, 1978 and Goldberg, 1990). In background/foreground separation, by stacking the pixels of each frame of a video in a column vector, the static background is a rank one component in the resulting matrix because it stays unchanged across frames, see Bouwmans et al.\enspace(2017) for a survey. In gene expression prediction, the gene expression values may form a low-rank matrix because genes act in groups and at the expression levels they are interdependent (Kapur et al.\enspace, 2016).
To fix ideas, in this paper, we assume the data matrix, $Y$, is of $N\times T$ dimension and consists of a low-rank component $L_{0}$, namely
where $rank(L_{0})=r$ is small but unknown; $L_{0}$ can be random or deterministic. The magnitude of the its elements are allowed to diverge to infinity with $N$ and $T$ grow; $Z_{0}$ is a random error matrix with median zero entries that have positive densities around $0$.
We propose to estimate $L_{0}$ using a variant of the Principal Component Pursuit (PCP), introduced and studied by Cand\`es et al.\enspace(2011), Chandrasekaran et al.\enspace (2011) etc. We show the estimator is consistent for $L_{0}$ in the Frobenius norm under certain conditions. This is the first time that consistency is established with continuous and potentially fat-tailed random errors. Formally, the estimator $\hat{L}$ is as follows
where $||\cdot||_{*}$ denotes the nuclear norm. $||\cdot||_{1}$ and $||\cdot||_{\infty}$ are the $\ell_{1}$ norm and the $\ell_{\infty}$ norm of a matrix seen as a long vector. Both $\lambda$ and $\alpha$ are $N,T$ dependent. In particular, $\alpha$ can grow with $N$ and $T$. We call $\hat{L}$ the Bounded Principal Component Pursuit (BPCP) as compared to the standard PCP, it has an additional constraint bounding the max entry of the estimator by $\alpha$.
As a preview of how the estimator works, first note that the nuclear norm is the convex envelope of the rank of a matrix because by definition, the nuclear norm is the $\ell_{1}$ norm of the vector of singular values while the rank is its $\ell_{0}$ norm, that is, the number of nonzero elements in the vector. Thus, minimizing the nuclear norm is a convex surrogate of rank minimization.
The other term in the objective function is the $\ell_{1}$ penalty for the residuals to induce robustness. The intuition is analogous to the LAD (least absolute deviation) estimator for linear regression; it is well-known that minimizing the sum of absolute deviations is robust to fat-tailed errors. It turns out for the BPCP estimator to work, the errors essentially only need to have zero median, just like the standard LAD estimator.
Finally, the constraint in (ref) is needed for technical reasons. We allow $\alpha$ to be $N,T$ dependent and can go to infinity as $N$ and $T$ increase. We will be precise about the allowed rate of such divergence. As it turns out, the restriction is actually very mild because in many stochastic models of $L_{0}$, $||L_{0}||_{\infty}$ diverges slower than the rate allowed with high probability. Therefore, imposing the constraint is without loss of generality in these models.
This paper adds to the theory of PCP and some of its variants, developed in Cand\`es et al.\enspace(2011), Chandrasekaran et al.\enspace (2011), Ganesh et al.\enspace(2010), Zhou et al.\enspace(2010), etc. In Cand\`es et al.\enspace(2011), it is assumed that $Z_{0}$ follows a Bernoulli model, i.e., each element in $Z_{0}$ is equal to $0$ with probability $1-\rho$. They show that when $1-\rho$ is large enough, PCP can \textit{exactly} recover both $L_{0}$ and $Z_{0}$ with high probability. The proof hinges on the existence of a matrix called \textit{dual certificate}, which relies on the sparsity of $Z_{0}$. Ganesh \textit{et al}.\enspace(2010) generalize the result by allowing for an arbitrarily small but fixed $1-\rho$ and show a dual certificate still exists. However, their results are not applicable to a continuously distributed error matrix because then for any entry $Z_{0,it}$, $P(Z_{0,it}=0)=0$ by definition, and thus $1-\rho=0$. Zhou \textit{et al}.\enspace(2010) study a variant of PCP called the \textit{Stable Principal Component Pursuit} (SPCP). In their model, $Y$ is known to consist of a low rank, a sparse, and a dense component. They minimize a similar objective function over both $L$ and $Z$ with a constraint bounding the difference between the sum of them and $Y$ in the Frobenius norm. They show that the Frobenius norm of the estimation error is bounded by the Frobenius norm of the dense component multiplied by $\max\{N,T\}$. This bound is evidently too large for consistency. Hsu \textit{et al}.\enspace(2011) change the objective function by adding a squared Frobenius norm penalty for the difference between $Y$ and the sum of the low-rank and the sparse components. They prove nuclear norm consistency for the low-rank matrix provided that the sparse component has an increasing fraction of zeros and the dense errors are Gaussian (see the second example in their section D). All the existing work in the above literature require $Z_{0}$ to have a certain fraction of entries to be $0$ with positive probability. In contrast, in this paper all entries in $Z_{0}$ can be nonzero almost surely. Agarwal \textit{et al}.\enspace(2012) study a broad class of models allowing $Y$ to be determined by a general operator of $L_{0}$ and $Z_{0}$, which are not necessarily to be exactly low-rank or sparse. Similar to Hsu \textit{et al}.\enspace(2011), their objective function also has the additional Frobenius norm penalty as they allow for the existence of an additional noise component whose operator norm is not too large. Their results are more comparable with ours because they also allow for a $Z_{0}$ with all entries nonzero. However, to obtain consistency, $Z_{0}$ needs to be approximately sparse, i.e., the fraction of the entries that are large in magnitude needs to be shrinking to $0$ and the sum of the absolute values of the rest entries is $o_{p}(NT)$. This condition rejects many random models for $Z_{0}$, especially if $Z_{0}$ has very fat tails. By contrast, this paper only focuses on the linear decomposition model (ref) and under stronger assumptions including a probabilistic model for $Z_{0}$, consistency is established even if all entries of $Z_{0}$ are nonzero and most of them are large in magnitude; we do not put any restrictions on the entries' tail distributions, so long as they have zero median and positive densities around $0$.
This paper also lies in the broader literature of estimating low-rank components in various settings. The following is only a small portion of many contributions in this literature. Tao and Yuan (2011), Xu et al.\enspace(2012), Wong and Lee (2017) and Brahma et al.\enspace(2018) study the case where both a sparse component and a dense noise component exist besides the low-rank component in the decomposition. Xu et al.\enspace(2012) replaces the $\ell_{1}$ norm in the objective function with the $\ell_{1,2}$ norm. Wong and Lee (2017) changes it to the Huber loss function. Brahma et al.\enspace(2018) allows other general forms of penalty, such as SCAD penalty. To achieve consistency, they all need the noise matrix to have small norms. Wright et al.\enspace(2013) study the noiseless case but $Y$ is a compressive measurement of $L_{0}+Z_{0}$. Bach et al.\enspace(2008) and Negahban and Wainwright (2011) study regression models with low-rank coefficient matrices. Bai and Li (2012) impose a factor structure on the low-rank component and estimate it using MLE. Xie and Xing (2014) consider the Principal Component Analysis explicitly assuming Cauchy noise under the MLE framework and their minimization problem is nonconvex. Cai \textit{et al}.\enspace(2010), Cand\`es and Recht (2009), and Cand\`es and Tao (2010), among others, study matrix completion problem for a low-rank matrix. Xia \textit{et al}.\enspace(2018) develop methods for estimating low-rank tensors.
The rest of the paper is organized as follows. Section (ref) introduces the main technical tool we developed for this problem, named Bernoulli device. This device decomposes $Z_{0}$ into a matrix $D_{0}$ only containing $Z_{0}$'s small entries and a matrix $S_{0}$ with large entries of $Z_{0}$. Section (ref) extends results in Ganesh et al.\enspace(2010) by showing a dual certificate exists even if the fraction of zero entries in $S_{0}$ decreases to $0$ slowly. Section (ref) presents the key condition for consistency derived from the optimality condition by using the dual certificate and exploiting the complementary structure of $D_{0}$ and $S_{0}$ endowed by the Bernoulli device. Section (ref) states and proves the main theorems of the paper. Simulation results are demonstrated in Section (ref). Section (ref) concludes. Proofs of some of the lemmas are contained in the Appendix.
Throughout, $||\cdot||_{*}$, $||\cdot||_{1}$, $||\cdot||_{F}$ and $||\cdot||_{\infty}$ denote the nuclear norm, the $\ell_{1}$ norm, the Frobenius norm, and the max norm of a matrix. $||\cdot||$ denotes the Euclidean norm of a vector, or the operator norm of a matrix or an operator. For the same matrix, $||\cdot||\leq||\cdot||_{F}\leq||\cdot||_{*}$ and $||\cdot||_{F}\leq||\cdot||_{1}$. For two generic scalars $a$ and $b$, denote $a\land b\equiv \min\{a,b\}$ and $a\lor b\equiv \max\{a,b\}$. For any positive sequences $a$ and $b$, $a\asymp b$ means there exist $0<c_{1}\leq c_{2}<\infty$ such that $c_{1}a\leq b\leq c_{2}a$. For any matrices $A$ and $B$ of the same size, $A\circ B$ is the componentwise product of A and B. For any two random objects $X$ and $Y$, denote independence between $X$ and $Y$ by $X\protect\mathpalette{\protect\independenT}{\perp} Y$. Finally, $C$, $C'$, $C_{1}$ and $C_{2}$ denote generic positive constants that may be different in different uses.
The object of interest in this paper is
where $\hat{L}$ is defined in (ref). To bound this quantity and obtain consistency, we follow the idea in Cand\`es et al.\enspace(2011), Ganesh et al.\enspace(2010) and Zhou et al.\enspace(2010) to use a dual certificate, a matrix which will be defined in the next section, to derive the optimality condition for (ref). This condition will then yield a bound for (ref). The main theoretical challenge is that in the first place, the existence of a dual certificate hinges on the existence of zero entries (with positive probability) in the error matrix, a luxury we do not have in this paper. The key idea is to decompose $Z_{0}$ into $S_{0}+D_{0}$ in such a way that i) a large enough fraction of entries in $S_{0}$ are $0$ with positive probability to guarantee the existence of a dual certificate, and ii) that fraction cannot be too large, on the other hand, so that (ref) can be bounded by a function of $||D_{0}||_{1}$ that converges to $0$ in probability. We begin with constructing this decomposition using a Bernoulli device.
We first introduce the following assumption.
Note that equicontinuity in c) can be replaced with continuity if we strengthen the independence condition in part b) to be i.i.d.
Under Assumption (ref), let $\{\delta\}$ be a positive sequence such that $\delta\to 0$ as $N,T\to\infty$. Let $(\underline{\gamma}_{it},\bar{\gamma}_{it})$ be a pair of constants satisfying
Assumption (ref) b) and c) guarantees the existence and uniqueness of such a pair for large enough $N,T$ while $\delta$ approaches $0$.
With $\underline{\gamma}_{it}$ and $\bar{\gamma}_{it}$, let $M$ be an $N\times T$ matrix whose entries are defined by
where $\mathbbm{1}(\cdot)$ is the indicator function. Then let $D_{0}=M\circ Z_{0}$ and $S_{0}=Z_{0}-M\circ Z_{0}$ and we have
Under Assumption (ref), $D_{0}$ and $S_{0}$ have the following properties:
The Bernoulii device $M$ thus delivers a pair $(D_{0},S_{0})$ that achieves the two goals described in the beginning of this section. First, as will be seen, items 2 and 3 guarantee the existence of a dual certificate under certain conditions. Second, by Hoeffding's inequality, the order of $||D_{0}||_{1}$ is no greater than $C'NT\delta^{2}$ in high probability provided that $\delta$ converges to $0$ at an appropriate rate. Note this holds regardless of how the distribution of the original error $Z_{0,it}$ behaves except the requirements for the zero-median and the positive and continuous density at $0$.
In this section, we treat $S_{0}$ as the “sparse” error matrix and show that although $\delta$ decreases to $0$, a dual certificate that is similar to Ganesh et al.\enspace(2010) exists. As mentioned in the introduction, Ganesh et al.\enspace(2010) show a dual certificate exists for any small yet fixed $\delta$. We extend their results by carefully choosing the rate of $\delta$, $\lambda$ and other constants in their proof. We closely follow their construction of the dual certificate but for completeness, we record it here and shall indicate where the construction needs to be modified to handle a shrinking $\delta$ by construction.
First we need an identification condition to guarantee $L_{0}$ to be non-sparse so that it is distinguishable from $S_{0}$. We adopt the incoherence condition in Cand\`es and Recht (2009), Cand\`es et al.\enspace(2011), Ganesh et al.\enspace(2010), etc. Besides it, as there is an additional constraint in (ref), we need a condition to guarantee $L_{0}$ to be a feasible solution. Let $U\Sigma V^{*}$ be a singular value decomposition of $L_{0}$, where $U$ and $V$ are $N\times r$ and $T\times r$ matrices of left and right singular vectors and $\Sigma$ is an $r\times r$ diagonal matrix with singular values in descending order on its diagonal.
Inequality (ref) in Assumption (ref) is the incoherence condition, stating that the singular vectors of $L_{0}$ are well-spread. A direct and useful consequence of (ref) is that
by noticing that $||UV^{*}||_{\infty}=\max_{it}|\sum_{k=1}^{r}U_{ik}V_{tk}|\leq\sqrt{\sum_{k=1}^{r}|U_{ik}|^{2}}\cdot \sqrt{\sum_{k=1}^{r}|V_{tk}|^{2}}\leq C\frac{\mu r}{\sqrt{NT}}$ where the first inequality follows from the Cauchy-Schwarz inequality and the second is from (ref). Here $\mu$ characterizes how coherent the singular vectors are with the canonical bases. It can be $N,T$ dependent and is allowed to diverge to $\infty$. Cand\`es and Recht (2009) provide examples where $\mu=O(\log(N\lor T))$ and one of them is the random orthogonal model in which the columns in $U$ and $V$ are sampled uniformly among all families of $r$ orthonormal vectors independently of each other. Fan et al.\enspace(2018) also give an example where $\mu=O\big(\sqrt{\log(N\lor T)}\big)$.
Inequality (ref) is an inclusion assumption which implies $L_{0}$ is a feasible solution with probability approaching 1. It restricts the magnitude of the maximal entry in $L_{0}$. Again $\alpha$ is allowed to increase to $\infty$ with $N$ and $T$. Note that (ref) and (ref) imply that $L_{0}$'s largest singular value $\sigma_{1}\leq \frac{\alpha}{C}\frac{\sqrt{NT}}{\mu r}$ because $L_{0}=\sum_{k=1}^{r}u_{k}v^{*}_{k}\sigma_{k}$ while $UV^{*}=\sum_{k=1}^{r}u_{k}v^{*}_{k}$, where $u_{k}$ and $v_{k}$ are the $k$th column of $U$ and $V$, respectively.
Before we define the dual certificate, it is useful to introduce some notations.
Let $\Phi$ be the linear space of matrices
and let its orthogonal complement be $\Phi^{\perp}$. Denote the linear projection onto $\Phi$ and $\Phi^{\perp}$ by $\mathcal{P}_{\Phi}$ and $\mathcal{P}_{\Phi^{\perp}}$, respectively. Then it can be shown that for any $N\times T$ matrix $R$ (e.g. Cand\`es and Recht (2009)),
and
Let $\Omega$ be the support of $S_{0}$, i.e., $\Omega=\{(i,t):S_{0,it}\neq 0\}$. With a slight abuse of notation, we also denote the linear space of matrices supported on $\Omega$ by $\Omega$. The projection onto this space is denoted by $\mathcal{P}_{\Omega}$. Specifically, for an $N\times T$ matrix $R$,
The complement of the support set $\Omega$ is denoted by $\Omega^{c}$. Let the linear space of matrices supported on it be $\Omega^{\perp}$ and the projection onto the space be $\mathcal{P}_{\Omega^{\perp}}$, defined similarly as (ref).
Finally, we characterize the subgradient of $||L||_{*}$ and $||S||_{1}$ using these notations. The subgradient of $||L||_{*}$ evaluated at $L_{0}$ is equal to $UV^{*}+W$ where $\mathcal{P}_{\Phi}W=0$ and $||W||\leq 1$. Meanwhile, recall that $E$ denotes the sign of $S_{0}$, so the subgradient of $||S||_{1}$ at $S_{0}$ is equal to $E+F$, where $\mathcal{P}_{\Omega}F=0$ and $||F||_{\infty}\leq 1$.
Now we are ready to define the dual certificate $W$ as any $N\times T$ matrix satisfying the following conditions:
Note this definition is very similar to equation (6) in Ganesh et al.\enspace(2010). The only important difference occurs on the right hand side of the third inequality.
Now we present the construction of $W$ that is similar to Ganesh et al.\enspace(2010) with necessary modifications to accommodate $\delta\to 0$.
Let $W=W_{L}+W_{S}$:
The following lemmas then provide sufficient conditions for $W=W_{L}+W_{S}$ to satisfy (ref). The proof of Lemma (ref) is omitted because it is immediately implied by Theorem 2.6 in Cand\`es et al.\enspace(2011), stated in the Appendix. The proof of Lemma (ref) follow Lemma 3 and 4 in Ganesh et al.\enspace(2010) closely, but are tailored in a way to allow $\delta\to 0$. For completeness, they are contained in the Appendix.
In this section, using the dual certificate $W$ and the complementarity of the support sets of $S_{0}$ and $D_{0}$, we derive the optimality condition for the minimization problem in (ref) and as it turns out, the condition induces an upper bound for (ref) that facilitates the analysis of consistency.
For any feasible solution $L$ to (ref), let $Z\equiv Y-L$. Let $H\equiv L-L_{0}$ so $Z=Z_{0}-H$. By the Bernoulli device, we have the following properties: i) if we consider the difference $H_{S}\equiv Z-S_{0}$, then by construction $H_{S}=D_{0}-H$, and ii) $\mathcal{P}_{\Omega^{\perp}}D_{0}=D_{0}$. Utilizing them, we have the following lemma.
From Lemma (ref), we obtain bounds for $||\mathcal{P}_{\Phi^{\perp}}H||_{*}$ and $||\mathcal{P}_{\Omega^{\perp}}H||_{1}$. The next lemma bounds the corresponding norms of their complements.
Given Lemmas (ref) and (ref), now we are ready to state and prove the main result in the paper.
In this section, we present two simulation experiments to illustrate the effectiveness of the BPCP estimator. In the first experiment, we generate random $L_{0}$ and $Z_{0}$ with $r=1,3,5$. In each case, we examine the performance of the estimator for both Gaussian and Cauchy error matrices by comparing $\frac{1}{NT}||\hat{L}-L_{0}||_{F}^{2}$ as well as the relative estimation error $||\hat{L}-L_{0}||_{F}^{2}/||L_{0}||_{F}^{2}$ in each case as we gradually increase the sample size. In the second experiment, we fix a picture as $L_{0}$, and superimpose it with Gaussian and Cauchy white noise. The purpose of this example is to visually show how the estimator performs when the error is continuously fat-tailed distributed.
To implement the estimator, note the minimization problem in (ref) is equivalent as
We first set $\alpha$ to be a large number and solve the problem without the inequality constraint, then verify whether the solution satisfies the inequality. For the first step, we adopt the Augmented Lagrangian Multiplier algorithm (ALM) studied in Lin et al.\enspace(2010) and Yuan and Yang (2013). The algorithm solves the following problem
where $\Lambda$ is the Lagrangian multiplier for the equality constraint and the last term is a penalty for deviating from the constraint. The algorithm solves the minimization problem iteratively and terminates if both $||Y-L_{k}-Z_{k}||_{F}$ and $\nu ||Z_{k}-Z_{k-1}||_{F}$ are small enough, where the subscript $k$ denotes the $k$th iteration. We set the stopping criteria to be $10^{-7}||Y||_{F}$ and $10^{-5}$ respectively. Following Yuan and Yang (2013) and Cand\`es et al.\enspace(2011), we set $\nu=\frac{NT}{4||Y||_{1}}$. Finally, $\lambda$ is set to be $0.7(\frac{\log(N\land T)}{NT})^{1/3}$ in the first experiment and $0.5(\frac{\log(N\land T)}{NT})^{1/3}$ in the second; both of them satisfy the conditions in the theory when $\mu$ is of the order $\log(N)$ and $N\asymp T$.
In this experiment, we set $N=T=200,300,400,500$ and $r=1,3,5$. In each case, we draw an $N\times r$ and an $r\times T$ matrix from $N(0,1)$ and their product is $L_{0}$. For Gaussian error, entries of $Z_{0}$ are independently drawn from $N(0,1)$ while for Cauchy error, they are drawn from the standard Cauchy distribution.
Table (ref) shows the value $\frac{1}{NT}||\hat{L}-L_{0}||_{F}^{2}$ in each case. It decreases as $N$ and $T$ increase. Also, consistent with Theorem (ref), fixing $N$ and $T$, the estimation error increases as the rank increases. However, the magnitude of the increase is smaller than indicated by Theorem (ref), where the upper bound is proportional to $r^{3}$. This suggests our bound may be further improved and deriving the lower bound may also be useful. Another observation is that the estimation error under Cauchy errors is systematically bigger than under Gaussian errors. To explain this phenomenon, note that the standard Cauchy density is lower than the standard normal density around $0$, so for any given $\delta$ that is small enough, $-\underline{\gamma}_{it}$ and $\bar{\gamma}_{it}$ introduced in Section (ref) are bigger for the Cauchy error, resulting in bigger $||D_{0}||_{1}$. This difference is negligible when $N,T\to\infty$ but leads to finite sample differences in the error bound.
Table (ref) shows the relative estimation error $\frac{||\hat{L}-L_{0}||_{F}^{2}}{||L||_{F}^{2}}$. We can see that in every column, this quantity is also decreasing. Under the same error distribution, the differences across ranks are now significantly smaller. This is because as the rank increases, $||L_{0}||_{F}^{2}$ also increases by construction. Specifically, since $L_{0}$ is the product of two independent matrices with i.i.d. standard normal entries, $||L_{0}||_{F}^{2}$ is of the order of $NTr$, so dividing it mitigates the effect of $r$.
In this experiment, we superimpose white noises drawn from the standard normal or the standard Cauchy distribution on a picture. The picture we use is directly downloaded from a built-in example in MATLAB (the file name is eight.tif). The picture has resolution $242\times 308$. We stack all columns into a long vector and duplicate it for 199 times, obtaining a $74536\times 200$ matrix $L_{0}$. By construction, $L_{0}$ has rank $1$ because all columns are equal. Then we draw a $74536\times 200$ matrix with i.i.d. entries from either of the two distributions as $Z_{0}$ and add it to $L_{0}$. This procedure simulates 200 frames from a video of the static picture interfered by white noise.
Figure (ref) and Figure (ref) display the results. In each figure, the northwest (NW) is the original picture. There are four coins in it, two heads and two tails. The southwest (SW) shows one of the 200 frames after the the picture is superimposed with Gaussian noises (Figure (ref)) or Cauchy noises (Figure (ref)). It can be seen that the details of the coins are no longer unrecognizable. The northeast (SE) quadrant shows the recovered picture. Differences in this picture between Figure (ref) and Figure (ref) are hardly to be seen, except the background in Figure (ref) is slightly darker. The southeast quadrant shows the same frame of residuals, in which we cannot see contours of the coins, indicating it contains very few information about the original picture.
This paper provides sufficient conditions under which the Bounded Principal Component Pursuit (BPCP) is consistent for the random or deterministic low-rank component, whose entries can go to infinity with $N$ and $T$. Built on the prior work on Principal Component Pursuit (PCP), consistency is shown by constructing a Bernoulli device which induces sufficiently tight bounds for the estimation error whose order does not depend on the moments of the error $Z_{0}$; the random error only needs to have zero median, analogous to the LAD estimator for linear models. The results indicate the estimator is highly robust to large outliers.
There are a few aspects of the estimator that remains to be studied. This paper only shows consistency in the Frobenius norm. It would be interesting to see whether consistency holds componentwisely or in the $\ell_{\infty}$ norm. These aspects will be studied in future work.