Low Rank Matrix Recovery with Simultaneous Presence of Outliers and Sparse Corruption

TL;DR

This paper introduces a robust PCA algorithm capable of simultaneously handling element-wise sparse corruptions and outlier columns in data matrices, with a scalable randomized implementation and theoretical guarantees.

Contribution

It proposes a novel robust PCA method that addresses both sparse and column-wise corruptions simultaneously, extending beyond existing algorithms.

Findings

The algorithm effectively distinguishes inliers from outliers using sparse approximation.

The approach is scalable with a randomized design.

Theoretical analysis guarantees accurate sparse representation despite corruptions.

Abstract

We study a data model in which the data matrix D can be expressed as D = L + S + C, where L is a low rank matrix, S an element-wise sparse matrix and C a matrix whose non-zero columns are outlying data points. To date, robust PCA algorithms have solely considered models with either S or C, but not both. As such, existing algorithms cannot account for simultaneous element-wise and column-wise corruptions. In this paper, a new robust PCA algorithm that is robust to simultaneous types of corruption is proposed. Our approach hinges on the sparse approximation of a sparsely corrupted column so that the sparse expansion of a column with respect to the other data points is used to distinguish a sparsely corrupted inlier column from an outlying data point. We also develop a randomized design which provides a scalable implementation of the proposed approach. The core idea of sparse approximation is analyzed analytically where we show that the underlying ell_1-norm minimization can obtain the representation of an inlier in presence of sparse corruptions.

Figures

Tables

Table 1. TABLE I: Recovery error in the sparse component using algorithm ( 4 ).

$r$	2	5	10	15
${\Vert {𝐒}^{{}^{\prime }}-{\widehat{𝐒}}^{{}^{\prime }}\Vert }_F/{\Vert {𝐒}^{{}^{\prime }}\Vert }_F$	0.04	0.26	0.51	0.65

Equations

$$\mathbf{D}=\mathbf{L}+\mathbf{S}+\mathbf{C},$$

$$\mathbf{D}=\mathbf{L}+\mathbf{S}.$$

$$\begin{array}{rlrl}& \underset{\dot{\mathbf{L}},\dot{\mathbf{S}}}{\min }& & \Vert \dot{\mathbf{L}}{\Vert }_{\ast }+\lambda \Vert \dot{\mathbf{S}}{\Vert }_1\\ & \text{subject to}& & \dot{\mathbf{L}}+\dot{\mathbf{S}}=\mathbf{D}.\end{array}$$

$$\mathbf{D}=\mathbf{L}+\mathbf{C}.$$

$$\mathbf{D}=\mathbf{L}+\mathbf{C}+\mathbf{S}.$$

$$\text{Log-Recovery Error}={\log }_{10}\left(\Vert \mathbf{U}-\hat{\mathbf{U}}{\hat{\mathbf{U}}}^T\mathbf{U}{\Vert }_F/\Vert \mathbf{U}{\Vert }_F\right),$$

$$\begin{array}{rlrl}& \underset{\dot{\mathbf{L}},\dot{\mathbf{S}}}{\min }& & \Vert \dot{\mathbf{L}}{\Vert }_{\ast }+\lambda \Vert \dot{\mathbf{S}}{\Vert }_1+\gamma \Vert \dot{\mathbf{C}}{\Vert }_{1,2}\\ & \text{subject to}& & \dot{\mathbf{L}}+\dot{\mathbf{S}}+\dot{\mathbf{C}}=\mathbf{D},\end{array}$$

$$\underset{\mathbf{z}}{\min }\Vert {\mathbf{A}}^{-i}\mathbf{z}-{\mathbf{a}}^i{\Vert }_0,$$

$$\underset{\mathbf{z}}{\min }\Vert {\mathbf{S}}^{-i}\mathbf{z}-{\mathbf{s}}^i{\Vert }_1\text{s.t.}{\mathbf{B}}^{-i}\mathbf{z}={\mathbf{b}}^i.$$

$$\begin{array}{rl}& {\mathcal{L}}_S=\{k\in [N_1]:{\mathbf{s}}_k=0\},\\ & {\mathcal{L}}_{{\mathbf{z}}_o^{\ast }}=\{k\in [N_1]:({\mathbf{s}}_k{)}^T{\mathbf{z}}_o^{\ast }=0\}.\end{array}$$

$$\xi =\sqrt{\frac{2}{\pi }}\frac{N_1-n_s}{\sqrt{r_b}}-2\sqrt{N_1-n_s}-t_1\sqrt{\frac{N_1-n_s}{r_b-1}},$$

$$\begin{array}{rl}& \frac{1}{2}\xi >n_s^{{}^{\prime }}+\sum _{i\in {\mathcal{L}}_{{\mathbf{z}}_o^{\ast }}}\Vert {\mathbf{s}}_i\Vert +\sqrt{n_s-n_s^{{}^{\prime }}}{t}_2\text{and}\\ & \frac{\Vert {\mathbf{e}}_i^T{\mathbf{P}}_b\Vert }{2\Vert {\mathbf{e}}_i^T{\mathbf{R}}_b\Vert }\xi >\sum _{i\in {\mathcal{L}}_{{\mathbf{z}}_o^{\ast }}}\Vert {\mathbf{s}}_i\Vert +\sqrt{n_s-n_s^{{}^{\prime }}}{t}_2,\end{array}$$

$$\underset{\mathbf{z}}{\min }\Vert {\mathbf{A}}^{-i}\mathbf{z}-{\mathbf{a}}^i{\Vert }_1$$

$$\mathbb{E}\left\{\frac{\Vert {\mathbf{e}}_i^T{\mathbf{P}}_b\Vert }{\Vert {\mathbf{e}}_i^T{\mathbf{R}}_b\Vert }\right\}=\frac{n-r_b}{r_b}.$$

$$\begin{array}{r}\underset{\mathbf{z}}{\min }\Vert \mathbf{Dz}{\Vert }_1+\lambda \Vert \mathbf{z}{\Vert }_1\text{s. t.}{\mathbf{z}}^T{\mathbf{e}}_i=1\end{array}$$

$$\begin{array}{rlrl}& \underset{{\dot{\mathbf{L}}}_m,{\dot{\mathbf{S}}}_m}{\min }& & \frac{1}{\sqrt{N_1}}\Vert {\dot{\mathbf{S}}}_m{\Vert }_1+\Vert {\dot{\mathbf{L}}}_m{\Vert }_{\ast }\\ & \text{subject to}& & {\dot{\mathbf{L}}}_m+{\dot{\mathbf{S}}}_m=\mathbf{M}\end{array}$$

$$\mathbf{D}=\mathbf{UQ}+\mathbf{S}+\mathbf{C},$$

$${\mathbf{D}}_{{\varphi }_1}={\mathbf{L}}_{{\varphi }_1}+{\mathbf{S}}_{{\varphi }_1}+{\mathbf{C}}_{{\varphi }_1},$$

$${\mathbf{d}}^i=\mathbf{U}{\mathbf{q}}^i+{\mathbf{s}}^i.$$

$$\underset{\dot{\mathbf{q}}}{\min }\Vert {\mathbf{\Phi }}_2^T{\mathbf{d}}^i-{\mathbf{\Phi }}_2^T\mathbf{U}\dot{\mathbf{q}}{\Vert }_1,$$

$${\mathbf{\Phi }}_2^T{\mathbf{d}}^i-{\mathbf{\Phi }}_2^T\mathbf{U}{\mathbf{q}}^{\ast }={\mathbf{\Phi }}_2^T{\mathbf{s}}^i.$$

