Robust principal component analysis with intra-block correlation

Abstract

This paper proposes a novel optimization program for solving the Robust Principle Component Analysis (RPCA) problem, which decomposes a data matrix into a conventional low-rank part plus a particular block-sparse residual. This kind of block-sparse residual often scattered in the source signal scope as contaminants, and often existed in many practical applications, such as an ordinary imaging system, a Hyper Spectral Imaging system, EEG and MEG, and types of physiological signals. Different from most currently existing approaches, the study perceived especially a highly spatial correlation among the inner structure of the neighbouring pixels in this contiguously block-sparse residual. The high intra-block correlation is then introduced as prior information to deal the governing optimization problem. In order to enhance the block-sparsity and maintain the local smoothness simultaneously, a localized low-rank promoting method is introduced with a theoretical guarantee. An efficient solving algorithm is designed accordingly with a convergence analysis by adopting the classical Alternating Direction Method of Multipliers (ADMM) framework. In addition to the theoretical model derivation, several synthetic simulations together with a real application on image denoising experiment have been conducted to validate the proposed model. As expected, the models outperforms significantly the compared state-of-the-arts.

Introduction

In the last decade, classical PCA [1] is the most widely used tool for dimensionality reduction, but is highly sensitive to grossly corrupted observations. Unfortunately, gross errors are now ubiquitous in practical applications such as face recognition [2], image inpainting [3], and video surveillance [4]. Robust Principle Component Analysis (RPCA) [5], [6], as an extension of PCA, have been proposed to seek the low-dimensional structure from highly corrupted measurements with arbitrary magnitude entries. This approach possessed of a polynomial-time algorithm receives much attention in many problems of Low-Rank and Sparse Decomposition (LRSD), including background subtraction [7], latent semantic indexing [8] and image/video denoising [9], etc.

Suppose we have a large data matrix $\mathbf{M}\in {\mathbb{R}}^{m\times n}$ which can generally be decomposed as $\mathbf{M}=\mathbf{L}+\mathbf{S},$ as shown in Fig. 1 (up row). The RPCA method aims to recover the low-rank matrix $\mathbf{L}\in {\mathbb{R}}^{m\times n}$ from the real observation matrix $\mathbf{M}\in {\mathbb{R}}^{m\times n}$ with the premise of modelling the corrupted noise parts with arbitrary distribution and magnitude as a sparse matrix $\mathbf{S}\in {\mathbb{R}}^{m\times n}$. Mathematically, it can be realized by solving the following optimization minimization [5], [10]:

$$\underset{\mathbf{L},\mathbf{S}}{\min }rank\left(\mathbf{L}\right)+\lambda {\parallel \mathbf{S}\parallel }_{0}s.t.\mathbf{M}=\mathbf{L}+\mathbf{S},$$

where rank( · ) works as a sparsity regularization of matrix singular values; ‖ · ‖₀ denotes the number of the nonzero elements, which is used as the regularization term to promote sparsity; and λ > 0 is denoted as the relatively balanced parameter between the two term. Unfortunately, minimizing the optimization function directly is an NP-hard problem, i.e., it is difficult to be solved within polynomial time. Hence, in order to deal with this issue, we consider a surrogate problem which can be formulated as computing the minimizer of

$$\underset{\mathbf{L},\mathbf{S}}{\min }{\parallel \mathbf{L}\parallel }_{*}+\lambda g\left(\mathbf{S}\right)s.t.\mathbf{M}=\mathbf{L}+\mathbf{S},$$

where ${\parallel \mathbf{L}\parallel }_{*}={\sum }_i{\delta }_i\left(\mathbf{L}\right)$ denotes the nuclear norm [11] and ${\delta }_i\left(\mathbf{L}\right)$ represents the ith singular value of $\mathbf{L}$ (sorted in decreasing order). In terms of the surrogate functions $g\left(\mathbf{S}\right)$ of the sparse penalty term ${\parallel \mathbf{S}\parallel }_0,$ the general option of classical RPCA is l₁-norm ${\parallel \mathbf{S}\parallel }_1={\sum }_{i=1}^m{\sum }_{j=1}^n\left|{\mathbf{S}}_{ij}\right|$ [5]. However, the l₁-norm is only a loose approximation of the l₀-norm, which will result in a suboptimal solution to the original l₀-minimization in reality. Several more efficient surrogate functions $g\left(\mathbf{S}\right)$ are hence proposed for extending the classical RPCA model, such as l_p-norm (0 < p < 1) [12], Smoothly Clipped Absolute Deviation (SCAD) [13], Exponential-Type Penalty (ETP) [14] and generalized Weighted-based Integral Convolution Penalty (g-WICP) [15], etc.

It is well known that the formula of l₀-norm treats each entry (pixel) in the sparse matrix independently, thus the classical RPCA method is extremely insensitive to column/row outliers. And we further find that the block-sparse structural information does exist abundantly in real circumstances. For example, the background (see Fig. 2(b)) of an image (see Fig. 2(a)) which randomly selected from Wallflower dataset [16] can be easily captured by the low-rank matrix, while the corresponding foreground object (see Fig. 2(c)) can be regarded as a block-sparse corruption for it is spatially contiguous and occupies a portion of the scene. Under the condition of block corruption, a l_2,0-norm which stands for the number of nonzero l₂-norm of the column was proposed to replace the original l₀-norm to manage the block corruption of the sparse matrix. Therefore the new optimization problem for decomposing an observation matrix into a low-rank part plus a block-sparse residual can be reformulated as:

$$\underset{\mathbf{L},\mathbf{S}}{\min }{\parallel \mathbf{L}\parallel }_{*}+\lambda {\parallel \mathbf{S}\parallel }_{2,0}s.t.\mathbf{M}=\mathbf{L}+\mathbf{S}.$$