$$\underset{\mathbf{z}}{\min }\Vert {\mathbf{D}}_{{\varphi }_1}\mathbf{z}{\Vert }_1+\lambda \Vert \mathbf{z}{\Vert }_1\text{s.t.}{\mathbf{z}}^T{\mathbf{e}}_i=1,$$

$$\begin{array}{rlrl}& \underset{{\dot{\mathbf{L}}}_{{\varphi }_1},{\dot{\mathbf{S}}}_{{\varphi }_1}}{\min }& & \frac{1}{\sqrt{N_1}}\Vert {\dot{\mathbf{S}}}_{{\varphi }_1}{\Vert }_1+\Vert {\dot{\mathbf{L}}}_{{\varphi }_1}{\Vert }_{\ast }\\ & \text{subject to}& & {\dot{\mathbf{L}}}_{{\varphi }_1}+{\dot{\mathbf{S}}}_{{\varphi }_1}={\mathbf{M}}_{{\varphi }_1},\end{array}$$

$$\underset{\dot{\mathbf{Q}}}{\min }\Vert {\mathbf{D}}_{{\varphi }_2}-{\mathbf{\Phi }}_2^T\hat{\mathbf{U}}\dot{\mathbf{Q}}{\Vert }_1.$$

$$\Vert (\mathbf{I}-\mathbf{U}{\mathbf{U}}^T)\hat{\mathbf{U}}{\Vert }_F\le 10^{-3}.$$

$$\underset{\mathbf{z}}{\min }\Vert \mathbf{Dz}{\Vert }_0\text{subject to}\Vert \mathbf{z}\Vert =1$$

$$\underset{\mathbf{z}}{\min }\Vert \mathbf{Sz}{\Vert }_0\text{subject to}\mathbf{Bz}=0\text{and}\Vert \mathbf{z}\Vert =1.$$

$$\underset{\mathbf{z}}{\min }\Vert \mathbf{Dz}{\Vert }_1\text{subject to}{\mathbf{z}}^T\mathbf{v}=1,$$

$$[{\mathbf{z}}^{\ast }(1:i-1),-1,{\mathbf{z}}^{\ast }(i:n-1){]}^T,$$

$$\underset{\mathbf{z}}{\min }\Vert \mathbf{Sz}{\Vert }_1\text{subject to}\mathbf{Bz}=0\text{and}{\mathbf{z}}^T\mathbf{v}=1.$$

Full text

[TABLE]

where $\mathbf{L}$ is a low rank matrix, $\mathbf{S}$ an element-wise sparse matrix and $\mathbf{C}$ a matrix whose non-zero columns are outlying data points. To date, robust PCA algorithms have solely considered models with either $\mathbf{S}$ or $\mathbf{C}$, but not both. As such, existing algorithms cannot account for simultaneous element-wise and column-wise corruptions. In this paper, a new robust PCA algorithm that is robust to simultaneous types of corruption is proposed. Our approach hinges on the sparse approximation of a sparsely corrupted column so that the sparse expansion of a column with respect to the other data points is used to distinguish a sparsely corrupted inlier column from an outlying data point. We also develop a randomized design which provides a scalable implementation of the proposed approach. The core idea of sparse approximation is analyzed analytically where we show that the underlying $\ell_{1}$-norm minimization can obtain the representation of an inlier in presence of sparse corruptions.

I Introduction

Standard tools such as Principal Component Analysis (PCA) has been routinely used to reduce dimensionality by finding linear projections of high-dimensional data into lower dimensional subspaces. The basic idea is to project the data along the directions where it is most spread out so that the residual information loss is minimized. This has been the basis for much progress in a broad range of data analysis problems, including problems in computer vision, communications, image processing, machine learning and bioinformatics [1, 2, 3, 4, 5].

PCA is notoriously sensitive to outliers, which prompted substantial effort in developing robust algorithms that are not unduly affected by outliers. Two distinct robust PCA problems were considered in prior work depending on the underlying data corruption model, namely, the low rank plus sparse matrix decomposition [6, 7] and the outlier detection problem [8].

I-A Low rank plus sparse matrix decomposition

In this problem, the data matrix $\mathbf{D}$ is a superposition of a low rank matrix $\mathbf{L}$ representing the low dimensional subspace, and a sparse component $\mathbf{S}$ with arbitrary support, whose entries can have arbitrarily large magnitude modeling element-wise data corruption [7, 6, 9, 10, 11, 12], i.e.,

[TABLE]

For instance, [7] assumes a Bernoulli model for the support of $\mathbf{S}$ in which each element of $\mathbf{S}$ is non-zero with a certain small probability. Given the arbitrary support, all the columns/rows of $\mathbf{L}$ may be affected by the outliers. The cutting-edge Principal Component Pursuit (PCP) approach developed in [6] and [7] directly decomposes $\mathbf{D}$ into its low rank and sparse components by solving a convex program that minimizes a weighted combination of the nuclear norm $\|\dot{\mathbf{L}}\|$(sum of singular values), and the $\ell_{1}$-norm $\|\dot{\mathbf{S}}\|_{1}$,

[TABLE]

If the column and row spaces of $\mathbf{L}$ are sufficiently incoherent and the non-zero elements of $\mathbf{S}$ sufficiently diffused, (2) can provably recover the exact low rank and sparse components [6, 7].

I-B Outlier detection

In the outlier detection problem, outliers only affect a portion of the columns of $\mathbf{L}$, i.e., corruption is column-wise. The given data is modeled as

[TABLE]

A set of the columns of the outlier matrix $\mathbf{C}$, the so-called outliers, are non-zero and they do not lie in the Column Space (CS) of $\mathbf{L}$. In this problem, it is required to retrieve the CS of $\mathbf{L}$ or locate the outlying columns.

Many approaches were developed to address this problem, including[13, 14, 15, 8, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 18]. In [8], it is assumed that $\mathbf{C}$ is a column sparse matrix (i.e., only few data points are outliers) and a matrix decomposition algorithm is proposed to decompose the data into low rank and column-sparse components. In [18], the $\ell_{2}$-norm in the PCA optimization problem is replaced with an $\ell_{1,2}$-norm to grant robustness against outliers. In [26], we leveraged the low mutual coherence between the outlier data points and the other data points to set them apart from the inliers. An alternative approach relies on the observation that small subsets of the columns of $\mathbf{L}$ are linearly dependent (since $\mathbf{L}$ lies in a low-dimensional subspace) while small subsets of the outlier columns are not given that outliers do not typically follow low-dimensional structures. Several algorithms exploit this feature to locate the outliers [23, 15].

I-C Notation

Capital and small letters are used to denote matrices and vectors, respectively. For a matrix $\mathbf{A}$, $\mathbf{a}_{i}$ is the $i^{\text{th}}$ row of $\mathbf{A}$, $\mathbf{a}^{i}$ its $i^{\text{th}}$ column, and $\mathcal{N}_{a}$ its null space, i.e., $\mathcal{N}_{a}$ is the complement of the row space of $\mathbf{A}$. $\|\mathbf{A}\|$ denotes its spectral norm, $\|\mathbf{A}\|_{*}$ its nuclear norm which is the sum of the singular values, and $\|\mathbf{A}\|_{1}$ its $\ell_{1}$-norm given by \|\mathbf{A}\|_{1}=\sum\limits_{i,j}\big{|}\mathbf{A}(i,j)\big{|}. In addition, the matrix $\mathbf{A}^{-i}$ denotes the matrix $\mathbf{A}$ with its $i^{\text{th}}$ column removed. In an $N$-dimensional space, $\mathbf{e}_{i}$ is the $i^{\text{th}}$ vector of the standard basis. For a given vector $\mathbf{a}$, $\|\mathbf{a}\|_{p}$ denotes its $\ell_{p}$-norm. The element-wise functions $\operatorname{sgn}(.)$ and $|\>.\>|$ are the sign and absolute value functions, respectively.

I-D The identifiability problem

In this paper, we consider a generalized data model which incorporates simultaneous element-wise and column-wise data corruption. In other words, the given data matrix $\mathbf{D}\in\mathbb{R}^{N_{1}\times N_{2}}$ can be expressed as $\mathbf{D}=\mathbf{L}+\mathbf{S}+\mathbf{C}$, where $\mathbf{L}$ is the low rank matrix, $\mathbf{C}$ contains the outliers, and $\mathbf{S}$ an element-wise sparse matrix with a arbitrary support. Without loss of generality, we assume that the columns of $\mathbf{L}$ and $\mathbf{S}$ corresponding to the non-zero columns of $\mathbf{C}$ are equal to zero. We seek a robust PCA algorithm that can exactly decompose the given data matrix $\mathbf{D}$ into $\mathbf{L}$, $\mathbf{S}$ and $\mathbf{C}$. Without any further assumptions, this decomposition problem is clearly ill-posed. Indeed, there are many scenarios where a unique decomposition of $\mathbf{D}$ may not exist. For instance, the low rank matrix can be element-wise sparse, or the non-zero columns of $\mathbf{C}$ can be sparse, or the sparse matrix can be low rank. In the following, we briefly discuss various identifiability issues.

1. Distinguishing $\mathbf{L}$ from $\mathbf{S}$: The identifiability of the low rank plus sparse decomposition problem [7, 6] in which $\mathbf{D}=\mathbf{L}+\mathbf{S}$ was studied in [6]. This problem was shown to admit a unique decomposition as long as the column and row spaces of $\mathbf{L}$ are sufficiently incoherent with the standard basis and the non-zero elements of $\mathbf{S}$ are sufficiently diffused (i.e., not concentrated in few columns/rows) These conditions are intuitive in that they essentially require the low rank matrix to be non-sparse and the sparse matrix not to be of low rank.

2. Distinguishing outliers from inliers: Consider the outlier detection problem in which $\mathbf{D}=\mathbf{L}+\mathbf{C}$. Much research was devoted to study different versions of this problem, and various requirements on the distributions of the inliers and outliers were provided to warrant successful detection of outliers. The authors in [8] considered a scenario where $\mathbf{C}$ is column-sparse, i.e., only few data columns are actually outliers, and established guarantees for unique decomposition when the rank of $\mathbf{L}$ and the number of non-zero columns of $\mathbf{C}$ are sufficiently small. The approach presented in [23] does not necessitate column sparsity but requires that small sets of outliers be linearly independent. Under this assumption on the distribution of the outliers, exact decomposition can be guaranteed even if a remarkable portion of the data columns are outliers. In this paper, we make the same assumption about the distribution of the outliers, namely, we assume that an outlier cannot be obtained as a linear combination of few other outliers.

3. Distinguishing a sparse matrix from an outlier matrix: Suppose $\mathbf{D}=\mathbf{C}+\mathbf{S}$ and assume that the columns of $\mathbf{S}$ corresponding to the non-zero columns of $\mathbf{C}$ are equal to zero. Thus, if the columns of $\mathbf{S}$ are sufficiently sparse and the non-zero columns of $\mathbf{C}$ are sufficiently dense, one should be able to locate the outlying columns by examining the sparsity of the columns of $\mathbf{D}$. For example, suppose the support of $\mathbf{S}$ follows the Bernoulli model with parameter $\rho$ and that the non-zero elements of $\mathbf{C}$ are sampled from a zero mean normal distribution. If $N_{1}$ is sufficiently large, the fraction of non-zero elements of a non-zero column of $\mathbf{S}$ concentrates around $\rho$ while all the elements of a non-zero column of $\mathbf{C}$ are non-zero with very high probability.

I-E Data model

To the best of our knowledge, this is the first work to account for the simultaneous presence of both sources of corruption. In the numerical examples presented in this paper, we utilize the following data model.

Data Model 1**.**

*The given data matrix follows the following model.

1. The data matrix $\mathbf{D}\in\mathbb{R}^{N_{1}\times N_{2}}$ can be expressed as*

[TABLE]

*2. $\operatorname{rank}(\mathbf{L})=r$.

3. Matrix $\mathbf{C}$ has $K$ non-zero columns. Define $\{\mathbf{g}_{i}\}_{i=1}^{K}$ as the non-zero columns of $\mathbf{C}$. The vectors $\{{\mathbf{g}_{i}}/{\|\mathbf{g}_{i}\|}\}_{i=1}^{K}$ are i.i.d. random vectors uniformly distributed on the unit sphere $\mathbb{S}^{N_{1}-1}$ in $\mathbb{R}^{N_{1}}$. Thus, a non-zero column of $\mathbf{C}$ does not lie in the CS of $\mathbf{L}$ with overwhelming probability.

4. The non-zero elements of $\mathbf{S}$ follow the Bernoulli model with parameter $\rho$, i.e., each element is non-zero independently with probability $\rho$.

5. Without loss of generality, it is assumed that the columns of $\mathbf{L}$ and $\mathbf{S}$ corresponding to the non-zero columns of $\mathbf{C}$, are equal to zero.*

Remark 1**.**

The uniform distribution of the outlying columns (the non-zero columns of $\mathbf{C}$) over the unit sphere is not a necessary requirement for the proposed methods. We have made this assumption in the data model to ensure the following requirements are satisfied with high probability:

• •

The non-zero columns of $\mathbf{C}$ do not lie in the CS of $\mathbf{L}$.

• •

The non-zero columns of $\mathbf{C}$ are not sparse vectors.

• •

A small subset of the non-zero columns of $\mathbf{C}$ is linearly independent.

Similarly, the Bernoulli distribution of the non-zero elements of the sparse matrix $\mathbf{S}$ is not a necessary requirement. This assumption is used here to ensure that the support is not concentrated in some columns/rows. This is needed for $\mathbf{S}$ to be distinguishable from the outlier matrix $\mathbf{C}$ and for ensuring that the sparse matrix is not low rank with high probability.

The proposed data model is pertinent to many applications of machine learning and data analysis. Below, we provide two scenarios motivating the data model in (3).

I. Facial images with different illuminations were shown to lie in a low dimensional subspace [1]. Now a given dataset consists of some sparsely corrupted face images along with few images of random objects (e.g., buildings, cars, cities, etc). The images of the random objects cannot be modeled as face images with sparse corruption, which calls for means to recover the face images while being robust to the presence of random images.

II. A users rating matrix in recommender systems can be modeled as a low rank matrix owing to the similarity between people’s preferences for different products. To account for natural variability in user profiles, the low rank plus sparse matrix model can better model the data. However, profile-injection attacks, captured by the matrix $\mathbf{C}$, may introduce outliers in the user rating databases to promote or suppress certain products. The model (3) captures both element-wise and column-wise abnormal ratings.

I-F Motivating scenarios

The low rank plus sparse matrix decomposition algorithms – which only consider the presence of $\mathbf{S}$ – are not applicable to our generalized data model given that $\mathbf{C}$ is not necessarily a sparse matrix. Also, when $\mathbf{C}$ is column-sparse, it may well be low rank which violates the identifiability conditions of the PCP approach for the low rank plus sparse matrix decomposition [6]. As an illustrative example, assume $\mathbf{D}\in\mathbb{R}^{300\times 500}$ follows Data model 1 with $\rho=0.01$. We apply the decomposition method (2) to $\mathbf{D}$ and learn the CS of $\mathbf{L}$ from the obtained low rank component. Define

[TABLE]

where $\mathbf{U}$ is an orthonormal basis for the CS of $\mathbf{L}$ and $\hat{\mathbf{U}}$ is the learned basis. Fig. 1 shows the log recovery error versus $K$. Clearly, (2) cannot yield correct subspace recovery in the presence of outliers.