An illustration for RPCA with a block-sparse residual is shown in the bottom row of Fig. 1 to contrast significantly to that with snow-sparse residual. For this block-sparse penalty term, the l_2,1-norm ${\parallel \mathbf{S}\parallel }_{2,1}={\sum }_{i=1}^m\sqrt{{\sum }_{j=1}^n|{\mathbf{S}}_{ij}^2}|$ was first adopted in [17] to approximate the l_2,0-norm for detecting outliers with column-wise sparsity. The l_2,p-norm has also been proved more exact than the l_2,1-norm when imposed on the reconstruction error to handle outlier pursuit [12]. Related works of low-rank and block-sparse matrix decomposition are currently following up with some researches. In [18], Tang and Nehorai proposed a RPCA-LBD detector to separate the primary components from the outliers, which conducted a good result when the sparse pattern involved in an entire column. With the potential for outlier detection applications, this approach then was wildly applied to many practical problems, such as salient motion detection [19] and imaging processing [20]. Moreover, a great progress in applying one alternative block-sparse RPCA to the topic of background modelling was achieved too. By Liu et al. [21] who proposed a novel approach, namely Low-rank and Structured sparse Decomposition (LSD), for foreground detection, and attained the algorithm successfully through the method of Augmented Lagrange Multiplier (ALM). The essential of such an issue is still regarded as the decomposition of an underlying image into a low-rank matrix and a block-sparse outlier matrix.

Although there were few previous works considered to include the prior information about the high-correlated structure in the outlier matrix. Apart from block sparsity, the images in real applications also often have the tendency to own extra useful structural compositions, which can be explored to improve the reconstruction performance. As shown in Fig. 3(a), it is obvious to find that there are several non-zero structural blocks in the actual foreground image (Fig. 3(b)). Upon this observation, a small non-zero region which owned such kind of structural block was randomly picked to examine the correlations among the compositions within the region. A correlation graph is hence plotted along the rows as shown in Fig. 3(c). The graph concludes visually that the inner structure among the neighbouring elements of these columns possesses highly spatial correlation. Such kind of observations are also found in many other practical scenarios. For example, when HyperSpectral Imaging (HSI) [22] technique is applied to collect approximately continuous spectral information of ground objects, the block-sparse structures with correlated coefficients are regularly appeared since the object generally covers a continuous portion of the scene. Source localization in EEG and MEG [23] falls into this category too, i.e., the coefficients in a certain source may be highly correlated. The other related examples involved several physiological signals, e.g. blood pressure, MECG, and temperature [24], which are all regarded to be typically spatially continuous signals with intra-block correlation. So far, there were several algorithms applied to solve the intra-block correlation problem. Typically, a method named Localized lOw rank Promoting (LOOP) was proposed in [25] for recovery of the block-sparse signals with intra-block correlated elements, in which the neighboring coefficients of the sparse signal were organized to form a number of 2 × 2 matrices. Clearly, these matrices were mostly not full-ranked. A log-determinant function was hence adopted to approximately suppress the matrices in low rank, and guaranteed its corresponding performance. The intra-block high-correlated entries act as the pieces of prior information, which guarantees the recovery of the block-sparse signals even without requiring the knowledge of the sizes and locations of the blocks.

In this work, a new approach of low-rank and block-sparse decomposition motivated by the LOOP method, named RPCA-LOOP, is proposed. Different from other RPCA methods, this new proposed model specially deals with the pixels of matrices in a correlative manner with their neighboring pixel values. To separate the principal components from the block-sparse outliers with the prior high correlation, an efficient optimization approach which can be mainly summarized as two steps, has been designed. In the first step, since the block-sparsity penalty term in the original optimization problem is replaced by a non-convex logarithmic function, many effective algorithms, applied to solve the traditional l₁-norm convex minimization problem, are not appropriate to be used here. Instead of directly minimizing the objective function, Majorization Minimization (MM) approach [26] was adopted to solve the non-convex minimization model. The MM approach suggests that iteratively minimizing a simple convex upper-bound can push the objective function downhill to the optimum. The related plausibility and operability also have been demonstrated in [27], [28]. Thus, the original model can be reformulated as a convex optimization surrogate. In the second step, to solve this convex minimization program, an effective iterative scheme based on the ADMM framework is employed to converge to the minimum of the original problem, and the optimization subproblems both have closed-form solutions. Moreover, a detailed convergence analysis of ADMM with respect to our new problem is presented to complete the study.

The main contributions of the paper are summarized as follows:

• A novel and effective model, called RPCA-LOOP, is proposed, which not only is an ordinary decomposition for both the low-rank component and the block-sparsity component, but also exploits spatial high-correlation structures existing in the columns of the block-sparsity component which can provide better performance in handling many practical applications.

• A MM approach is employed to approximate the original non-convex minimization model by a convex surrogate optimization problem, which can be solved with the unique closed-form solution by the polynomial-time iterative algorithms, i.e. ADMM. Particularly, the theoretical convergence analysis is also discussed by a clear proof in detail.

• A series of numerical experiments together with image denoising applications are taken to demonstrate the superiority of the proposed RPCA-LOOP.

The remainder of paper is organized as follows. In Section 2, we give some notations and preliminaries which are used in later section. Section 3 presents an efficient algorithm to solve RPCA problem with a block-sparse residual by adopting the ADMM. In Section 4, the convergence analysis of the new proposed algorithm is given in detail. Experimental results are presented in Section 5 to demonstrate the effectiveness of our method. Section 6 summarizes the conclusion and future work.