On the other hand, robust PCA algorithms that solely consider the column-wise corruption are bound to fail in the presence of the sparse corruption $\mathbf{S}$ since a crucial requirement of such algorithms is that a set of the columns of $\mathbf{D}$ lies in the CS of $\mathbf{L}$. However, in presence of the sparse corruption matrix $\mathbf{S}$, even the columns of $\mathbf{D}$ corresponding to the zero columns of $\mathbf{C}$ might not lie in the CS of $\mathbf{L}$. For instance, assume $\mathbf{D}\in\mathbb{R}^{100\times 400}$ follows Data model 1 with $K=200$ and $r=5$. Fig. 2 shows the log recovery error versus $\rho$. In this example, the robust PCA algorithm presented in [14] is utilized for subspace recovery. It is clear that the algorithm cannot yield correct subspace recovery for $\rho\geq 0.01$. The work of this paper is motivated by the preceding shortcomings of existing approaches.

On a first thought, one may be able to tackle the simultaneous presence of sparse corruption and outliers by solving

[TABLE]

where $\gamma$ and $\lambda$ are regularization parameters. This formulation combines the norms used in the algorithms in [8] and [6]. However, this method requires tuning two parameters but, more importantly, inherits the limitations of [8]. Specifically, the PCP approach in [8] requires the rank of $\mathbf{L}$ to be substantially smaller than the dimension of the data, and fails when there are too many outliers. Also, our experiments have shown that (4) does not yield an accurate decomposition of the data matrix. For illustration, consider $\mathbf{D}\in\mathbb{R}^{200\times 400}$ following Data model 1 with $\rho=0.01$, and $K=100$. Let $\mathbf{S}^{{}^{\prime}}$ be the columns of $\mathbf{S}$ indexed by the complement of the column support of $\mathbf{C}$ and $\hat{\mathbf{S}}^{{}^{\prime}}$ the corresponding sparse component recovered by (4). Table I shows the error in recovery of the sparse component ${\|\mathbf{S}^{{}^{\prime}}-{\hat{\mathbf{S}}}^{{}^{\prime}}\|_{F}}/{\|\mathbf{S}^{{}^{\prime}}\|_{F}}$ versus $r$. As shown, the convex program in (4) which combines (2) and [8] cannot yield an accurate decomposition knowing that in the absence of column outliers (i.e., when $K=0$) and setting $\gamma=0$, (2) does recover the sparse component with recovery error below 0.01 for all values of $r$ in Table I.

I-G Summary of contributions

In this paper, we develop a new robust PCA approach, dubbed the sparse approximation approach, which can account for both types of corruptions simultaneously. Below, we provide a summary of contributions.

• •

The sparse approximation approach: In this approach, we put forth an $\ell_{1}$-norm minimization formulation that allows us to find a sparse approximation for the columns of $\mathbf{D}$ using a sparse representation. This idea is used to locate the outlying columns of $\mathbf{D}$ which, once identified, reduces the primary problem to one of a low rank plus sparse matrix decomposition.

• •

We develop a new randomized design that provides a scalable implementation of the proposed method. In this design, the CS of $\mathbf{L}$ is learned using few randomly sampled data columns. Subsequently, the outliers are located using few randomly sampled rows of the data.

• •

We provide a mathematical analysis of the sparse approximation idea underlying our approach where we prove that the $\ell_{1}$-norm minimization can yield the linear representation of a sparsely corrupted inlier.

I-H Paper organization

The rest of the paper is organized as follows. In section II, the idea of sparse approximation is explained. Section III presents the proposed robust PCA method and Section IV exhibits the numerical experiments. The proofs of all the theoretical results are provided in the appendix along with additional theoretical investigations of the sparse approximation problem.

II Sparse approximation of sparsely corrupted data

Suppose the vector $\mathbf{h}$ lies in the CS of a matrix $\mathbf{H}$, i.e., $\mathbf{h}=\mathbf{H}\mathbf{b}$, where vector $\mathbf{b}$ is the linear representation of $\mathbf{h}$ with respect to the columns of $\mathbf{H}$. In a least-square sense, this representation can be obtained as the optimal point of $\underset{{\hat{\mathbf{b}}}}{\min}\>\>\|\mathbf{H}\hat{\mathbf{b}}-\mathbf{h}\|_{2}$. The main question that this section seeks to address is whether we can recover such a representation using a convex optimization formulation when both $\mathbf{H}$ and $\mathbf{h}$ are sparsely corrupted. In this section, we propose an $\ell_{1}$-norm minimization problem and prove that it can yield the underlying representation. We refer to this approach as “sparse approximation”.

II-A Theoretical result

The following definition presents the notion of sparse approximation.

Definition 1**.**

Suppose $\mathbf{A}\in\mathbb{R}^{N_{1}\times n}$ can be expressed as $\mathbf{A}=\mathbf{B}+\mathbf{S}$, where $\mathbf{S}$ is a sparse matrix, $n\ll N_{1}$, and $r_{b}<n$ ($r_{b}$ is the rank of $\mathbf{B}$). We say that $\mathbf{A}^{-i}\mathbf{z}$, where $\mathbf{z}\in\mathbb{R}^{n-1}$, is a sparse approximation of $\mathbf{a}^{i}$ if $\mathbf{b}^{i}={\mathbf{B}}^{-i}\mathbf{z}$. Thus, if $\mathbf{A}^{-i}\mathbf{z}$ is a sparse approximation of $\mathbf{a}^{i}$, then $\mathbf{A}^{-i}\mathbf{z}-\mathbf{a}^{i}={\mathbf{S}}^{-i}\mathbf{z}-\mathbf{s}^{i}$.

The reason that we refer to $\mathbf{A}^{-i}\mathbf{z}$ as the sparse approximation of $\mathbf{a}^{-i}$ is that $\mathbf{A}^{-i}\mathbf{z}-\mathbf{a}^{i}={\mathbf{S}}^{-i}\mathbf{z}-\mathbf{s}^{i}$ is a sparse vector if $\rho$ and $n$ are sufficiently small.

Assume that the CS of $\mathbf{B}$ does not include sparse vectors, i.e., the CS of $\mathbf{B}$ is not coherent with the standard basis. According to Definition 1, if $\rho$ and $n$ are small enough and $\mathbf{b}^{i}$ lies in the CS of ${\mathbf{B}}^{-i}$, then $\mathbf{A}^{-i}\mathbf{z}^{*}$ is a sparse approximation of $\mathbf{a}^{i}$, where $\mathbf{z}^{*}$ is the optimal point of

[TABLE]

because the span of $\mathbf{B}$ does not contain sparse vectors so the only way to obtain a sparse linear combination is to cancel out $\mathbf{b}^{i}$. Moreover, we could show that if the $\ell_{0}$-norm of (5) is relaxed to an $\ell_{1}$-norm, we will still able to obtain the sparse approximation. The following lemma establishes that the sparse approximation can be recovered by solving a convex $\ell_{1}$-norm minimization problem if $\mathbf{S}$ is sufficiently sparse. In order to obtain concise sufficient conditions, Lemma 1 assumes a randomized model for the distribution of the rows of $\mathbf{B}$ in its row space. In the appendix, we present deterministic sufficient conditions for a more general optimization problem.

Assumption 1**.**

The rows of the matrix $\mathbf{B}$ are i.i.d. random vectors uniformly distributed on the intersection of the row space of $\mathbf{B}$ and the unit sphere $\mathcal{S}^{n-1}$.

Before we state the lemma, we define $\mathbf{z}_{o}^{*}$ as the optimal point of the following oracle optimization problem

[TABLE]

In addition, define $n_{s}$ and $n_{s}^{{}^{\prime}}$ as the cardinalities of $\mathcal{L}_{S}^{c}$ and $\mathcal{L}_{\mathbf{z}_{o}^{*}}\cap\mathcal{L}_{S}^{c}$, respectively, where $\mathcal{L}_{S}$ and $\mathcal{L}_{\mathbf{z}_{o}^{*}}$ are defined as

[TABLE]

Thus, $n_{s}$ is the number of non-zero rows of $\mathbf{S}$ and $n_{s}^{{}^{\prime}}$ is the number of non-zero rows of $\mathbf{S}$ which are orthogonal to $\mathbf{z}_{o}^{*}$. If the support of $\mathbf{S}$ follows the random model, $n_{s}^{{}^{\prime}}$ is much smaller than $n_{s}$. In addition, if $n\ll N_{1}$ and $\rho$ is small enough, $n_{s}$ will be much smaller than $N_{1}$.

Lemma 1**.**

Suppose matrix $\mathbf{A}\in\mathbb{R}^{N_{1}\times n}$ is a full rank matrix which can be expressed as $\mathbf{A}=\mathbf{B}+\mathbf{S}$, where $\mathbf{B}$ follows Assumption 1. Define

[TABLE]

and define $\mathbf{R}_{b}$ and $\mathbf{P}_{b}$ as orthonormal bases for the row space and null space of $\mathbf{B}$, respectively. If

[TABLE]

then $\mathbf{z}_{o}^{*}$ (the optimal point of (6)) is the optimal point of

[TABLE]

with probability at least $1-\exp(-t_{1}^{2}/2)-\exp\left(-\frac{r_{b}}{2}(t_{2}^{2}-\log t_{2}^{2}-1)\right)$, for all $t_{1}\geq 0\>,\>t_{2}>1$.

Remark 2**.**

Suppose the support of $\mathbf{S}$ follows the Bernoulli model with parameter $\rho$. If $\rho$ and $n$ are small enough, $n_{s}\ll N_{1}$. Thus, the order of $\xi$ is roughly $\frac{N_{1}}{\sqrt{r_{b}}}$. Define $\kappa$ such that $\kappa\geq\underset{i,j}{\max}\>{\mathbf{S}(i,j)}$. The vector $\mathbf{z}_{o}^{*}$ cannot be simultaneously orthogonal to too many non-zero rows of $\mathbf{S}$. Thus, $n_{s}^{{}^{\prime}}$ is much smaller than $n_{s}$. Therefore, the order of the RHS of (8) is roughly \mathcal{O}\Big{(}\sqrt{n_{s}}+n_{s}^{{}^{\prime}}\kappa(\sqrt{\rho n}+1)\Big{)}. If we assume that the row space of $\mathbf{B}$ is a random $r_{b}$-dimensional subspace in the $n$-dimensional space, then

[TABLE]

Accordingly, the sufficient conditions in Lemma 1 amount to a requirement that $\frac{(n-r_{b})N_{1}}{r_{b}\sqrt{r_{b}}}$ is sufficiently larger than \Big{(}\sqrt{n_{s}}+n_{s}^{{}^{\prime}}\kappa(\sqrt{\rho n}+1)\Big{)}. In our problem, $r_{b}\ll N_{1}$ and $n_{s}\ll N_{1}$. Thus, the sufficient conditions of Lemma 1 are naturally satisfied.

III Proposed sparse approximation method

In this section, the sparse approximation (SA) method is presented. We also present a randomized design which can reduce the complexity of the proposed method from $\mathcal{O}(N_{2}^{3})$ to $\mathcal{O}(r^{2}\max(N_{1},N_{2}))$. The table of Algorithm 1 presents the proposed algorithm. The $\ell_{1}$-norm functions are utilized to enforce sparsity to both the representation vector $\mathbf{z}$ and the residual vector $\mathbf{D}\mathbf{z}$ [27, 28, 29, 30, 31]. The main idea is to first locate the non-zero columns of $\mathbf{C}$ to reduce the problem to that of a low rank plus sparse matrix decomposition. In order to identify the outliers, we attempt to find a sparse approximation for each data column $\mathbf{d}^{i},i=1,\ldots,N_{2}$, using a sparse linear combination of the columns of $\mathbf{D}^{-i}$. If for certain columns such an approximation cannot be found, we identify these columns as outliers.

The key idea underlying this approach is that the sparsity rate of $\mathbf{D}\mathbf{z}_{i}^{*}$ in (10) can be used to certify the identity of $\mathbf{d}^{i}$ being an outlier or a sparsely corrupted inlier. Before providing some insight, we consider an illustrative example in which $\mathbf{D}\in\mathbb{R}^{200\times 500}$ follows Data model 1, $r=5$, $\rho=0.01$ and the first 50 columns of $\mathbf{C}$ are non-zero ($K=50$). Fig. 3 shows $\mathbf{D}\>\mathbf{z}_{50}^{*}$ and $\mathbf{D}\>\mathbf{z}_{51}^{*}$. The outlying column is clearly distinguishable.

Insight for the SA method: To gain more insight, consider the scenario where the $i$-th column of $\mathbf{C}$ is zero, i.e., $\mathbf{c}^{i}=\mathbf{0}$, so that the $i$-th data column is a sparsely corrupted inlier. In this case, if the regularization parameter $\lambda$ is chosen appropriately, (10) can identify a sparse vector $\mathbf{z}_{i}^{*}$ (sparsity of $\mathbf{z}$ is promoted by the $\ell_{1}$-norm regularizer) such that $\mathbf{D}\mathbf{z}_{i}^{*}$ is also sparse. The $\ell_{1}$-norm functions forces (10) to put the non-zero values of $\mathbf{z}_{i}^{*}$ on the columns of $\mathbf{D}$ such that a linear combination of their low rank components cancel out the low rank component of $\mathbf{d}^{i}$ (i.e. they provide a sparse approximation for $\mathbf{d}_{i}$) and the linear combination of their sparse component yields a sparse vector. In other word, $\mathbf{L}\mathbf{z}_{i}^{*}=0$ and $\mathbf{S}\mathbf{z}_{i}^{*}$ is a sparse vector since it is a linear combination of few sparse vectors and the algorithm automatically puts the non-zero values of $\mathbf{z}_{i}^{*}$ on the columns such that $\mathbf{S}\mathbf{z}_{i}^{*}$ is roughly as sparse as possible.

On the other hand, if $\mathbf{c}^{i}\neq\mathbf{0}$, i.e., $\mathbf{c}^{i}$ is an outlying column, $\mathbf{D}\mathbf{z}_{i}^{*}$ is not likely to be sparse for a sparse $\mathbf{z}_{i}^{*}$ since small subsets of outlying columns are linearly independent, and an outlying column is unlikely to admit a sparse representation in the sparsely corrupted columns of $\mathbf{L}$, to say that linear combinations of few sparsely corrupted columns of $\mathbf{L}$ are unlikely to approximate an outlying column.

III-A Randomized implementation of the proposed method

Randomized techniques are utilized to reduce the sample and computational complexities of robust PCA algorithms [32, 33, 34, 35, 36, 37, 38, 39]. Algorithm 1 solves an $N_{2}\times N_{2}$ dimensional optimization problem to identify all the outlying columns. However, here we show that this problem can be simplified to a low-dimensional subspace learning problem. Let $\mathbf{L}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$ be the compact singular value decomposition (SVD) of $\mathbf{L}$, where $\mathbf{U}\in\mathbb{R}^{N_{1}\times r}$, $\mathbf{\Sigma}\in\mathbb{R}^{r\times r}$ and $\mathbf{V}\in\mathbb{R}^{N_{2}\times r}$. We can rewrite (3) as

[TABLE]

where $\mathbf{Q}=\mathbf{\Sigma}\mathbf{V}^{T}$. We call $\mathbf{Q}$ the representation matrix. The table of Algorithm 2 details the randomized implementation of the proposed SA method. In the randomized implementation, first the CS of $\mathbf{L}$ is obtained using a random subset of the data columns. Matrix $\mathbf{D}_{\phi_{1}}\in\mathbb{R}^{N_{1}\times m_{1}}$ consists of $m_{1}$ randomly sampled columns. The proposed outlier detection approach is applied to $\mathbf{D}_{\phi_{1}}$ to identify the outlying columns of $\mathbf{D}_{\phi_{1}}$. The matrix $\mathbf{D}_{\phi_{1}}$ can be expressed as

[TABLE]

where $\mathbf{L}_{\phi_{1}}$, $\mathbf{S}_{\phi_{1}}$ and $\mathbf{C}_{\phi_{1}}$ are the corresponding columns sampled from $\mathbf{L}$, $\mathbf{S}$ and $\mathbf{C}$, respectively. If $m_{1}$ is sufficiently large, $\mathbf{L}_{\phi_{1}}$ and $\mathbf{L}$ will have the same CS. Thus, if we remove the outlying columns of $\mathbf{D}_{\phi_{1}}$ and decompose the resulting matrix, $\mathbf{M}_{\phi_{1}}$, the obtained low rank component can yield the CS of $\mathbf{L}$.

Suppose the CS of $\mathbf{L}$ is learned correctly and assume that the $i$-th column of $\mathbf{C}$ is equal to zero. Thus, the $i$-th column of $\mathbf{D}$ can be represented as

[TABLE]

It was shown in [33] that $\mathbf{q}^{i}$ can be obtained as the optimal point of

[TABLE]

if $m_{2}$ is sufficiently large and some mild sufficient conditions are satisfied. Define $\mathbf{q}^{*}$ as the optimal point of (15). If $\mathbf{q}^{*}=\mathbf{q}^{i}$, then

[TABLE]

The vector $\mathbf{\Phi}_{2}^{T}\mathbf{s}^{i}$ consists of $m_{2}$ randomly sampled elements of the sparse vector $\mathbf{s}^{i}$. Accordingly, if $\mathbf{c}^{i}$ (the $i$-th column of $\mathbf{C}$) is equal to zero, then $\mathbf{\Phi}_{2}^{T}\mathbf{d}^{i}-\mathbf{\Phi}_{2}^{T}\mathbf{U}{\mathbf{q}}^{*}$ is a sparse vector. However, if $\mathbf{c}^{i}$ is not equal to zero and $m_{2}$ is sufficiently large, it will be highly unlikely that $\mathbf{\Phi}_{2}^{T}\mathbf{d}^{i}-\mathbf{\Phi}_{2}^{T}\mathbf{U}{\mathbf{q}}^{*}$ is a sparse vector since $\mathbf{\Phi}_{2}^{T}\mathbf{U}{\mathbf{q}}^{*}$ cannot cancel out the component of $\mathbf{\Phi}_{2}^{T}\mathbf{c}^{i}$ that does not lie in the CS of $\mathbf{\Phi}_{2}^{T}\mathbf{U}$.

Remark 3**.**

In the CS learning step, we identify the outlying columns via the sparsity of $\mathbf{D}_{\phi_{1}}\mathbf{z}_{i}^{*}$. If $\mathbf{z}_{i}^{*}$ lies in the null space of $\mathbf{L}_{\phi_{1}}$, then $\mathbf{z}_{i}^{*}$ has $c\>r$ non-zero elements where $c>1$. According to our investigations, if $\lambda$ is chosen appropriately, $c$ is a small number mostly smaller than 3. Thus, $\mathbf{D}\mathbf{z}_{i}^{*}$ has at most $2r\rho N_{1}$ non-zero elements. In practice, it is much smaller than $2r\rho N_{1}$ since the optimization searches for the most sparse linear combination. In step 3 of the algorithm, the outlying columns are located by examining the sparsity of the columns of $\mathbf{H}=\mathbf{D}_{\phi_{2}}-\mathbf{\Phi}_{2}^{T}\hat{\mathbf{U}}\hat{\mathbf{Q}}$. If the $i$-th column of $\mathbf{C}$ is equal to zero, $\hat{\mathbf{U}}$ and $\hat{\mathbf{Q}}$ are correctly recovered, and $m_{2}$ is sufficiently large, then the $i$-th column of $\mathbf{H}$ is a sparse vector with roughly $\rho m_{2}$ non-zero elements. Accordingly, if we set an appropriate threshold for the number of dominant non-zero elements, the outlying columns are correctly identified.

IV Numerical Simulations

In this section, we present a set of numerical experiments to study the performance of the proposed approach. First, we validate the idea of sparse approximation for outlier detection and study its requirements. Second, we provide a set of phase transition plots to demonstrate the requirements of the randomized implementation, i.e., the sufficient number of randomly sampled columns/rows. Finally, we study the sparse approximation approach for outlier detection with real world data.

IV-A The idea of outlier detection

Suppose the given data follows Data model 1 and $\mathbf{D}\in\mathbb{R}^{100\times 200}$. The first 20 columns of $\mathbf{C}$ are non-zero. The matrix $\mathbf{S}$ follows the Bernoulli model with $\rho=0.01$. The rank of $\mathbf{L}$ is equal to 5. We solve (10) with the constraint vector set equal to $\{\mathbf{e}_{i}\}_{i=1}^{N_{2}}$ and define $\{\mathbf{z}_{i}^{*}\}_{i=1}^{N_{2}}$ as the corresponding optimal point. Define $\mathbf{h}_{i}=|\mathbf{D}\mathbf{z}^{*}_{i}|/\max(|\mathbf{D}\mathbf{z}^{*}_{i}|)$. We also, define a vector $\mathbf{g}\in\mathbb{R}^{N_{2}\times 1}$ whose $i$-th entry is set equal to the number elements of $\mathbf{h}_{i}$ greater than $0.1$. Thus, the $i$-th element of $\mathbf{g}$ is the number of dominant non-zero elements of $\mathbf{D}\mathbf{z}_{i}^{*}$. Fig. 4 shows the elements of $\mathbf{g}/N_{2}$. As shown, the indices corresponding to the outlying columns are clearly distinguishable.

IV-B Phase transition

In the presented theoretical analysis, we have shown that if the rank of the low rank component is sufficiently small and the sparse component is sufficiently sparse, the $\ell_{1}$-norm optimization problem can yield the sparse approximation (Lemma 1 and Theorem 2). In this section, we assume the data $\mathbf{D}\in\mathbb{R}^{400\times 400}$ follows Data model 1 and study the phase transition of Algorithm 1 (which uses sparse approximation for outlier detection) in the 2D-plane of $r$ and $\rho$. The first 200 columns of the given data are outlying columns, i.e., $K=200$. We define a vector $\mathbf{h}_{i}=|\mathbf{D}\mathbf{z}^{*}_{i}|/\max(|\mathbf{D}\mathbf{z}^{*}_{i}|)$ as before corresponding to each data column, and we classify the $i$-th column as an outlier if more than 40 percent of the elements of $\mathbf{h}_{i}$ are greater than 0.1. Fig. 5 shows the phase transition in the plane of $r$ and $\rho$. For each pair $(r,\rho)$, we generate 10 random realizations. In this figure, white designates that all outliers are detected correctly and no inlier is misclassified as an outlier. One can observe that if $\rho<0.07$, we can correctly identify all the outliers even with $r=30$.

In practice, the proposed method can handle larger values of $r$ and higher sparsity levels (i.e., more non-zero elements) because the columns of the low rank matrix typically exhibit additional structures such as clustering structures [28, 40, 41]. In the simulation corresponding to Fig. 5, the low rank matrix was generated as $\mathbf{L}=\mathbf{U}\mathbf{Q}$, where the elements of $\mathbf{U}\in\mathbb{R}^{N_{1}\times r}$ and $\mathbf{Q}\in\mathbb{R}^{r\times N_{2}}$ are drawn from a zero mean normal distribution. Thus, the columns of $\mathbf{L}$ are distributed randomly in the CS of $\mathbf{L}$. Accordingly, the elements of $\mathbf{z}_{i}^{*}$ must have at least $r+1$ non-zero elements to yield the sparse approximation. However, if the columns of $\mathbf{L}$ lie in a union of, say, $n$ $r/n$-dimensional subspaces, $r/n+1$ non-zero elements can be sufficient. Thus, if the data exhibits a clustering structure, the algorithm can bear with higher rank and sparsity levels. As an example, assume $\mathbf{D}\in\mathbb{R}^{400\times 120}$ follows Data model 1 with $r=10$, $\rho=0.1$, and $K=40$ (the last 40 columns are outlying columns). Fig. 6 shows the sorted elements of $\mathbf{D}\mathbf{z}_{1}^{*}$. In the left plot, the columns of $\mathbf{L}$ lie in one 10-dimensional subspace, whereas in the right plot, the columns of $\mathbf{L}$ lie in a union of ten 1-dimensional subspaces. As shown, the proposed method yields a better output if the data admits a clustering structure.

IV-C Phase transition for the randomized implementation

In this section, the requirements of Algorithm 2 are studied. The data matrix $\mathbf{D}\in\mathbb{R}^{1000\times 1000}$ follows Data model 1. The phase transition shows the performance of Algorithm 2 as function of $m_{1}$ and $m_{2}$. A trial is considered successful if the rank of $\hat{\mathbf{U}}$ is equal to $r$, step 3.3 identifies the outlying columns correctly, and

[TABLE]

Fig. 7 shows the phase transition plots with different values of $r$, $\rho$ and $K/N_{2}$. One can see that the required values for $m_{1}$ and $m_{2}$ increase if the $r$ or $K/N_{2}$ increase. The required value for $m_{1}$ is approximately linear in $r\frac{N_{2}}{N_{2}-K}$ and the required value for $m_{2}$ is linear in $r$ [33, 39]. Fig. 8 Shows the phase transition for different dimensions of $\mathbf{D}$, where $r=5$, $\rho=0.02$ and $\frac{K}{N_{2}}=0.5$. Interestingly, the required values for $m_{1}$ and $m_{2}$ are nearly independent of the size of $\mathbf{D}$.

IV-D The proposed approach with real data

The Extended Yale Face Database [42] consists of face images of 38 human subjects under different illumination. We select 50 images of a human subject. Fig. 9 shows few samples of the face images. It has been observed that these images roughly follow the low rank plus sparse model, and the low rank plus sparse matrix decomposition algorithms were successfully applied to such images to remove shadows and specularities [7]. In addition, we add a sparse matrix with $\rho=0.01$ to the face images. We form a data matrix $\mathbf{D}\in\mathbb{R}^{32256\times 100}$ (32256 is the number of pixels per image) consisting of these 50 sparsely corrupted face images plus 50 randomly sampled images from the Caltech101 database [43] as outlying data points (50 images from random subjects). Fig. 10 shows a subset of the images sampled from the Caltech101 database. The first 50 columns are the face images and the last 50 columns are the random images. We have found that, for $1\leq i\leq 50$, the average number of non-zero elements of $\mathbf{D}\mathbf{z}^{*}_{i}/\max(|\mathbf{D}\mathbf{z}_{i}^{*}|)$ with absolute values greater than 0.1 is $0.05N_{2}$ with standard deviation $0.04N_{2}$. For $51\leq i\leq 100$, the average number of non-zero elements of $\mathbf{D}\mathbf{z}^{*}_{i}/\max(|\mathbf{D}\mathbf{z}_{i}^{*}|)$ with absolute value greater than 0.1 is $0.37N_{2}$ with standard deviation $0.14N_{2}$. Thus, the non-face images can be identified with a proper threshold on the number of non-zero elements of $\mathbf{D}\>\mathbf{z}^{*}_{i}$.

V Appendix

In this section, we study a more general theoretical problem, dubbed null space learning problem, of which the sparse approximation problem is a special case. Similar to the model used in Section II, assume that matrix $\mathbf{A}\in\mathbb{R}^{N_{1}\times n}$ is a full rank matrix that can be expressed as $\mathbf{A}=\mathbf{B}+\mathbf{S}$, where $\mathbf{S}$ is a sparse matrix and the columns of $\mathbf{B}$ are linearly dependent. If the CS of $\mathbf{B}$ is not coherent with the standard basis, we expect the optimal point of

[TABLE]

to lie in the null space of $\mathbf{B}$ (i.e., $\mathbf{B}\mathbf{z}=\mathbf{0}$) because the CS of $\mathbf{B}$ does not contain sparse vectors and the optimal point should cancel out the component corresponding to $\mathbf{B}$. Accordingly, the optimal point of (20) is equal to the optimal point of

[TABLE]

Therefore, the $\ell_{0}$-norm minimization problem can learn a direction in the null space of $\mathbf{B}$. In fact, the optimization problem (20) finds the most sparse vector in the CS111Interestingly, finding the most sparse vector in a linear subspace has bearing on, and has been effectively used in, other machine learning problems, including dictionary learning and spectral estimation [44, 45, 46]. of $\mathbf{D}$.

The optimization problem (20) is non-convex. We relax the cost function using an $\ell_{1}$-norm and replace the quadratic constraint with a linear constraint as

[TABLE]

where $\mathbf{v}\in\mathbb{R}^{n}$ is a fixed vector. We refer to (22) as the convex null space learning optimization problem. If we set the constraint vector $\mathbf{v}$ equal to $\mathbf{e}_{i}$, the null space learning optimization problem is equivalent to the sparse approximation optimization problem (9), hence the latter is a special case of the former. Note that in (9), $\mathbf{z}\in\mathbb{R}^{n-1}$ and in (22) $\mathbf{z}\in\mathbb{R}^{n\times 1}$, yet the equivalence stems from the fact that if $\mathbf{z}^{*}$ is the optimal point of (9), then the optimal point of (22) with $\mathbf{v}=\mathbf{e}_{i}$ will be equal to

[TABLE]

where $\mathbf{z}(i:j)$ denotes the elements of a vector $\mathbf{z}$ from index $i$ to $j$.

The following Theorem establishes sufficient conditions for the optimal point of the $\ell_{1}$-norm minimization problem (22) to lie in the null space of $\mathbf{B}$. Before we state the theorem, let us define $\mathbf{z}_{o}^{*}$ as the optimal point of

[TABLE]

The sets $\mathcal{L}_{S}$ and $\mathcal{L}_{\mathbf{z}_{o}^{*}}$ are defined similar to (7).

Theorem 2**.**

Suppose matrix $\mathbf{A}\in\mathbb{R}^{N_{1}\times n}$ is a full rank matrix that can be expressed as $\mathbf{A}=\mathbf{B}+\mathbf{S}$. Define $\mathbf{z}_{o}^{*}$ as the optimal point of (23), and define

[TABLE]

If

[TABLE]

where $\mathbf{b}_{i}$ and $\mathbf{s}_{i}$ are the $i^{\text{th}}$ row of $\mathbf{B}$ and $\mathbf{S}$, respectively, and the subspace $\mathcal{R}_{b}$ is the row space of $\mathbf{B}$, then $\mathbf{z}_{o}^{*}$ is the optimal point of (22).

The sufficient conditions (25) reveal some interesting properties which merit some intuitive explanation provided next. According to Theorem 2, the following are important factors to ensure that the optimal point of (22) lies in the null space of $\mathbf{B}$:

1. The CS of $\mathbf{B}$ should not be coherent with the standard basis: Recall that we assume that $n\ll N_{1}$. Thus, if $\mathbf{S}$ follows the Bernoulli model and $\rho$ and $n$ are sufficiently small, then $|\mathcal{L}_{S}^{c}|\ll|\mathcal{L}_{S}|$, in which case the LHS of (25) will approximate the permeance statistic[25] – a measure of how well the rows of $\mathbf{B}$ are distributed in the row space of $\mathbf{B}$. The permeance statistic increases if the rows are more uniformly distributed in the row space. But, if they are aligned along some specific directions, the permeance statistic tends to be smaller, wherefore the CS of $\mathbf{B}$ will be more coherent with the standard basis. This is in agreement with our initial intuition since linear combinations of the columns of $\mathbf{B}$ are more likely to form sparse vectors when $\mathbf{B}$ is highly coherent. In other words, the coherence of $\mathbf{B}$ would imply that $\mathbf{A}\mathbf{z}$ could be sparse even if $\mathbf{z}$ does not lie in the null space of $\mathbf{B}$.

2. The matrix $\mathbf{S}$ should be sufficiently sparse: Per the first inequality of (25), $\mathbf{S}$ should be sufficiently sparse otherwise the cardinality of $\mathcal{L}_{S}$ will not be sufficiently large. This requirement also confirms our initial intuition because the optimal point of (22) cannot yield a sparse linear combination even if it lies in the null space of $\mathbf{B}$ unless $\mathbf{S}$ is sparse.

3. The vector $\mathbf{v}$ should not be too incoherent with $\mathcal{N}_{b}$: Recalling that $\mathbf{P}_{b}$ and $\mathbf{R}_{b}$ are orthonormal bases for the null space and row space of $\mathbf{B}$, respectively, the factor

[TABLE]

on the LHS of the second inequality of (25) unveils that the vector $\mathbf{v}$ should be sufficiently coherent with $\mathcal{N}_{b}$ (the null space of $\mathbf{B}$) in order to ensure that the optimal point of (22) lies in $\mathcal{N}_{b}$. The intuition is that if $\mathbf{v}$ has a very small projection on $\mathcal{N}_{b}$ and the optimal point of (22) lies in $\mathcal{N}_{b}$, then the optimal point of (22) should have a large Euclidean norm to satisfy the linear constraint of (22). In that sense, points lying in $\mathcal{N}_{b}$ would be unlikely to attain the minimum of the objective function in (22). In addition, this coherency requirement implies that the optimal point of (22) is more likely to lie in $\mathcal{N}_{b}$ when the rank of $\mathbf{B}$ is smaller because the null space would have higher dimension, which makes $\mathbf{v}$ more likely to be coherent with $\mathcal{N}_{b}$.

V-A Proof of Lemma 1

Lemma 1 is a special case of Theorem 2. In order to prove Lemma 1, we make use of the following Lemmas from [47, 25] to lowerbound $\underset{\delta\in\mathcal{R}_{b}\atop\|\delta\|=1}{\inf}\>\>\sum_{i\in\mathcal{L}_{S}}|\mathbf{b}_{i}^{T}\delta|$ and upperbound $\|\alpha\|$ in (25).

Lemma 3**.**

(Lower-bound on the permeance statistic from [25]) Suppose that $\mathbf{g}_{1},...,\mathbf{g}_{n}$ are i.i.d. random vectors uniformly distributed on the unit sphere $\mathbb{S}^{r-1}$ in $\mathbb{R}^{r}$. When $r=1$,

[TABLE]

When $r\geq 2$, for all $t\geq 0$,

[TABLE]

with probability at least $1-exp(-t^{2}/2)\>$.

Lemma 4**.**

If $\mathbf{g}_{1},...,\mathbf{g}_{n}$ are i.i.d. random vectors uniformly distributed on the unit sphere $\mathbb{S}^{r-1}$ in $\mathbb{R}^{r}$, then

[TABLE]

for all $t>1$.

V-B Proof of Theorem 2

We want to show that

[TABLE]

Define $g(\delta)$ as

[TABLE]

Since (22) is a convex optimization problem, it suffices to check that $g(\delta)\geq 0$ for every sufficiently small non-zero perturbation $\delta$ such that

[TABLE]

The conditions on $\delta$ is to ensure that $\mathbf{z}_{o}^{*}-\mathbf{\delta}$ is a feasible point of (22). If $\mathbf{z}_{o}^{*}$ is the optimal point of (23), then the cost function of (23) is increased when we move from the optimal point along a feasible perturbation direction. Observe that $\mathbf{z}_{o}^{*}-\delta_{n}$ is a feasible point of (23), if and only if the perturbation $\delta_{n}$ satisfies

[TABLE]

where $\mathcal{N}_{b}$ is the null space of $\mathbf{B}$. Therefore, for any non-zero $\mathbf{\delta}_{n}$ which satisfies (32)

[TABLE]

When $\delta_{n}\to 0$, we can rewrite (33) as

[TABLE]

where the last identity follows from the Taylor expansion of the square root. Thus,

[TABLE]

has to be greater than zero for small $\delta_{n}$ which satisfies (32). Therefore,

[TABLE]

To simplify $g(\delta)$, we decompose $\delta$ into

[TABLE]

where $\delta_{1}\in\mathcal{R}_{b}$ and $\delta_{2}\in\mathcal{N}_{b}$, where $\mathcal{R}_{b}$ is the row space of $\mathbf{B}$. The vectors $\mathbf{z}_{o}^{*}$ and $\delta_{2}$ lie in $\mathcal{N}_{b}$. Therefore, we can expand $g(\delta)$ as

[TABLE]

Thus, it suffices to ensure that

[TABLE]

is non-negative since (39) is a lower bound on $g(\delta)$. In addition, as $\delta\to 0$,

[TABLE]

Therefore, according to (30), (40) and (39), it is enough to show that

[TABLE]

for every $\delta\neq 0$ which satisfies (31). Define

[TABLE]

Therefore, to show that $g(\delta)$ is non-negative, it suffices to ensure that

[TABLE]

For the first inequality of (43), it is enough to ensure that a lower bound on the LHS is greater than an upper bound on the RHS. Thus, it suffices to have

[TABLE]

Thus, the first inequality of (25) guarantees the first inequality of (43).

Observe that the second inequality of (43) is homogeneous in $\delta$ since

[TABLE]

We scale $\delta$ such that $\delta_{1}^{T}\mathbf{v}=-\delta_{2}^{T}\mathbf{v}=1$. To ensure that the second inequality of (43) is satisfied, it is enough to show that

[TABLE]

Let us decompose $\delta_{2}$ into

[TABLE]

where

[TABLE]

and $\mathbf{v}^{{}^{\prime}}$ is defined as

[TABLE]

where $\mathbf{P}_{b}$ was defined as an orthonormal basis for $\mathcal{N}_{b}$. For the second inequality, it is enough to show that the LHS of (46) is greater than

[TABLE]

According to the definition of $\mathbf{v}^{{}^{\prime}}$ and $\delta_{2p}$,

[TABLE]

Therefore, according to (36),

[TABLE]

Thus, it is enough to show that the LHS of (46) is greater than

[TABLE]

In addition, the LHS of (46) can be simplified as

[TABLE]

According to (46), (52) and (53), the second inequality of (25) guarantees that the second inequality of (43) is satisfied.