Thierry Bouwmans
email: tbouwman@univ-lr.fr

Andrews Sobral

Sajid Javed

Soon Ki

El-Hadi Zahzah

S Jung

S ||w ||l

L = Lw

No Bian

] Rosl

A = Dα

Decomposition into Low-rank plus Additive Matrices for Background/Foreground Separation: A Review for a Comparative Evaluation with a Large-Scale Dataset

Keywords: Background Subtraction, Foreground Detection, Robust Principal Component Analysis, Robust Non-negative Matrix Factorization, Robust Matrix Completion, Subspace Tracking, Low Rank Minimization Robust Problem Formulations Categories Sub-categories Authors -Dates

Background/foreground separation is the first step in video surveillance system to detect moving objects. Recent research on problem formulations based on decomposition into low-rank plus sparse matrices shows a suitable framework to separate moving objects from the background. The most representative problem formulation is the Robust Principal Component Analysis (RPCA) solved via Principal Component Pursuit (PCP) which decomposes a data matrix in a low-rank matrix and a sparse matrix. However, similar robust implicit or explicit decompositions can be made in the following problem formulations: Robust Non-negative Matrix Factorization (RNMF), Robust Matrix Completion (RMC), Robust Subspace Recovery (RSR), Robust Subspace Tracking (RST) and Robust Low-Rank Minimization (RLRM). The main goal of these similar problem formulations is to obtain explicitly or implicitly a decomposition into low-rank matrix plus additive matrices. These formulation problems differ from the implicit or explicit decomposition, the loss function, the optimization problem and the solvers. As the problem formulation can be NP-hard in its original formulation, and it can be convex or not following the constraints and the loss functions used, the key challenges concern the design of efficient relaxed models and solvers which have to be with iterations as few as possible, and as efficient as possible. In the application of background/foreground separation, constrainst inherent to the specificities of the background and the foreground as the temporal and spatial properties need

to be taken into account in the design of the problem formulation. Practically, the background sequence is then modeled by a low-rank subspace that can gradually change over time, while the moving foreground objects constitute the correlated sparse outliers. Although, many efforts have been made to develop methods for the decomposition into low-rank plus additive matrices that perform visually well in foreground detection with reducing their computational cost, no algorithm today seems to emerge and to be able to simultaneously address all the key challenges that accompany real-world videos. This is due, in part, to the absence of a rigorous quantitative evaluation with synthetic and realistic large-scale dataset with accurate ground truth providing a balanced coverage of the range of challenges present in the real world. In this context, this work aims to initiate a rigorous and comprehensive review of the similar problem formulations in robust subspace learning and tracking based on decomposition into low-rank plus additive matrices for testing and ranking existing algorithms for background/foreground separation. For this, we first provide a preliminary review of the recent developments in the different problem formulations which allows us to define a unified view that we called Decomposition into Low-rank plus Additive Matrices (DLAM). Then, we examine carefully each method in each robust subspace learning/tracking frameworks with their decomposition, their loss functions, their optimization problem and their solvers. Furthermore, we investigate if incremental algorithms and real-time implementations can be achieved for background/foreground separation. Finally, experimental results on a large-scale dataset called Background Models Challenge (BMC 2012) show the comparative performance of 32 different robust subspace learning/tracking methods.

Introduction

The detection of moving objects is the basic low-level operation in video analysis. This detection is usually done by using foreground detection. This basic operation consists of separating the moving objects called "foreground" from the static information called "background". Many foreground detection methods have been developed [30][26][27][28][31][269], and several implementations are available in the BGS Library [290]. Several foreground detection methods are based on subspace learning models such as Principal Component Analysis (PCA) [START_REF] Bouwmans | Subspace learning for background modeling: A survey[END_REF]. In 1999, Oliver et al. [233] were the first authors to model the background by PCA. Foreground detection is then achieved by thresholding the difference between the generated background image and the current image. PCA provides a robust model of the probability distribution function of the background, but not of the moving objects while they do not have a significant contribution to the model. Although there are several PCA improvements [310][311] that address the limitations of classical PCA with respect to outlier and noise, yielding to the field of robust PCA via outliers suppression, these methods do not possess the strong performance guarantees provided by the recent works on robust PCA via decomposition into low-rank plus sparse matrices [START_REF] Candes | Robust principal component analysis? International[END_REF][335] [START_REF] Chandrasekaran | Ranksparsity incoherence for matrix decomposition[END_REF] [START_REF] Abdel-Hakim | Distortion impact on low-dimensional manifold recovery of high-dimensional data[END_REF]. The idea of this recent RPCA approach is that the Fig. 1 RPCA via decomposition into low-rank plus sparse matrices in foreground/background separation: Original image (309), low-rank matrix L (background), sparse matrix S (foreground), foreground mask (Sequences from BMC 2012 dataset [315]).

data matrix A can be decomposed into two components such that A = L+S, where L is a low-rank matrix and S is a matrix that can be sparse. The decomposition into low-rank plus additive matrices are used in similar problem formulations such as Robust Non-negative Matrix Factorization (RNMF), Robust Matrix Completion (RMC), Robust Subspace Recovery (RSR), Robust Subspace Tracking (RST) and Robust Low-Rank Minimization (RLRM) [375]. RNMF assumes that the matrix L is a non-negative matrix. Sparsity constraints are applied on S in RPCA and not in LRM. Furthermore, changes can be tracked in the subspace, that is the field of subspace tracking. Applying RPCA via decomposition into low-rank plus sparse matrices in video-surveillance, the background sequence is modeled by the low-rank subspace that can gradually change over time, while the moving foreground objects constitute the correlated sparse outliers. For example, Fig. 1 shows original frames of sequences from the BMC dataset [315] and its decomposition into the low-rank matrix L and sparse matrix S. We can see that L corresponds to the background whereas S corresponds to the foreground. The fourth image shows the foreground mask obtained by thresholding the matrix S. So, the different advances in the different problem formulations of the decomposition into low-rank plus additive matrices are fundamental and can be applied to background modeling and foreground detection in video surveillance [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF] [START_REF] Bouwmans | Traditional and recent approaches in background modeling for foreground detection: An overview[END_REF].

The rest of this introduction is organized as follows. Firstly, we provide a preliminary overview of the different problem formulations for the robust subspace learning/tracking frameworks which used the decomposition into low-rank plus additive matrices. Then, we present a unified view of decomposition into low-rank plus additive matrices with a discussion about its adequation for the application of background/foreground separation. Thus, we review quickly similar decompositions such as sparse and mixed decompositions. Finally, we introduce our motivations to provide this review for a comparative evaluation in the application of background/foreground separation.

Problem Formulations based on Decomposition into Low-rank plus Additive Matrices: A Preliminary Overview

The aim of this section is to allow the reader to have a quick preliminary overview of the different robust problem formulations that are reviewed in details in the different sections of this paper. These different problem formulations based on an implicit or explicit decomposition into low-rank plus additive matrices are classified in the following categories: robust PCA, robust non-negative matrix factorization, robust subspace recovery, robust subspace tracking, robust matrix completion and robust low-rank minimization. 1.1.1 Robust Principal Component Analysis (RPCA) Recent research in robust PCA are based on the explicit decomposition into lowrank plus sparse matrices which differs from the decomposition, the loss functions, the optimization problem and the solvers used. These different approaches can be classified as follows:

1. RPCA via Principal Component Pursuit (RPCA-PCP): The first work on RPCA-PCP developed by Candes et al. [40][335] and by Chandrasekharan et al. [48] proposed the robust PCA problem as one of separating a low-rank matrix L (true data matrix) and a sparse matrix S (outliers' matrix) from their sum A (observed data matrix). Thus, a convex optimization allowed them to address the robust PCA problem. Under minimal assumptions, this approach called Principal Component Pursuit (PCP) perfectly recovers the low-rank and the sparse matrices. The background sequence is then modeled by a low-rank subspace that can gradually change over time, while the moving foreground objects constitute the correlated sparse outliers. Therefore, Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] showed visual results on foreground detection that demonstrated encouraging performance but PCP presents several limitations for foreground detection. The first limitation is that it required algorithms to be solved that are computational expensive. The second limitation is that PCP is a batch method that stacked a number of training frames in the input observation matrix.

In real-time application such as foreground detection, it would be more useful to estimate the low-rank matrix and the sparse matrix in an incremental way quickly when a new frame comes rather than in a batch way. The third limitation is that the spatial and temporal features are lost as each frame is considered as a column vector. The fourth limitation is that PCP imposed the low-rank component being exactly low-rank and the sparse component being exactly sparse but the observations such as in video surveillance are often corrupted by noise affecting every entry of the data matrix. The fifth limitation is that PCP assumed that all entries of the matrix to be recovered are exactly known via the observation and that the distribution of corruption should be sparse and random enough without noise. These assumptions are rarely verified in the case of real applications because of the following main reasons:

(1) only a fraction of entries of the matrix can be observed in some environments, [START_REF] Abdel-Hakim | Distortion impact on low-dimensional manifold recovery of high-dimensional data[END_REF] the observation can be corrupted by both impulsive and Gaussian noise, and (3) the outliers i.e moving objects are spatially localized. Many efforts have been recently concentrated to develop low-computational algorithms for solving PCP [188][189] [38][368][352][281][335][221][194][190][205] [START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex function[END_REF].

Other authors investigated incremental algorithms of PCP to update the lowrank and sparse matrix when a new data arrives [242][244][243][245]. Real-time implementations [8][9][306] have been developed too. Moreover, other efforts have addressed problems that appear specifically in real applications as background/foreground separation:

(a) Presence of noise: Noise in image is due to a poor quality image source such as images acquired by a web cam or images after compression. (b) Quantization of the pixels: The quantization can induce at most an error of 0.5 in the pixel value. (c) Spatial and temporal constraints of the foreground pixels: Low-rank and sparse decomposition is based on the condition that the outlier/noise can be considered as sparsity patterns and are uniformly located at the scene, which is not realistic in realworld applications as foreground moving objects are located in a connexed area. Furthermore, foreground moving objects present a continuous motion through the sequence. These two points need to introduce spatial and temporal constraints on the detection. (d) Local variations in the background: Variations in the background may be due to a camera jitter or dynamic backgrounds.

To address (a), Zhou et al. [393] proposed a stable PCP (SPCP) that guarantees accurate recovery in the presence of entry-wise noise. Becker et al. [20] proposed a inequality constrained version of PCP to take into account the quantization error of the pixel's value (b). To address (c), Tang and Nehorai [306] proposed a PCP method via a decomposition that enforces the low-rankness of one part and the block sparsity of the other part. Wohlberg et al. [333] used a decomposition corresponding to a more general underlying model consisting of a union of low-dimensional subspaces for local variation in the background (d). Furthermore, RPCA is generally applied in the pixel domain by using intensity or color features but other features can be used such as depth [141] and motion (optical flow [266]) features. Furthermore, RPCA can been extended to the measurement domain, rather than the pixel domain, for use in conjunction with compressive sensing [329][330][169][152][153][350][395][181][180][156][155][252].

Although experiments show that moving objects can be reliably extracted by using a small amount of measurements, we have limited the investigation and the comparative evaluation in this paper to the pixel domain to compare with the classical background subtraction methods.

2. RPCA via Outlier Pursuit (RPCA-OP): Xu et al. [342] proposed a robust PCA via Outlier Pursuit to obtain a robust decomposition when the outliers corrupted entire columns, that is every entry is corrupted in some columns. Moreover, Xu et al. [342] proposed a stable OP (SOP) that guarantee stable and accurate recovery in the presence of entry-wise noise.

3. RPCA via Sparsity Control (RPCA-SpaCtrl): Mateos and Giannakis [213][214] proposed a robust PCA where a tunable parameter controls the sparsity of the estimated matrix, and the number of outliers as a by-product.

4. RPCA via Sparse Corruptions (RPCA-SpaCorr): Even if the matrix A is exactly the sum of a sparse matrix S and a low-rank matrix L, it may be impossible to identify these components from the sum. For example, the sparse matrix S may be low-rank, or the low-rank matrix L may be sparse. To address this issue, Hsu et al. [130] imposed conditions on the sparse and low-rank components in order to guarantee their identifiability from A.

5. RPCA via Log-sum heuristic Recovery (RPCA-LHR): When the matrix has high intrinsic rank structure or the corrupted errors become dense, the convex approaches may not achieve good performances. Then, Deng et al. [START_REF] Deng | Low-rank structure learning via log-sum heuristic recovery[END_REF] used the log-sum heuristic recovery to learn the low-rank structure.

6. RPCA via Iteratively Reweighted Least Squares (IRLS): Guyon et al. [113] proposed to solve the RPCA problem by using an Iteratively Reweighted Least Squares (IRLS) alternating scheme for matrix low-rank decomposition. Furthermore, spatial constraint can be added in the minimization process to take into account the spatial connexity of pixels [112]. The advantage of IRLS over the classical solvers is its fast convergence and its low computational cost. Furthermore, Guyon et al. [111] improved this scheme by addressing in the minimization the spatial constraints and the temporal sparseness of moving objects.

7. RPCA via Stochastic Optimization (RPCA-SO): Goes et al. [START_REF] Goes | Robust stochastic principal component analysis[END_REF] proposed a robust PCA via a stochastic optimization. Feng et al. [83] developed an online Robust PCA (OR-PCA) that processes one sample per time instance and hence its memory cost is independent of the number of samples, significantly enhancing the computation and storage efficiency. The algorithm is equivalent to a reformulation of the batch RPCA [START_REF] Goes | Robust stochastic principal component analysis[END_REF]. Therefore, Javed et. al [144] modified OR-PCA via stochastic optimization method to perform it on background subtraction. An initialization scheme is adopted which converges the algorithm very fastly as compared to original OR-PCA. Therefore, OR-PCA was further improved to enhance the foreground segmentation using the continuous constraints such as Markov Random Filed (MRF) [147] and using dynamic feature selection [146].

8. RPCA with Dynamic Mode Decomposition (RPCA-DMD): Grosek et al. [START_REF] Grosek | Dynamic mode decomposition for real-time background/foreground separation in video[END_REF] introduced the use of dynamic mode decomposition (DMD) for robustly separating videoframes into background and foreground components in realtime. DMD [START_REF] Grosek | Robust real-time image processing through dynamic mode decomposition[END_REF] is a technique used for characterizing nonlinear dynamical systems in an equation-free manner by decomposing the state of the system into low-rank terms whose Fourier components in time are known. DMD terms with Fourier frequencies near the origin (zero-modes) are interpreted as background portions of the given video frames, and the terms with Fourier frequencies bounded away from the origin are their sparse counterparts. An approximate low-rank/sparse separation is achieved at the computational cost of one singular value decomposition and one linear equation solve, thus producing results orders of magnitude faster than RPCA. Kutz et al. [167][166] [START_REF] Grosek | Dynamic mode decomposition for robust pca with applications to foreground/background subtraction in video streams[END_REF] improved this approach for robustly separating background and foreground component into a hierarchy of multi-resolution time-scaled component.

9.

Bayesian RPCA (BRPCA): Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] proposed a Bayesian framework which infers an approximate representation for the noise statistics while simultaneously inferring the low-rank and sparse components. Furthermore, Markov dependency is introduced spatially and temporarily between consecutive rows or columns corresponding to image frames. This method has been improved in a variational Bayesian framework [17] and a factorized variational Bayesian framework [5]. In a similar manner, Zhao et al. [385] developed a generative RPCA model under the Bayesian framework by modeling data noise as a mixture of Gaussians (MoG).

10. Approximated RPCA: Zhou and Tao [389] proposed an approximated lowrank and sparse matrix decomposition. This method called Go Decomposition (GoDec) produces an approximated decomposition of the data matrix whose RPCA exact decomposition does not exist due to the additive noise, the predefined rank on the low-rank matrix and the predefined cardinality of the sparse matrix. GoDec is significantly accelerated by using bilateral random projection. Furthermore, Zhou and Tao [389] proposed a Semi-Soft GoDec which adopts soft thresholding to the entries of S, instead of GoDec which imposes hard thresholding to both the singular values of the low-rank part L and the entries of the sparse part S.

11. Sparse Additive Matrix Factorization: Nakajima et al. [222] [223] developed a framework called Sparse Additive Matrix Factorization (SAMF). The aim of SAMF is to handle various types of sparse noise such as row-wise and columnwise sparsity, in addition to element-wise sparsity (spiky noise) and low-rank sparsity (low-dimensional). Furthermore, their arbitrary additive combination is also allowed. In the original robust PCA [START_REF] Candes | Robust principal component analysis? International[END_REF], row-wise and column-wise sparsity can capture noise observed only in the case when some sensors are broken or their outputs are unreliable. SAMF due to its flexibility in sparsity design incorporate side information more efficiently. In background/foreground separation, Nakajima et al. [222] [223] induced the sparsity in SAMF using image segmentation.

12. Variational Bayesian Sparse Estimator: Chen et al. [START_REF] Chen | Variational Bayesian methods for multimedia problems[END_REF] proposed a variational Bayesian Sparse Estimator (VBSE) based algorithm for the estimation of the sparse component of an outlier corrupted low-rank matrix, when linearly transformed composite data are observed. It is a generalization of the original robust PCA [START_REF] Candes | Robust principal component analysis? International[END_REF]. VBSE can achieved background/foreground separation in blurred and noisy video sequences.

1.1.2 Robust Non-negative Matrix Factorization (RNMF)

Non-negative matrix factorization (NMF) approximates a non-negative matrix A by a product of two non-negative low-rank factor matrices W and H. Classical NMF methods minimize either the Euclidean distance or the Kullback-Leibler divergence between X and W T H to model the Gaussian noise or the Poisson noise.

Practically, these methods do not perform well when the noise distribution is heavy tailed as in the background/foreground separation. To address this problem, Guan et al. [99] proposed to minimize the Manhattan distance between X and W T H. This method called Manhattan NMF (MahNMF) robustly estimates the low-rank part and the sparse part of a non-negative matrix and performs effectively when data are contaminated by outliers. MahNMF shows similar qualitative performance as the RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF]. In another way, Kumar et al. [164] proposed a family of conical hull finding procedures called Xray for Near-separable NMF (NS-NMF) problems with Frobenius norm loss. However, in the presence of outliers or different noise characteristics, the use of Frobenius norm approximations is not optimal. Then, Kumar and Sindhwani [163] improved Xray to provide robust factorizations with respect to the 1 1 -loss function, and approximations with respect to the family of Bregman divergences. This algorithm is called RobustXray. Quantitative results [163] show that RobustXray outperforms the RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] in presence of noise. Another approach developed by Woo and Park [334] used a formulation called l∞-norm based robust asymmetric nonnegative matrix factorization (RANMF) for the grouped outliers and low nonnegative rank separation problems. The main advantage of RANMF is that the denseness of the low nonnegative rank factor matrices can be controlled. Furthermore, RANMF is not sensitive to the nonnegative rank constraint parameter due to the soft regularization method.

Robust Matrix Completion (RMC)

The matrix completion aims at recovering a low-rank matrix from partial observations of its entries. Robust matrix completion RMC, also called RPCA plus matrix completion problem can also be used for background/foreground separation. RPCA via principal component pursuit [START_REF] Candes | Robust principal component analysis? International[END_REF] can be considered as RMC using l 1 -norm loss function. Following this idea, Yang et al. [360] proposed a nonconvex relaxation approach to the matrix completion problems when the entries are contaminated by non-Gaussian noise or outliers. A nonconvex loss function based on the lσ-norm instead of the l 1 -norm is used with a rank constrained as well as a nuclear norm regularized model. This method can be solved via two algorithms based on iterative soft thresholding (IST) and iterative hard thresholding (IHT). A nonconvex loss function used in robust statistics is used with a rank constrained as well as a nuclear norm regularized model. This method called RMC-lσ-IHT is also faster than RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF]. In another way, Shang et al. [275][273] proposed a scalable, provable structured low-rank matrix factorization method to recover low-rank plus sparse matrices from missing and grossly corrupted data. A scalable robust bilinear structured factorization (RBF) method recovered low-rank plus sparse matrices from incomplete, corrupted data or a small set of linear measurements. In a similar way, Shang et al. [274] proposed a scalable convex model (RMC with convex formulation) and a non-convex model solved with matrix factorization (RMC-MF) in which the desired low-rank matrix L is factorized into two much smaller matrices. In another way, Mansour and Vetro [212] proposed a factorized robust matrix completion (FRMC) algorithm with global motion compensation. The algorithm decomposes a sequence of video frames into the sum of a low-rank background component and a sparse motion component. The algorithm alternates between the solution of each component following a Pareto curve tra-jectory for each subproblem. For videos with moving background, Mansour and Vetro [212] used the motion vectors extracted from the coded video bitstream to compensate for the change in the camera perspective. This approach is faster than state-of-the-art solvers and results in highly accurate motion segmentation. In a similar way, Yang et al. [351] proposed a Motion-Assisted Matrix Completion (MAMC) which used a dense motion field for each frame. This field is then mapped into a weighting matrix to assign the reliability of pixels that belong to the background. This method is robust to slowly moving objects and camouflage. Yang et al. [351] extended MAMC to a robust MAMC model (RMAMC) which is robust to noise. Sobral et al. [291] provided a comparison of several matrix completion algorithms on the SBMI dataset [210].

1. 1.4 Robust Subpace Recovery (RSR) This category contains the robust decompositions other than RPCA and RNMF decompositions. First, Wang et al. [325] studied the problem of discovering a subspace in the presence of outliers and corruptions. In this context, additional knowledge is added to relax this problem as a convex programming problem. Thus, Wang et al. [325] provided a Robust Subspace Discovery (RSD) method solved via an efficient and effective algorithm based on the Augmented Lagrangian Multiplier.

Since high dimensional data is supposed to be distributed in a union of low dimensional subspaces, Bian and Krim [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] proposed a bi-sparse model as a framework to take into account that the underlying structure may be affected by sparse errors and/or outlier. So, Bian and Krim [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] provided an algorithm called Robust Subspace Recovery via bi-sparsity pursuit (RoSuRe) to recover the union of subspaces in presence of sparse corruptions. Experimental results [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] show robustness in the case of camera jitter. Conventional robust subspace recovery models address the decomposition problem by iterating between nuclear norm and sparsity minimization. However, this scheme is computationally prohibitive to achieve real time requirements. To solve this problem, Shu et al. [287] proposed a Robust Orthogonal Subspace Learning (ROSL) method to achieve efficient low-rank recovery. A rank measure on the low-rank matrix is introduced that imposes the group sparsity of its coefficients under orthonormal subspace. Furthermore, an efficient sparse coding algorithm minimizes this rank measure and recovers the low-rank matrix at quadratic complexity of the matrix size. Finally, Shu et al. [287] developed a random sampling algorithm to further speed up ROSL such that its accelerated version (ROSL+) has linear complexity with respect to the matrix size. Experiments [287] demonstrate that both ROSL and ROSL+ provide more efficiency against RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] with the same detection accuracy. In a different manner, She et al. [277] proposed a robust orthogonal complement principal component analysis (ROC-PCA). The aim is to deal with orthogonal outliers that are not necessarily apparent in the original observation space but could affect the principal subspace estimation. For this, She et al.

[277] introduced a projected mean-shift decomposition and developed a fast alternating optimization algorithm on the basis of Stiefel manifold optimization and iterative nonlinear thresholdings. Subspace tracking aims to address the problem when new observations come in asynchronously in the case of online streaming application. The algorithm cannot store in general all the input data in memory. The incoming observations must be immediate processed and then discarded. Furthermore, since the subspace can be identified from incomplete vectors, it can be subsampled in order to improve on computational efficiency, and it still retain subspace estimation accuracy. The involved subspaces can have low-rank and/or sparse structures as in the previous decomposition problem formulations. In this idea, He et al. [119][120] proposed an incremental gradient descent on the Grassmannian, the manifold of all d-dimensional subspaces for fixed d. This algorithm called Grassmannian Robust Adaptive Subspace Tracking Algorithm (GRASTA) uses a robust l 1 -norm cost function in order to estimate and track non-stationary subspaces when the streaming data vectors, that are image frames in foreground detection, are corrupted with outliers, that are foreground objects. This algorithm allows to separate background and foreground online. GRASTA shows high-quality visual separation of foreground from background. Following the idea of GRASTA, He et al. [122][123] proposed (transformed-GRASTA) which iteratively performs incremental gradient descent constrained to the Grassmannian manifold of subspaces in order to simultaneously estimate a decomposition of a collection of images into a low-rank subspace, a sparse part of occlusions and foreground objects, and a transformation such as rotation or translation of the image. t-GRASTA is four times faster than state-of-the-art algorithms, has half of the memory requirement, and can achieve alignment in the case of camera jitter. Although the l 1 -norm in GRASTA leads to favorably conditioned optimization problems it is well known that penalizing with non-convex l 0 -surrogates allows reconstruction even in the case when l 1 -based methods fail. Therefore, Hage and Kleinsteuber [115][268] proposed an improved GRASTA using l 0 -surrogates solving it by a Conjugate Gradient method. This method called pROST [268] outperforms GRASTA in the case of multi-modal backgrounds. Another approach developed by Xu et al. [343] consists of a Grassmannian Online Subspace Updates with Structured-sparsity (GOSUS), which exploits a meaningful structured sparsity term to significantly improve the accuracy of online subspace updates. Their solution is based on Alternating Direction Method of Multipliers (ADMM), where most key steps in the update procedure are reduced to simple matrix operations yielding to real-time performance. Finally, Ahn [4] proposed a fast adapted subspace tracking algorithm which shares the procedure of separating frames into background and foreground with GRASTA, but it uses a recursive least square algorithm for subspace tracking, which makes it fast adapted to dynamic backgrounds. Low-rank minimization (approximation or representation) is a minimization problem, in which the cost function measures the fit between a given data matrix A and an approximating matrix L, subject to a constraint that the approximating matrix L has reduced rank. In the application of background/foreground separation, Zhou et al. [391] proposed a framework called Detecting Contiguous Outliers in the Low-Rank Representation (DECOLOR) and formulated outlier detection in the robust low-rank representation, in which the outlier support and the low-rank matrix are estimated. This approach integrates object detection and background learning into a single process of optimization solved by an alternating algorithm. In a different manner, Xiong et al. [341] proposed a direct robust matrix factorization (DRMF) assuming that a small portion of the matrix A has been corrupted by some arbitrary outliers. The aim is to get a reliable estimation of the true low-rank structure of this matrix and to identify the outliers. To achieve this, the outliers are excluded from the model estimation. Furthermore, Xiong et al. [341] proposed an extension of DRMF to deal with the presence of outliers in entire columns. This method is called DRMF-Row (DRMF-R). In another way, Wang et al. [321] proposed a probabilistic method for robust matrix factorization (PRMF) based on the l 1 -norm loss and l 2 -regularizer, which bear duality with the Laplace error and Gaussian prior, respectively. But, PRMF treats each pixel independently with no clustering effect but the moving objects in the foreground usually form groups with high within-group spatial or temporal proximity. Furthermore, the loss function is defined based on the l 1 -norm, and its results are not robust enough when the number of outliers is large. To address these limitations, Wang et al. [366] proposed a full Bayesian formulation called Bayesian Robust Matrix Factorization (BRMF). BRMF used a Laplace mixture with the generalized inverse Gaussian distribution as the noise model to further enhance model robustness. Furthermore, BRMF contained a Markov extension (MBRMF) which assumes that the outliers exhibit spatial or temporal proximity. In a different manner, Zheng et al.

[386] added a convex nuclear-norm regularization term to improve convergence of LRM, without introducing too much heterogenous information. This method is called Practical Low-Rank Matrix Factorization (PLRMF). The previous low-rank factorization used loss functions such as the l 2 -norm and l 1 -norm losses. l 2 -norm is optimal for Gaussian noise, while l 1 -norm is for Laplacian distributed noise. Since videos are often corrupted by an unknown noise distribution, which is unlikely to be purely Gaussian or Laplacian, Meng et al. [216] used a low-rank matrix factorization problem with a Mixture of Gaussians (LRMF-MoG) noise model. Since the MoG model is a suitable approximator for any continuous distribution, it is able to model a wider range of noise.

Table 1 and Table 2 show an overview of the different problem formulations based on decomposition into low-rank plus additive matrices. The first column indicates the name of the different problem formulations and the second column shows the different categories of each problem formulations. The third column indicates the different methods of each category with their corresponding acronym. The fourth column gives the name of the authors and the date of the related publication. The previous surveys in the field are indicated in bold and the reader can refer to them for more references on the corresponding category or sub-category. Furthermore, we present in different tables some quick comparisons on the different key characteristics of these different problem formulations based on the decomposition into low-rank plus additive matrices. Thus, Table 3, Table 4, Table 5 and Table 6 show an overview of the different decompositions into low-rank plus additive matrices in terms of minimization, constraints and convexity to allow us to define a unified view that we called Decomposition into Low-rank plus Additive Matrices (DLAM).

A = L + S min L,S ||L|| * + λ||S|| l 1 A -L -S = 0 Yes
Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] Modified PCP (Fixed Rank)

A = L + S min L,S
||L|| * + λ||S|| l 1 rank(L) = known r Yes Leow et al. [174] Modified PCP (Nuclear Norm Free)

A = u1 T + S min u ||A -u1 T || l 1 rank(u1 T ) = 1 Yes Yuan et al. [367] Modified PCP (Capped Norm) A = L + S min L,S rank(L) + λ||S|| l 0 ||A -L -S|| 2 F ≤ σ 2 No Sun et al. [299]
Modified PCP (Inductive)

A = P A + S min P,S ||P || * + λ||S|| l 1 A -P A -S = 0 Yes
Bao et al. [18] Modified PCP (Partial Subspace Knowledge) Zhan and Vaswani [369] p,q-PCP (Schatten-p norm, lq norm) Wang et al.[320] Modified p,q-PCP (Schatten-p norm, Lq seminorm) Zhou and Tao [390] Modified SPCP (Nuclear Norm Free) Yuan et al. [367] Modified SPCP (Nuclear Norm Free for blur in video) Yuan et al. [367] Modified SPCP (Undercomplete Dictionary)

A = L + S min L,S ||L|| * + λ||S|| l 1 L + P Γ ⊥ S = P Γ ⊥ A Yes
A = L + S min L,S ||L|| p Sp + λ||S|| lq A = L + S No
A = L + S min L,S ||L|| p Sp + λ||S|| q Lq A = L + S No Shao et al.[276] Modified PCP (2D-PCA) A = L + S min U,V 1 T T i=1 ||Ai -UΣiV T || 2 F U T U = Ir×r , V T V = Ic×c No Sun et al.[301] Modified PCP (Rank-N Soft Constraint) A = L + S min L,S min(m,n) i=N +1 |σi(L)| + λ||S|| l 1 A = L + S Yes
A = UV + S + E min U,V,S λ||S||l 1 + ||A -UV -S|| 2 F rank(U) = rank(S) ≤ r Yes
A = U1 T + S + E min S∈Rm×n ,u∈Rm ||S||l 1 + µ 2 ||A -u1 T -S|| 2 F rank(U1 T ) = 1 Yes
A = U1 T + S + E min S∈Rm×n ,u∈Rm ||S||l 1 + µ 2 ||A -H(u1 T + S)|| 2 F rank(U1 T ) = 1 Yes
A = UT + S + E min U,T ,S λ||S||l 1 + λ1 2 (||U|| 2 F + ||T || 2 F ) + 1 2 ||A -UT -S)|| 2 F rank(UT ) ≤ r Yes Sprechman et al. [295]
Variational SPCP (Huber penalty) [START_REF] Aravkin | A variational approach to stable principal component pursuit[END_REF] Modified SPCP (Three Term Low-rank Optimization) [START_REF] Mansour | Video background subtraction using semi-supervised robust matrix completion[END_REF] [212] Motion-Assisted Matrix Completion Yang et al. [351] Table 5 Decompositions in low-rank plus additive matrices: An homogeneous overview (Part 3)

A = L + S + E min Φ(L, S) ρ(L + S -Y ) < ǫ Yes Aravkin et al.
A = L + S + E min A,L,S λ||L|| * + λ1||fπ(S)|| + λ2||E|| 2 F A = L + S + E Yes Oreifej et al. [236] QPCP A = L + S min L,S ||L|| * + λ||S||l 1 ||A -L -S||∞ < 0.5 Yes Becker et al. [20] BPCP A = L + S min L,S ||L|| * + κ(1 -λ)||L||l 2,1 + κλ||S||l 2,1 A -L -S = 0 Yes Tang and Nehorai [306] LPCP A = AU + S min U,S α||U||l 1 + β||U||l 2,1 + β||S||l 1 A -AU + S = 0 Yes Wohlberg et al. [333] RPCA-OP OP A = L + S min L,S ||L|| * + λ||S||l 1,2 A -L -S = 0 Yes Xu et al. [342] SOP A = L + S + E min L,S ||L|| * + λ||S||l 1,2 ||A -L -S||F < δ Yes Xu et al. [342] RPCA-SpaCtrl Sparsity Control A = M + U T P + S + E min U,S ||X + 1N M T -P U T -S|| 2 F + λ||S||l 2(r) UU T = Iq Yes Mateos et Giannakis [213][214] RPCA-SpaCorr Sparse Corruptions (case 1) A = L + S min L,S ||L|| * + λ||S||l 1 ||A -L -S||l 1 ≤ ǫ1 Yes Hsu et al. [130] ||A -L -S|| * ≤ ǫ * ||L||∞ ≤ b Sparse Corruptions (case 2) A = L + S min L,S ||L|| * + λ||S||l 1 + 1 2µ ||A -L -S|| 2 F ||A -L -S||l 1 ≤ ǫ1 Yes Hsu et al. [130] ||A -L -S|| * ≤ ǫ * ||A -S||∞ ≤ b RPCA-LHR LHR A = L + S min X∈ D 1 2 (||Diag(Y )||L + ||Diag(Z)||L) + λ||S||L X = {Y, Z, L, S} No Deng et al. [66] D = (Y, Z, L, S) : Y L L T Z ≥ 0, (L, S) ∈ C RPCA-IRLS IRLS A = UV + S min U ∈R n×p ,V ∈R p×m µ||UV || * + ||(A -UV ) • W1||l α,β A -UV -S = 0 Yes Guyon et al.[113] RPCA-SO OR-RCPA A = LR + S + E min L∈R n×p ,R∈Rn×r 1 2 ||A -LR T -S|| 2 F + λ1 2 (||L|| 2 F + ||R|| 2 F ) + λ2||S||l 1 A -LR -S = 0 Yes Feng et al. [83]
A = UV T + RS + E ρ(A, U, V, S, γ, α, β) Distribution constraints - Chen et al. [223] Robust Non-negative Matrix Factorization Manhattan NMF A = W T H + S min W ≥0,H≥0 f (W, H) = ||A -W T H||M r ≪ min(m, n) No (MahNMF) Guan et al.[99] NS-NMF (Xray-l2) A = W H = AB H + S min A B ≥0,H≥0 ||A -AB H|| 2 F W = AB ≥ 0, H ≥ 0 Yes (Xray-l2) Kumar et al. [164] NS-NMF (RobustXray) A = W H + E = AB H + S + E min A B ≥0,H≥0 ||A -AB H||l 1 W = AB ≥ 0, H ≥ 0 Yes Kumar and Sindhwani [163] RANMF A = L + S = W ΛH + S min L,S,Φ Φ(S) + alpha 2 ||A -L -S|| 2 F + R(L) ≤ τ, 0 < L < BL Yes Woo and Park [334] β(Ψ (S, Φ) + γT V (Q(Φ)) Robust Matrix Completion RMC-lσ -IHT A = L + S min S∈R m×n lσ (L) + λ||L|| * rank(L) ≤ r No Yang et al.[360] RMC-RBF A = L + S = UV T + S min U,V,S ||PΩ (S)||1 + λ||V || * PΩ (A) = PΩ (UV T + S), U T U = I Yes Shang et al. [273] RMC (convex formulation) A = L + S min L,S ||L|| * + λ||PΩ (S)||l 1 PΩ (L + S) = PΩ (A), E Ω C = 0 Yes Shang et al. [274] RMC-MF (non-convex formulation) A = L + S = UV T + S min U,V,S ||V || * + λ||PΩ (S)||l 1 A = UV T + S, U T U = I No Shang et al. [274] Factorized Robust Matrix Completion A = L + S = LLL T R + S ||L|| * = min L L ∈R m,r ,L R ∈R n,r 1
2 (||LL|| 2 F + ||LR|| 2 F ) LLL T R = L Yes (FRMC)
A = L + S min L,S ||L|| * + λ||S||l 1 W3 • A = W3 • (L + S) Yes (MAMC)
||S||l 1 AΩ t • τ = UΩ t w + S , U ∈ G(d, n) No He et al. [122][123] GASG21 A = UW + S min w ||Uw -A||l 2,1 -m j=1 ||Uj wj -Aj ||l 2 A = Uw + S , U ∈ G(d, n) No He and Zhang [124] pROST A = UW + S + E min Rank(L)≤k ||UW -A||L p A -UW -S = 0 No Hage and Kleinsteuber [115][268] GOSUS A = UW + S + E min U T U =I d ,W,S l i=1 µi||D i S||l 2 + λ 2 ||UW + S -A|| 2 l2 A -UW -S = 0 No Xu et al. [343] FARST A = UW + S + E min W ||UW -A||l 1 A -UW -S = 0 No
Ahn [4] Robust Low-Rank Minimization Contiguous Outlier Detection Zheng et al. [386] Low Rank Matrix Factorization with MoG noise

A = L + S + E min L,S α||L|| * + β||F ||l 1 + γ||Cvec(F )|| l1 rank(L) ≤ K No (DECOLOR) Zhou et al. [391] + 1 2 ||P F (A -L)|| 2 F Direct Robust Matrix Factorization A = L + S min L,S ||A -S -L||F rank(L) ≤ r, ||S||0 ≤ p No (DRMF) Xiong et al. [341] Original formulation PCP [40] Direct Robust Matrix Factorization-Row A = L + S min L,S ||A -S -L||F rank(L) ≤ r, ||S||2,0 ≤ p No (DRMF-R) Xiong et al. [341] Original formulation OP [342] Probabilistic Robust Matrix Factorization A = UV ′ + S log (p(U, V |A, λ, λU , λV )) Distribution constraints No (PRMF) Wang et al. [321] Bayesian Robust Matrix Factorization A = UV ′ + S log (p(U, V |A, λ, λU , λV )) Bayesian distribution constraints No (BRMF) Wang et al. [366] Practical Low-Rank Matrix Factorization A = UV + S min U,V ||W5 ⊙ (A -UV )||l 1 + λ||V || * U T U = Ir Yes (PLMR)
A = UV T + S max U,V,Π,Σ i,j∈Ω K k=1 πkN(xij |(u i ) T v j , σ
= UV T + S min L,U,V ||W5 ⊙ (A -L)||l 1 + λ 2 (||U|| 2 F + λ||V || 2 F ) L = UV T Yes (UNN-BF) Cabral et al. [37] Robust Rank Factorization A = BX + S + E min B min S,X ||A -BX -S|| 2 l2 + λ 2 ||S||l 1 λ > 0, A = BX + S + E Yes (RRF) Sheng et al. [282]
Table 6 Decompositions in low-rank plus additive matrices: An homogeneous overview (Part 4)

1.2 A Unified View of Decomposition into Low-rank plus Additive Matrices

Notations

To provide to the readers an easy comparison, we homogenized all the different notations found on all the papers as follows:

1. Matrices: For the common matrices, A stands for the observation matrix, L is the low-rank matrix, S is the unconstrained (residual) matrix or sparse matrix, and E is the noise matrix. I is the identity matrix. For the specific matrices, the notations are given in the section of the corresponding method.

2. Indices : m and n are the number of columns and rows of the observed data matrix A. In the case of background/foreground separation, m corresponds to the number of pixels in a frame, and n corresponds to the number of frames in the sequence. n is taken usually to 200 due to computational and memory limitations. i and j stand for the current indices of the matrix. r is the estimated or fixed rank of the matrix L. p stands for the p th largest value in truncated matrix.

Norms:

The different norms used in the paper for vectors and matrices can be classified as follows:

-Vector l α -norm with 0 ≤ α ≤ 2 : ||V || l 0 is the l 0 -norm of the vector V , and it corresponds to the number of non-zero entries. ||V || l 1 = i v i is the l 1 -norm of the vector V , and it corresponds to the sum of the entries.

||V || l 2 = i (v i ) 2
is the l 2 -norm of the vector V , and it corresponds to the Euclidean distance [375].

-Matrix lα-norm with 0 ≤ α ≤ 2 : ||M || l0 is the l 0 -norm of the matrix M , and it corresponds to the number of non-zero entries [375]. 375], and its corresponds to the Manhattan distance [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF]. ||M || l2 = i,j M 2 i,j is the l 2 -norm of the matrix M also known as the Frobenius norm.

||M || l1 = i,j |M ij | is the l 1 -norm of the matrix M [
-Matrix l∞-norm: 374] is the l∞-norm of the matrix M . It can be used to capture the quantization error of the observed value of the pixel as in Becker et al. [20]. It is equivalent to the max-norm [280]. -Matrix Lα-seminorm with 0 < α ≤ 2:

||M || l∞ = m ij ax |M ij | [
-Matrix l α,β -norm with 0 ≤ α, β ≤ 2: ||M || l α,β is the l α,β -
||M || Lα = ( i,j |M ij | α ) 1/α is the Lα-seminorm of the matrix M [276]. The L 1 -seminorm is equivalent to the l 1 -norm .
-Matrix Lα-quasi-norm with 0 < α < 1: Lα-quasi-norm is defined by

Lα(M ) = m i=1 (M 2 i + µ) 1 α [115][268].
-Matrix Frobenius norm:

||M || F = i,j M 2
i,j is the Frobenius norm [375]. The Frobenius norm is sometimes also called the Euclidean norm which may cause confusion with the vector l 2 -norm which is also sometimes known as the Euclidean norm too. (σ k (M)) α ) 1/α where σ k (M), denotes the k th singular values of M , can also be used as a surrogate of the nuclear norm as in [320][276]. The Schatten-1-norm is equivalent to the nuclear norm [280].

-Matrix Log-sum norm: The Log-sum norm ||M || L is defined as ij log (|M i j|+ δ) with δ > 0 is a small regularization constant [START_REF] Deng | Low-rank structure learning via log-sum heuristic recovery[END_REF]. Frobenius norm [START_REF] Candes | Robust principal component analysis? International[END_REF] Inequality Most of the time, proxy loss functions are used as surrogate of the original loss function (rank(.)) loss function for the low-rank constraint and l 0 -loss function for the sparsity constraint) to obtain a solvable problem. From the homogenized overview, we can see that all the decompositions in the different problem formulations can be grouped in a unified view that we called Decomposition into Low-rank plus Additive Matrices (DLAM).

(||A -L -S|| 2 F ≤ σ 2 ) Equality (||A -L -S||F = 0) Frobenius norm [299] Modified SPCP Inequality (ρ(A -L -S) ≤ ǫ): Huber penalty [10] QPCP Inequality (|||A -L -S||l ∞ ≤ 0.5) l∞-norm [20] SPCP SPCP Inequality (||A -L -S||F < δ) Equality (||A -L -S -E||F = 0) Frobenius norm [393] + Frobenius norm [236] on E

Decomposition Problem

All the decompositions can be written in a general formulation as follows:

A = K k=1 M k = M 1 + M 2 + M 3 = L + S + E (1)
with K ∈ {1, 2, 3}. The matrix A can be corrupted by element-wise outliers/noise, column or row wise outliers/noise or missing data [208] as shown in Figure 2. The characteristics of the matrices M k are as follows:

-The first matrix M 1 = L is a low-rank matrix. In some decompositions, L is decomposed as follows:

(1) a product of two matrices U V obtained by bilateral projection [390] or by matrix factorization [383] in the RPCA framework, (2) a product of two matrices U V T [275][273][274] obtained by robust matrix factorization in the RMC framework, and (3) a product of two matrices W H [99][164][163] or three matrices W ΛH [334] with the constraints to be positive in the RNMF framework.

-The second matrix M 2 is an unconstrained (residual) matrix in implicit decomposition (LRM, RMC). M 2 is a sparse matrix S in explicit decomposition (RPCA, RNMF, RSR, RST). This second matrix can be decomposed as follows:

(1) a sum of two matrices S 1 and S 2 which correspond to the foreground and the dynamic backgrounds [START_REF] Cao | Total variation regularized RPCA for irregularly moving object detection under dynamic background[END_REF], and (2) a product of two matrices RT where R is an unconstrained matrix and T is a sparse matrix [57][58][56] in the framework of RPCA.

-The third matrix M 3 is generally the noise matrix E. The noise can be modeled by a Gaussian, a mixture of Gaussians (MoG) or a Laplacian distribution. M 3 can capture turbulences in the case of background/foreground separation [236].

Practically, the decomposition is A ≈ L when K = 1. This decomposition is called implicit as the second matrix S is not explicitly used in the problem formulation but it can be obtained by the difference A -L. It corresponds to basic formulations (LRM, MC, PCA and NMF). The decomposition is A = L + S when K = 2. It is called explicit because S is explicitly determined and used in the problem formulation. It corresponds to robust formulations (RLRM, RMC, RPCA, RNMF,RSR and RST). In the case of K = 3, the decomposition is A = L + S + E. This decomposition is called "stable" explicit decomposition as it separates the outliers in S and the noise in E. It corresponds to stabble robust formulations. Table 8 shows an overview of the matrix M k in the different problem formulations. Thus, we can see that the decomposition is implicit in LRM and MC because the decomposition is made implicitly, i.e A ≈ L. S is the residual matrix in RLRM Fig. 2 Illustration of different types of corruption fo the matrix A: a) A without noise b) A with element-wise outliers/noise and c) A with both element-wise outliers/noise and missing data (Illustration from the slides of Ma et al. [208] ).

and RMC whereas S is a sparse matrix in RPCA, RNMF, RSR and RST with constraints on E for the stable versions.

Problem Formulations Matrix M1 = L Matrix M2 = S Matrix M3 = E Implicit Decomposition (K = 1) - LRM L Low-rank Not used - MC L Low-rank Not used - Explicit Decomposition (K = 2) - RLRM L Low-rank S Residual matrix - RMC L Low-rank S Residual matrix - RPCA L Low-rank S Sparse (outlier+noise) - RNMF L = W H T Low-rank (positive) S Sparse (outlier+noise) - RSR L Low-rank S Sparse - RST L Low-rank S Sparse - Stable Explicit Decomposition (K = 3) - Stable RPCA L Low-rank S Sparse (outlier) E Noise Stable RNMF L = W H T Low-rank (positive) S Sparse (outlier) E Noise Stable RSR L Low-rank S Sparse E Noise Stable RST L Low-rank S Sparse E Noise
Table 8 Decomposition into Low-rank plus Additive Matrices (DLAM): The different matrices in the different problem formulations.

Minimization problem

The corresponding minimization problem of Equation 1 can be formulated in a general way as follows:

min Mi K k=1 λ i f i (M i ) subj C i (2)
where the λ i are the regularization parameters. f i () are the loss functions with

f 1 () = f low (), f 2 () = fsparse() and f 3 () = f noise . C i is a constraint on L, S and E with A which varies following the value of K, that is C 1 is ||A -L|| 2 = 0, C 2 is ||A -L -S|| 2 = 0, C 3 is ||A -L -S -E|| 2 = 0. C i can
be expressed in an inequality form too.

Case K = 1: It is the degenerated case for the implicit decompositions (LRM, MC) which are not robust because there are no constraints on the matrix S = A-L. In this case, Equation 2 with K = 1 can be written as follows:

min L λ 1 f low (L) subj C 1 (3)
Case K = 2: Equation 2 with K = 2 for the explicit decompositions becomes as follows: min

L,S λ 1 f low (L) + λ 2 fsparse(S) subj C 2 (4)
where 7. This minimization problem can be N P -hard, and convex or not following the constraints and the loss functions used. Practically, when the problem is N P -hard and/or not convex, the constraints are relaxed and the loss functions are changed to obtain a tractable and convex problem. For example, the original formulation in RPCA [START_REF] Candes | Robust principal component analysis? International[END_REF] used the rank(.) and the l 0 -norm as loss functions for L and S, respectively as shown in Equation 21. As this problem is N P -hard, this formulation is relaxed with the nuclear norm and the l 1 -norm as shown in Equation 22. To minimize confusion, the models that minimize rank functions and nuclear norms are named the original model and the relaxed model, respectively. Thus, the corresponding minimization problem of Equation 4 can be formulated with norms to be convex and solvable as follows: min

L,S λ 1 ||L|| p1 norm1 + λ 2 ||S|| p2 norm2 subj C 2 (5)
where λ 1 and λ 2 are the regularization parameters. [START_REF] Abdel-Hakim | FRPCA: fast robust principal component analysis[END_REF] an equality such as ||A -L -S|| p0 norm0 = 0 or rank(L) = r, or (2) an inequality such as ||A -L -S|| p0 norm0 ≤ q or rank(L) ≤ r. ||.||norm 0 is a norm taken in the set of norms previously defined. Moreover, the minimization problem formulated in Equation 5 can be written in its Lagrangian form as follows:

min L,S λ 0 2 ||A -L -S|| p0 norm0 + λ 1 ||L|| p1 norm1 + λ 2 ||S|| p2 norm2 subj C L2 (6)
where λ 0 is regularization parameter. C L2 is the constraint similar to the constraint C 2 .

Case K = 3: Equation 2 with K = 3 for the stable explicit decomposition is written as follows: min

L,S λ 1 f low (L) + λ 2 fsparse(S) + λ 3 f noise (E) subj C 3 (7)
where the λ i are the regularization parameters. C 3 is a constraint on A, L, S and E. f noise (E) is a function that contrains E. Thus, the corresponding minimization problem of Equation 7 can be formulated with norms to be convex and solvable as follows:

min

L,S λ 1 ||L|| p1 norm1 + λ 2 ||S|| p2 norm2 + λ 3 ||E|| p3 norm3 subj C 3 ( 8
)
where p 3 is taken in the set {1, 2}. ||.||norm 3 could be any previous norms. Then, the minimization problem formulated in Equation 8 can be written in its Lagrangian form as follows:

min

L,S λ 0 2 ||A -L -S -E|| p0 norm0 + λ 1 ||L|| p1 norm1 + λ 2 ||S|| p2 norm2 λ 3 ||E|| p3 norm3 subj C L3 (9
) where λ 0 is regularization parameter. C L3 is the constraint similar to the constraint C 3 . Finally, the minimization problem seeks to the following optimization problem:

F (X) = f (x) + g(x)
where we have:

f (x) = λ 0 2 ||A -L -S -E|| p0 norm0 (10) g(x) = λ 1 ||L|| p1 norm1 + λ 2 ||S|| p2 norm2 + λ 3 ||E|| p3 norm3 ( 11
)
where we have:

-f (x) is smooth, convex, and has Lipschitz continuous gradients [349].

-g(x) can be a nonsmooth, nonconvex function [349].

In general, solving a nonsmooth, nonconvex objective function is difficult with weak convergence guarantees [349].

Algorithms for solving the optimization problem

Algorithms which are called solvers are then used to solve the minimization problem in its original form or in its Lagrangian form. Furthermore, instead of directly solving the original convex optimizations, some authors use their strongly convex approximations in order to design efficient algorithms. Zhang et al. [373] proved that these strongly convex programmings guarantee the exact low-rank matrix recovery as well. Moreover, solvers have different characteristics in terms of complexities: complexity per iteration, complexity to reach an accuracy of ǫ precision (ǫ-optimal solution), and convergence rate complexity following the number of iterations. The key challenges related to the solvers are the following ones [207]:

1. Choice of the solver to make the iterations as few as possible.

2. Choice of the SVD algorithm to make the iterations as efficient as possible.

The solvers can be broadly classified into two categories as developed by Chen [START_REF] Chen | Multidimensional signal processing for sparse and low-rank problems[END_REF]:

-Regularization based approaches: The decomposition is formulated as regularized fitting, where the regularizers are convex surrogates for rank and sparsity. Examples of algorithms in this category include the following solvers: Singular Value Thresholding (SVT) [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], the Accelerated Proximal Gradient (APG) [189], and the Augmented Lagrange Multiplier (ALM) [188]. All the

Solvers for PCP Complexity

Basic solvers Singular Value Threshold (SVT1 ) Oiter (mnmin(mn)), Opre=unknown, Oconv=unknown Cai et al. [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] Iterative Thresholding (IT)

Oiter (mnmin(mn)), Opre( L/ǫ), Oconv = 1/T Oiter =Unknown, Opre=Unknown, Oconv =Unknown Chai et al. [START_REF] Chai | An improved ADM algorithm for RPCA optimization problem[END_REF] Table 9 Solvers for RPCA-PCP: An overview of their complexity per iteration at running time O iter , their complexity Opre to reach an accuracy of ǫ precision and their convergence rate Oconv for T iterations. "Unknown" stands for not indicated by the authors.

solvers for the different problem formulations are grouped in Table 9, Table 11, Table 12, andTable 13. -Statistical inference based approaches: Hierarchical statistical models are used to model the data generation process and prior distributions are selected to capture the low-rank and sparse properties of the respective terms. The joint distribution involving the observations, unknown variables and hyperparameters can be determined from the priors and conditional distributions.

Posterior distributions of the unknowns are approximated using Bayesian inference approaches. Representative algorithms in this category can be found in [71][17][5][385].

For the SVD algorithms, approximated SVD solutions exist to avoid full SVD such as partial SVD [188], linear time SVD [357], limited memory SVD [195], symmetric low-rank product-Gauss-Newton [196], and Block Lanczos with Warm Start (BLWS) [191].

Solvers for PCP Complexity

Linearized solvers Linearized Augmented Lagrangian Method (LALM) Oiter (mnmin(mn)), Opre=Unknown, Oconv =Unknown Yang and Yuan [352] Linearized Alternating Direction Method (LADM)

Oiter (mnmin(mn)), Opre=Unknown, Oconv =Unknown Yang and Yuan [352] LADM with Adaptive Penalty (LADMAP 5 )

Oiter (rmn), Opre(1/ǫ), Oconv =Unknown Lin et al. 2011 [190] Accelerated version Linearized Symmetric Alternating Direction Method (LSADM 5 )

Oiter =Unknown, Opre(1/ǫ), Oconv =Unknown (ALM) Ma [205], Goldfarb et al. [90] Fast Linearized Symmetric Alternating Direction Method (Fast-LSADM 5 )

Oiter =Unknown, Opre( 1/ǫ), Oconv =Unknown (FALM) Ma [205] Oiter (mk 2 ) with k ≪ n, Opre=Unknown, Oconv =Unknown Abdel-Hakim and El-Saban [START_REF] Abdel-Hakim | FRPCA: fast robust principal component analysis[END_REF] Full SVD Inexact Fast Robust Principal Component Analysis (IFRPCA)

Oiter (mk 2 ) with k ≪ n, Opre=Unknown, Oconv =Unknown Abdel-Hakim and El-Saban [START_REF] Abdel-Hakim | FRPCA: fast robust principal component analysis[END_REF] Partial SVD Matrix Tri-Factorization (MTF)

Oiter (n 3 + (r 3 + r 2 n + mn 2 + rn 2 )) Liu et al. [199] Opre=Unknown, Oconv =Unknown Fast Tri-Factorization(FTF)

Oiter (r Oiter =Unknown, Opre=Unknown, Oconv =Unknown Aybat et al. [12] First-order Augmented Lagrangian algorithm for Composite norms (FALC9 )

Oiter =Unknown, Oconv =Unknown, Opre(1/ǫ) Aybat et al. [START_REF] Aybat | An augmented Lagrangian method for conic convex programming[END_REF] Augmented Lagragian method for Conic Convex (ALCC 9 )

Oiter =Unknown, Oconv =Unknown, Opre(log(1/ǫ)) Aybat et al. [START_REF] Aybat | A unified approach for minimizing composite norms[END_REF] Partially Smooth Proximal Gradient (PSPG 9 )

Oiter =Unknown, Oconv =Unknown, Opre(1/ǫ) Aybat et al. [START_REF] Aybat | Efficient algorithms for robust and stable principal component pursuit[END_REF] Alternating Direction Method -Increasing Penalty (ADMIP 9 )

Oiter =Unknown, Oconv =Unknown, Opre=Unknown Aybat et al. [16] Inexact Alternating Minimization -Matrix Manifolds (IAM-MM)

Oiter (rmn),Oconv=Unknown, Opre=Unknown (R2PCP Oiter (r(m + n) + mn) Robert and Cassela [256] Oconv =Unknown, Opre=Unknown Variational Bayesian Inference (VB)

Oiter =Unknown, Oconv =Unknow, Opre=Unknow Beal [19] Variational Bayesian RPCA (VB-RPCA 16 )

Babacan et al. [17] Approximate Bayesian Inference (AB)

Oiter (min(n 3 , r 3 ) + min(m 3 , r 3 )) Beal [19] Oconv =Unknown, Opre=Unknown Approximated RPCA (A-RPCA) (GoDec17 ) Zhou and Tao [389] Naive GoDec Linear convergence Zhou and Tao [389] Fast Godec via Bilateral Random Projection Linear convergence Zhou and Tao [389] Table 11 Solvers for RPCA (excepted PCP): An overview of their complexity per iteration at running time O iter , their complexity Opre to reach an accuracy of ǫ precision and their convergence rate Oconv for T iterations. "Unknown" stands for not indicated by the authors.

Methods

Solvers Complexity

Robust Non-negative Matrix Factorization (RNMF) Manhattan Non-negative Matrix Factorization (MahNMF Table 13 Solvers for ST and LRM: An overview of their complexity per iteration at running time O iter , their complexity Opre to reach an accuracy of ǫ precision and their convergence rate Oconv for T iterations. "Unknown" stands for not indicated by the authors.

Adequation for the background/foreground separation

For each problem formulations, we investigated its adequation with the application of background/foreground separation in their corresponding section in terms of following criteria: (1) its robustness to noise, (2) its spatial and temporal constraints, (3) the existence of an incremental version, and (4) the existence of a real-time implementation. The following observations can be made:

1. Robustness to noise: Noise is due to a poor quality image source such as images acquired by a web cam or images after compression. It affects the entries of the matrix A. In each problem formulation, assumptions are made to assure the exact recovery of the decomposition. PCP assumed that all entries of the matrix to be recovered are exactly known via the observation and that the distribution of corruption should be sparse and random enough without noise. These assumptions are rarely verified in the case in real applications because only a fraction of entries of the matrix can be observed and the observation can be corrupted by both impulsive and Gaussian noise. BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF], spatio-temporal IRLS [111], DECOLOR [391] and MBRMF [366] address both the spatial and temporal constraints. The different strategies used to take into account the spatial and/or temporal coherence can be classified as follows:

-For the regularization based approaches, the main strategies to take into account the spatial coherence consist of using a mixed norm (||.|| 2,1 [306][109][113][112] [111]) on the matrices L and/or S, using a structured sparsity norm [197] on the matrix S, and adding a term in the minimization problem such as a Total Variation penalty [112][111][334][108][43] or a gradient [113][112][111][333] on the matrix S. For the temporal coherence, optical flow is used as in the RPCA framework ([87][138]) and in RMC framework (MAMC [351], RMAMC [351]). Thus, the motion information can be used in several ways as follows :

(1) an adaptive λ [87][197] which is a function of the motion consistency to ensure that all the changes caused by the foreground motion will be entirely transfered to the matrix S, and (2) a weighting matrix W

[351] which is constructed from the optical flow to suppress slowly-moving objects, to enforce the recovery of the background that appears at only a few frames and to eliminate the influence of light conditions, camouflages, and dynamic backgrounds. Thus, the minimization problem expressed in Equation 2 can be extended for stable explicit decomposition to: min

L,S,E λ 1 f low (L) + λ 2 fsparse(Π(S)) + λ 3 f noise (E) Decomposition + λ 4 f back (L) + λ 5 f f ore (S) Application subj C 3 (12)
where λ 4 and λ 5 are regularization parameters. f back (L) and f f ore (S) are loss functions that allow the minimization to take into account the characteristics of the background and the foreground, respectively. f back (L) can be a mixed norm. f f ore (S) can be the gradient or the Total Variation on S. The function Π() allows to add a confidence map [236][292] on S. C 3 contains the constraints which are only on the recovery such as A = L+S+E,or both on the recovery and the spatial/temporal aspects such as A= W • (L + S + E) [351][365][292]. Note that the first part of Equation 12with f low (L), fsparse(S) and f noise (E) concerns mainly the decomposition into low-rank plus additive matrices, and the second part with f back (L) and f f ore (S) concerns mainly the application to background/foreground separation. Thus, the minimization problem can be formulated as follows:

min

L,S,E λ 1 ||L|| p1 norm1 + λ 2 ||Π(S)|| p2 norm2 + λ 3 ||E|| p3 norm3 Decomposition + λ 4 ||L|| l2,1 + δ 1 ||grad(S)|| l1 + δ 2 T V (S) + δ 3 Ω(S) Application subj ||A -L -W • S -E|| = 0 (13)
where δ 1 , δ 2 and δ 3 are regularization parameters. norm 2 is usually taken to force spatial homogeneous fitting in the matrix S, that is for example the norm l 2,1 with p 2 = 1 [306][109][113][112][111]. ||grad(S)|| 1 , T V (S) and Ω(S) are a gradient [113][112][111][333], a total variation [112][111][334][108] [START_REF] Cao | Total variation regularized RPCA for irregularly moving object detection under dynamic background[END_REF] and a structured norm [197] applied on the matrix S, respectively. Furthermore, S can be processed with a linear operator Π() that weights its entries according to their confidence of corresponding to a moving object such that the most probable elements are unchanged and the least are set to zero [236][292]. Note that the term λ 4 ||L|| l2,1 ensures the recovered L has exact zero columns corresponding to the outliers. Morevover, a weighting matrix W [351][365][292] can be used in the constraints C 3 to enforce the recovery of the background that appears at only a few frames and to eliminate the influence of light conditions, camouflages, and dynamic backgrounds.

-For the statistical inference based approaches, Markov Random Fields (MRF) are used to extract temporally and spatially localized moving objects as in BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF], DECOLOR [391] and MBRMF [366]. Statistical total variations can also be used as in the approach based on smoothness and arbitrariness constraints [108]. [4]). Thus, the decomposition can be written as follows:

A t = L t + S t + E t ( 14
)
where t is the indice for the time. L t , S t , E t are determined from L t-1 , S t-1 , E t-1 and the current observation.

4.

Real-time implementations: As background/foreground separation needs to be achieved in real-time, several strategies have been developed and are generaly based on submatrices computation [240] or GPU implementations [8][9]. Real-time implementations can be found for PCP [8][9][240][106] and for SPCP [209] 5. Strategies: Differents strategies can be used to applied DLAM for background/foreground separation. For example, Gao et al. [START_REF] Gao | Block-sparse RPCA for consistent foreground detection[END_REF] developed a two-pass RPCA process for consistent foreground detection. For objects or people which remain immobile for a certain period of time. Tepper et al. [308] proposed a method which detects foreground objects at different timescales, by exploiting the theoretical and practical properties of RPCA.

All these key challenges need to be addressed in the different problem formumations based on the decomposition into low-rank plus sparse matrices to be applied adequately to background modeling and foreground detection in video taken by a static camera.

Sparse decompositions

Sparse decompositions are similar to low-rank decompositions except that the first matrix is considered to be sparse instead of low-rank. Sparse decompositions are achieved in the different following problem formulations:

-Sparse Dictionary learning: Sparse dictionary learning (DL) builds data representation by decomposing each datum into a linear combination of a few components selected from a dictionary of basic elements, called atoms [START_REF] Chen | Compressed sensing and dictionary learning[END_REF]. Sparse dictionary learning is also called sparse coding in the literature [297][230]. Thus, the observation matrix is decomposed as follows:

A = X + N = Dα + N ( 15
)
where A is the matrix which contains the observations, X is a sparse noiseless matrix and N is the noise matrix. X is the product between D which is a dictionary, and α which is a sparse vector. Thus, the assumption of the sparse decomposition is that the observed image is an approximated linear combination over a dictionary D and a vector coefficients α. In order to recover the noiseless image (background), the decomposition problem seeks to the following minimization problems:

min D,α ||A -X|| 2 l2 + ||Dα -X|| 2 l2 + ||α|| 0 (16)
The first term minimizes the error between the recovered image and the observed version. The second term ensures that the denoised image is an approximated linear combination over the dictionary D and coefficients α. Finally, the third term determines the degree of sparsity of the coefficients, in fact 0 counts the null coefficients. Thus, the recovered image is represented with the smallest possible number of vectors from the dictionary. The minimization problem is solved iteratively in three steps. First, a solver such as a matching pursuit type algorithm is used to estimate the coefficients of the linear decomposition of the denoised image over the dictionary Second, the dictionary is updated. And finally, the last step updates the denoised image X. Applied to background/foreground separation, a dictionary learning method considers that 1) The background has a sparse linear representation over a learned dictionary, and 2) the foreground is sparse in the sense that majority pixels of the frame belong to the background. Learning the dictionary is a key step to the success of background modeling. The different approaches developed in the literature differs from 1) the algorithm used to learn the dictionary (K-SVD [3][63][61][62][288][139][267], RDL [380], MOD [388]), BPFA [388]), 2) the decomposition (Two terms [63][61][380], three terms [START_REF] David | Sparse coding and gaussian modeling of coefficients average for background subtraction[END_REF]), 3) the minimization problem with a different norm on the sparse error (l- [288], IRLS [203] and Least Angle Regression (LARS) [79][339]. Furthermore, online dictionary learning algorithms are developed as in Lu et al. [203], and Zhang et al. [377] with Symmetric Positive Definite (SPD) matrices.

-Sparse Linear Approximation/Regression: This formulation problem is similar to sparse dictionnary learning and leads to the same decomposition. First, Dikmen et al. [68][69][67] refer to linear approximation of the sparse error estimation, and basis selection (i.e the dictionary). This method viewed foreground objects as sparse corruption signals and estimated them by the sparse recovering method. Second, other authors [347][346][131] refer to sparse outlier estimation in a linear regression model regarding foreground objects as outliers and consider that the observation error is composed of foreground outlier and background noise. Thus, the foreground detection task has been converted into a outlier estimation problem.

-Compressive sensing: The CS theory states that a signal can be reconstructed from a small number of measurements with high probability, provided that the signal is sparse in the spatial domain or some transform domains [START_REF] Chen | Compressed sensing and dictionary learning[END_REF]. Assume that a signal X can be represented as X = Ψ Θ , where Ψ denotes a basis and Θ is the coefficients corresponding to the basis. The signal is said to be k-sparse if all other elements in Θ vanish except for nonzero coefficients. According to CS, for a sparse signal, compressive measurements can be collected by the following random projections:

A = ΦX + N (17)
where Φ ∈ R m×n is the measurement matrix with m ≤ n, A contains m measurements, and N is the measurement noise. Specifically, a high dimensional X vector is converted into a much lower dimensional measurement vector A.

Moreover, the compressive measurements in contain almost all the information of the sparse vector X. This means that CS works with data of significantly lower dimension so as to achieve computation efficiency as well as accuracy. In order to recover the noiseless signal, the decomposition problem seeks to the following minimization problems: min

D,α ||Θ|| l0 + 1 2 ||A -ΦX|| 2 (18)
Because Equation 18is an NP hard problem, the sparse solution can be obtained by replacing the nom l 0 by the norm l 1 as done in [45] [44][121][176][323][327]. Thus, the background/foreground separation problem can be viewed as a sparse approximation problem where convex optimization and greedy methods can be applied. It is not necessary to learn the background itself to detect the changes and the foreground objects which can be directly detected on the compressive samples. Hence, no foreground reconstruction is done until a detection is made to save computation. The different approaches developed in the literature differs mainly by the minimization problem (l 1 -l 1 minimization [219][218]) and the solvers. Examples of solvers include Basis Pursuit (BP) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][64], Basis Pursuit Denoising (BPDN) [39][45], Orthogonal Matching Pursuit (OMP) [START_REF] Davies | The effect of recovery algorithms on compressive sensing background subtraction[END_REF], Stagewise OMP (StOMP) [START_REF] Donoho | Sparse solution of under-determined linear equations by stagewise orthogonal matching pursuit[END_REF][344], Lattice Matching Pursuit (LaMP) [START_REF] Cevher | Sparse signal recovery using markov random fields[END_REF],

Compressive sampling matching pursuit (CoSaMP) [225][121], and Gradient Projection for Sparse Reconstruction (GPSR) [323]. Furthermore, structured sparsity [135][136][133][44] can be used to exploit a priori spatial information on coefficient structure in addition to signal sparsity as the foreground objects are usually not only sparse but also clustered in a distinct way. Dynamic Group Sparsity (DGS) [134][338][255][327] can also be used to exploit both temporal and spatial information. In another way, an adaptive algorithms called Adaptive Rate Compressive Sensing (ARCS) [328] [328] allow to choose the number of measurements so as to limit the data rate of the sensor while simultaneously maintaining enough information such that to be able to robustly detect the foreground objects.

As we defined the unified view DLAM, it is possible to define a unified view of the sparse decompositions called Decomposition in sparse plus additive matrices (DSAM) but it is out of the scope of this paper.

Mixed decompositions

Mixed decompositions stand at the intersection of the previous problem formulations. There are two main approaches in the literature:

-RPCA-CS: Waters et al. [329][330] proposed to recover the entries of a matrix A in terms of a low-rank matrix L and sparse matrix S from a small set of compressive measurements y = A(L+S) where A is an underdetermined linear operator. The optimization problem that unites the above two problem classes above is: Jiang et al. [152] extented this model to deal with the joint reconstruction of multiple color components. Jiang et al. [153] improved this model by adding low latency. In another way, Yang et al. [350] developed an online algorithm in which the background is learned adaptively as the compressive measurements are processed.

min ||y -A(L + S)|| l2 subj rank(L) < r, ||vec(S)|| l0 < K (19)
-Sparse Dictionary Learning-CS: Huang et al. [139] proposed an algorithm of moving object detection via the sparse representation and learned dictionary. First, compress image in order to reduce data redundancy and bandwidth. Then, data dictionary with CS measurement values and sparse basis is initialized, trained and updated through the K-SVD. Finaly, moving object detection is achieved via PCP. In another approach, Jiang et.

[154] used a spatial-temporal image patch (bricks) as atomic unit for sparse dictionary representation. Furthermore, Random Projection emerged from Compressive Sensing theory is used to reduce the dimension of the bricks so as to speed up the algorithm.

Motivations and Contributions

Since the works of Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] and Chandrasekaran et al. [START_REF] Chandrasekaran | Ranksparsity incoherence for matrix decomposition[END_REF], the last five years witnessed very significant publications on problem formulations based on the decomposition into low-rank plus additive matrices, and applications in computer vision generate new developments as developed in the handbook [START_REF] Bouwmans | Handbook on robust low-rank and sparse matrix decomposition: Applications in image and video processing[END_REF]. Furthermore, the different robust problem formulations based on the decomposition into lowrank plus additive matrices often outperform state-of-the-art methods in several computer vision applications [START_REF] Candes | Discussion of latent variable graphical model selection via convex optimization[END_REF]. Indeed, as this decomposition is nonparametric and does not make many assumptions, it is widely applicable to a large scale of problems ranging from:

-Latent variable model selection: Chandrasekaran et al. [START_REF] Chandrasekaran | Latent variable graphical model selection via convex optimization[END_REF] proposed to discover the number of latent components, and to learn a statistical model over the entire collection of variables by only observing samples of a subset of a collection of random variables. The geometric properties of the decompostion of low-rank plus sparse matrices play an important role in this approach [47][206].

-Image processing: Sometimes, it is needed to separate information from noise or outliers in image processing. RPCA framework was applied with success in image analysis [392] such as image denoising [START_REF] Gu | Weighted nuclear norm minimization with application to image denoising[END_REF], image composition [START_REF] Bhardwaj | Robust PCA-based solution to image composition using augmented lagrange multiplier (alm)[END_REF] [START_REF] Chen | Incremental low-rank and sparse decomposition for compressing videos captured by fixed cameras[END_REF], key frame extraction [START_REF] Dang | RPCA-KFE: key frame extraction for consumer video based robust principal component analysis[END_REF], hyperspectral video processing [START_REF] Gerhart | Convex optimization techniques and their application in hyperspectral video processing[END_REF], video restoration [149], video stabilization [START_REF] Ebadi | Approximated RPCA for fast and efficient recovery of corrupted and linearly correlated images and video frames[END_REF], and in background and foreground separation [6][236][241].

-3D Computer Vision: Structure from Motion (SfM) refers to the process of automatically generate a 3D structure of an object by its tracked 2D image frames. Practically, the goal is to recover both 3D structure, namely 3D coordinates of scene points, and motion parameters, namely attitude (rotation) and position of the cameras, starting from image point correspondences. Then, finding the full 3D reconstruction of this object can be posed as a low-rank matrix recovery problem [192][11][337].

In this context, the aim of this survey is then to provide a first complete overview of all the decomposition in low-rank plus additive matrices for (1) novices who could be students or engineers beginning in the field of computer vision, ( 2) experts as we put forward the recent advances that need to be improved, and (3) reviewers to evaluate papers in journals, conferences, and workshop such as RSL-CV 2015 29 . So, this survey is intended to be a reference for researchers and developers in industries, as well as graduate students, interested in robust decomposition applied to computer vision.

Here, we decide to focus on the application of background and foreground separation due to the following reasons:

1. This application witnessed very numerous papers (more than 190) since 2009. 2. Background/foreground separation is the most representative and demanding application as it needs to take into account both spatial and temporal constraints with incremental and real-time constraints [33][32].

However, no algorithm today seems to emerge and to be able to simultaneously address all the key challenges that accompany real-world videos. This is due, in part, to the absence of a rigorous quantitative evaluation with large-scale datasets with accurate ground truth providing a balanced coverage of the range of challenges present in the real world. Indeed, in the first publications, the authors usually compared qualitatively their method to RSL [310] or PCP [START_REF] Candes | Robust principal component analysis? International[END_REF]. Recent quantitative evaluations in foreground detection using the performance metrics have been made but they are limited to one algorithm [326][355][110][109]. In a more recent work, Guyon et al. [114] compared five algorithms RSL [310], RPCA-PCP solved via EALM [188], RPCA-PCP solved via IALM [188], QPCP [20] and BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF] with the Wallflower dataset [312], the I2R dataset [177] and Shah dataset [279]. Experimental results show that BRPCA that address spatial and temporal constraints outperforms the other methods. However, this evaluation is limited to five methods and it is not made on large datasets that present a coverage of the range of challenges. A similar study made by Rueda et al. [265] compared RPCA-PCP solved via EALM, BRPCA and GoDec. The authors also concluded that the BRPCA offers the best results in dynamic and static scenes by exploiting the existing correlation between frames of the video sequence using Markov dependencies. In a more complete survey, Bouwmans and Zahzah [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF] evaluated ten RPCA-PCP algorithms on the BMC dataset but this evaluation is limited to the framework of RPCA solved via PCP.

Moreover, we believe that we are living in a key transition in the field of background subtraction as we are progressively migrating from the conventional statistical models as MOG [298] [START_REF] Bouwmans | Background modeling using mixture of Gaussians for foreground detection -a survey[END_REF], KDE [80][27] and naive subspace learning models [START_REF] Bouwmans | Subspace learning for background modeling: A survey[END_REF] to models based on robust decomposition into low-rank plus additive matrices (RPCA, RNMF, RMC, RSR, RST, LRM) which can achieve at least the same performance in terms of precision than the conventional statistical models [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF]. Thus, the aim of this survey is to review and evaluate the robust decomposition into low-rank plus additive matrices for the application of background/foreground separation. For this, it reviews all the models since the first works of Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] and Chandrasekaran et al. [48] to the recent ones. By reviewing both existing and new ideas, this survey gives a complete overview of the decompositions, solvers, and applications related to background/foreground separation. Moreover, an accompanying website called the DLAM Website 30 is provided. It allows the reader to have a quick access to the main resources, and codes in the field. Finally, with this survey, we aim to bring a one-stop solution, i.e., access to a number of different decompositions, solvers, implementations and benchmarking techniques in a single paper. Considering all of this, we present a comprehensive review of different methods based on decomposition into low-rank plus additive matrices for testing and ranking existing algorithms for foreground detection. Contributions of this paper can be summarized as follows:

-A unified view of the decomposition into low-rank plus additive matrices: After a preliminary overview on the different robust problem formulations in Section 1.1, we provided in Section 1.2 a unified view of the different decompositions into low-rank plus additive matrices. Figure 3 shows an overview of this unified view.

-A review regarding different decomposition methods in low-rank plus additive matrices: RPCA models are reviewed in Section 2. For each method, we investigate how they are solved, and if incremental and real-time versions are available for foreground detection. Furthermore, their advantages and drawbacks are discussed in the case of outliers due to dynamic backgrounds or illumination changes. In the same manner, we review the RNMF models in Section 3, RMC models in Section 4, RSR models in Section 5, and the robust subspace tracking models in Section 6. Finally, robust low-rank minimization models are reviewed in Section 7.

-A systematic evaluation and comparative analysis: We compare and evaluate different decomposition methods in low-rank and additive matrices on a large-scale dataset in Section 8. This dataset is the Background Models Challenge (BMC 2012) dataset 31 [315] and we used the provided quantitative evaluation framework which allows us to do a fair and complete comparison.

The rest of this paper is organized as follows. Firstly, we review each original method in its section (Section 2 to Section 7). For each method, we investigate how they are solved, and if incremental and real-time versions are available for background/foreground separation. Then, the performance evaluation using quantitative metrics over the BMC dataset is given in Section 8. Finally, we conclude with promising research directions in Section 10.

2 Robust Principal Component Analysis

RPCA via Principal Component Pursuit

RPCA via PCP proposed by Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] in 2009 is currently the most investigated method. In the following sub-sections, we reviewed this method and all these modifications in terms of decomposition, solvers, incremental algorithms and real time implementations. Table 14 and Table 15 show an overview of the Principal Component Pursuit methods and their key characteristics.

Table 14 Principal Component Pursuit: A Complete Overview (Part 1). The first column indicates the concerned category and the second column the name of each method. Their corresponding acronym is indicated in the first parenthesis. The third column gives the name of the authors and the date of the related publication.

Table 15 Principal Component Pursuit: A Complete Overview (Part 2). The first column indicates the concerned category and the second column the name of each method. Their corresponding acronym is indicated in the first parenthesis. The third column gives the name of the authors and the date of the related publication.

Principal Component Pursuit

Candes et al. [40][335] proposed a convex optimization to address the robust PCA problem. The observation matrix A is assumed represented as:

A = L + S ( 20
)
where L is a low-rank matrix and S must be sparse matrix with a small fraction of nonzero entries. The straightforward formulation is to use l 0 -norm to minimize the energy function:

min L,S rank(L) + λ||S|| l0 subj A -L -S = 0 (21)
where λ > 0 is an arbitrary balanced parameter. But this problem is N P -hard, typical solution might involve a search with combinatorial complexity. This research seeks to solve for L with the following optimization problem: min

L,S ||L|| * + λ||S|| l1 subj A -L -S = 0 ( 22
)
where ||.|| * and ||.|| l1 are the nuclear norm (which is the l 1 -norm of singular value) and l 1 -norm, respectively, and λ > 0 is an arbitrary balanced parameter. Usually,

λ = 1 √ max(m,n)
. Under these minimal assumptions, this approach called Principal Component Pursuit (PCP) solution perfectly recovers the low-rank and the sparse matrices. Candes et al. [393] showed results on face images and background modeling that demonstrated encouraging performance. The low-rank minimization concerning L offers a suitable framework for background modeling due to the correlation between frames. So, minimizing L and S implies that the background is approximated by a low-rank subspace that can gradually change over time, while the moving foreground objects constitute the correlated sparse outliers which are contained in S.

To obtain the foreground mask, S needs to be thresholded. The threshold is determined experimentally. rank(L) influences the number of modes of the background that can be represented by L: If rank(L) is to high, the model will incorporate the moving objects in its representation; if the rank(L) is to low, the model tends to be uni-modal and then the multi-modality which appears in dynamic backgrounds will be not captured. The quality of the background/foreground separation is directly related to the assumption of the low-rank and sparsity of the background and foreground, respectively. In this case, the best separation is then obtained only when the optimization algorithm has converged.

Essentially, the nuclear-norm term corresponds to low-frequency components along the temporal while the l 1 norm describes the high-frequency components. However, the low-frequency components can leak into extracted background images for areas that are dominated by moving objects. The leakage as ghost artifacts which appear in extracted background cannot be well handled by adjusting the weights between the two regularization parameters. Practically, RPCA-PCP present other several limitations developed in the Section 1.1.1 and an overview of their solutions is given in the following sections.

Algorithms for solving PCP

Several algorithms called solvers have been proposed for solving PCP. An overview of these solvers as well as their complexity (when they are available) are grouped in Table 9 andTable 10. For an m × n input matrix A with estimated rank r, the complexity per iteration at running time is formulated as O iter (f iter (m, n, r, ..)) where f iter () is a function. The complexity to reach an accuracy of ǫ precision (ǫ-optimal solution) is formulated as Opre(fpre(ǫ)) where fpre() is a function. The convergence rate is formulated as Oconv(fconv(T )), where fconv() is a function of T which is the number of iterations. All these algorithms require solving the following type of subproblem in each iteration: min L,S η||L|| p1 norm1 + λ||S|| p2 norm2 [START_REF] Bian | Sparse and low-rank modeling on high dimensional data: A geometric perspective[END_REF] The above problem can have a closed form solution or not following the application. So, several solvers can be found in the the litterature:

-Basic solvers: When the problem is supposed to have a closed form solution, PCP can be reformulated as a semidefinite program and then be solved by standard interior point methods [START_REF] Chandrasekharan | Rank-sparsity incoherence for matrix decomposition[END_REF]. However, interior point methods have difficulty in handling large matrices because the complexity of computing the

step direction is O((mnmin(m, n)) 2 )
, where m×n is the size of the data matrix. If m = n, then the complexity is O(n 6 ). So the generic interior point solvers are too limited for many real applications where the number of data are very large.

To overcome the scalability issue, only the first-order information can be used. Cai et al. [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] showed that this technique, called Singular Value Thresholding (SVT), can be used to minimize the nuclear norm for matrix completion. As the matrix recovery problem in Equation ( 22) needs minimizing a combination of both the l 1 -norm and the nuclear norm, Wright et al.

[335] adopted a iterative thresholding technique (IT) to solve it and obtained similar convergence and scalability properties than interior point methods. However, the iterative thresholding scheme converges extremely slowly with Opre = L/ǫ where L is the Lipschitz constant of the gradient of the objective function. To alleviate this slow convergence, Lin et al. [189] proposed two algorithms: the accelerated proximal gradient (APG) algorithm and the gradient-ascent algorithm applied to the dual of the problem in Equation ( 22). However, these algorithms are all the same to slow for real application. More recently, Lin et al. [188] proposed two algorithms based on augmented Lagrange multipliers (ALM). The first algorithm is called Exact ALM (EALM) method that has a Q-linear convergence speed, while the APG is in theory only sub-linear. The second algorithm is an improvement of the EALM that is called Inexact ALM (IALM) method, which converges practically as fast as the EALM, but the required number of partial SVDs is significantly less. The IALM is at least five times faster than APG, and its precision is also higher [188]. However, the direct application of ALM treats the Equation ( 22) as a generic minimization problem and ignores its separable structure emerging in both the objective function and the constraint [188]. Hence, the variables S and L are minimized simultaneously. Yuan and Yang [368] [START_REF] Chai | An improved ADM algorithm for RPCA optimization problem[END_REF] is an improved alternating direction method (IADM) algorithm with a block based SVD approach. Experimental results [START_REF] Chai | An improved ADM algorithm for RPCA optimization problem[END_REF] on the I2R dataset [177] show that IADM outperforms SVT [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], APG [189], IALM [188] and ADM [368] with less computation time.

-Linearized solvers: When the resulting subproblems do not have closed-form solutions, Yang and Yuan [352] proposed to linearize these subproblems such that closed-form solutions of these linearized subproblems can be easily derived. Global convergence of these Linearized ALM (LALM) and ADM (LADM) algorithms are established under standard assumptions. Recently, Lin et al. [190] improved the convergence for the Linearized Alternating Direction Method with an Adaptive Penalty (LADMAP). They proved the global convergence of LADM and applied it to solve Low-Rank Representation (LRR). Furthermore, the fast version of LADMAP reduces the complexity O(mnmin(m, n)) of the original LADM based method to O(rmn), where r is the rank of the matrix to recover, which is supposed to be smaller than m and n. In a similar way, Ma [205] and Goldfarb et al. [START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex function[END_REF] proposed a Linearized Symmetric Alternating Direction Method (LSADM) for minimizing the sum of two convex functions. This method requires at most O(1/ǫ) iterations to obtain an ǫ-optimal solution, and its fast version called Fast-LSADM requires at most O(1/ √ ǫ) with little change in the computational effort required at each iteration.

-Fast solvers: All the previous solvers require computing SVDs for some matrices, resulting in O(mnmin(m, n)) complexity. Although partial SVDs are used to reduce the complexity to O(rmn) such a complexity is still high for large data sets. Therefore, recent researches focus on the reduction of the complexity by avoiding computation of SVD. Shen et al. [281] presented a method where the low-rank matrix is decomposed in a product of two low-rank matrices and then minimized over the two matrices alternatively. Although, they do not require nuclear norm minimization and so the computation of SVD, the convergence of the algorithm is not guaranteed as the problem is non-convex. Furthermore, both the matrix-matrix multiplication and QR decomposition based rank estimation technique require O(rmn) complexity. So, this method does not essentially reduce the complexity. In another way, Mu et al. [221] reduced the problem scale by random projections (linear or bilinear projection) but different random projections may lead to radically different results. Furthermore, additional constraints to the problem slow down the convergence. The complexity of this method is O(pmn) where p ×m is the size of the random projection matrix with p ≪ m, p ≪ n and p > r. So, this method is still nor linear complexity with respect to the matrix size. Its convergence needs more iterations than IALM but it requires less computation time. In another way, Liu et al. [194][193] proposed an algorithm, called l 1 -filtering for exactly solving PCP with an O(r 2 (m + n)) complexity. This method is a truly linear cost method to solve PCP problem when the date size is very large while the target rank is small. Moreover, l 1 -filtering is highly parallelizable. It is the first algorithm that can exactly solve a nuclear norm minimization problem in linear time.

Numerical experiments [194][193] show the great performance of l 1 -filtering in speed compared to the previous algorithms for solving PCP -2D solvers: Sun et al. [301][303] developed an iterative algorithm for robust 2D-PCA via alternating optimization which learns the projection matrices by bi-directional decomposition. To further speed up the iteration, Sun et al. [301][303] proposed an alternating greedy approach or l 0 -norm regularization, minimizing over the low-dimensional feature matrix and the sparse error matrix.

Algorithms for incremental PCP

PCP is an offline method which treats each image frame as a column vector of the matrix A. In real-time application such as foreground detection, it would be more useful to estimate the sparse matrix in an incremental way quickly as new frame comes rather than in batch way. Furthermore, the sparsity structure may change slowly or in a correlated way, which may result in a low-rank sparse matrix. In this case, PCP assumption will not get satisfied and S can't be separated from L. Moreover, the principal directions can change over time. So, the rank of the matrix L will keep increasing over time thus making PCP infeasible to do after time. This last issue can be solved by not using all frames of the sequences. So, several algorithms for incremental algorithms are available in literature.

A t = L t + S t = [UI] [x t S t ] t (24)
where Lois and Vawani [200]. However, ReProCS requires the support x t to be fixed and quite small for a given support size S t , but this does often not hold. So, ReProCS could not handle large outliers support sizes.

x t = U T L t and
2. Support-Predicted Modified-CS : Qiu and Vaswani [244] address the large outliers support sizes by using time-correlation of the outliers. This method called Support-Predicted Modified-CS RR-PCP [244] and Support-Predicted Modified-CS [243] is also an incremental algorithm and outperforms the Re-ProCS. However, this algorithm is only adapted for specific situation where there are only one or two moving objects that remain in scene. But, this is not applicable to real videos where multiple and time-varying number of objects can enter of leave the scene. Moreover, it requires knowledge of foreground motion.

Automated Recursive Projected CS (A-ReProCS):

To address the limitation of the Support-Predicted Modified-CS, Qiu and Vaswani [245] proposed a method called automated Recursive Projected CS (A-ReProCS) that ensures robustness when there are many nonzero foreground pixels, that is, there are many moving objects or large moving objects. Furthermore, A-ReProCS outperforms the previous incremental algorithms when foreground pixels are correlated spatially or temporally and when the foreground intensity is quite similar to the background one.

4. ReProCS with cluster-PCA (ReProCS-cPCA): To address the fact that the structure that we require is that L t is in a low dimensional subspace and the eigenvalues of its covariance matrix are "clustered", Qiu and Vaswani [247][249] introduced a Recursive Projected Compressive Sensing with cluster-PCA (ReProCS-cPCA). Under mild assumptions, ReProCS-cPCA with high probability can exactly recover the support set of S t at all times. Furthermore, the reconstruction errors of both S t and L t are upper bounded by a timeinvariant and small value. Despite the efforts to reduce the time complexity, the corresponding algorithms have prohibitive computational time for real application such as foreground detection. The main computation in PCP is the singular value decomposition of the large matrix A. Instead of computing a large SVD on the CPU, Anderson et al. [8][9] proposed an implementation of the communication-avoiding QR (CAQR) factorization that runs entirely on the GPU. This implementation achieved 30× speedup over the implementation on CPU using Intel Math Kernel Library (Intel MKL).

In another way, Pope et al. [240] proposed a variety of methods that significantly reduce the computational time of the ALM algorithm. These methods can be classified as follows:

-Reducing the computation time of SVD: The computation of the SVD is reduced using the Power method [240] that enables to compute the singular values in a sequential manner, and to stop the procedure when a singular value is smaller than a threshold. The use of the Power method by itself results in 4.32× lower runtime. Furthermore, the gain is improved by a factor of 2.02× speedup if the procedure is stopped when the singular value is smaller than the threshold. If the rank of L is fixed and the Power SVD is stopped when the number of singular value is equal to rank(L), the additional speedup is 17.35.

-Seeding the PCP algorithm: PCP operates on matrices consisting of blocks of contiguous frames acquired with a fixed camera. So, the low-rank matrix does not change significantly from one block to the next. Thus, Pope et al.

[240] use the low-rank component obtained by the ALM algorithm from the previous block as a starting point for the next block. This method allows an additional speedup of 7.73. -Partioning into subproblems: Pope et al. [240] proposed to partition the matrix A into P smaller submatrices. The idea is to combine the solutions of the P corresponding PCP subproblems to recover the solution of the full matrix A at lower computational complexity.

In this way, Pope et al. [240] demonstrated that the PCP algorithm can be in fact suitable for real-time foreground/background separation for video-surveillance application using off-theshelf hardware.

In a similar manner, Guo et al. [106] proposed a low-rank matrix recovery scheme, which splits the original RPCA into two small ones: a low-rank submatrix recovery and a low-rank submatrix reconstruction problems. This method showed a speedup of the ALM algorithm by more than 365 times compared to a C implementation with less requirement of both time and space. In addition, this method significantly cuts the computational load for decomposing the remaining frames.

2.1.5 Methods for finding the optimal PCP solution PCP recovers the true underlying low-rank matrix when a large portion of the measured matrix is either missing or arbitrarily corrupted. However, in the absence of a true underlying signal L and the deviation S, it is not clear how to choose a value of λ that produces a good approximation of the given data A for a given application.A typical approach would involve some cross-validation step to select λ to maximize the final results of the application. The issue with crossvalidation in this situation is that the best model is selected indirectly in terms of the final results, which can depend in unexpected ways on later stages in the data processing chain of the application. Instead, Ramirez and Sapiro [253][254] addressed this issue via the Minimum Description Length (MDL) principle [117] and so proposed a MDL-based low-rank model selection. The principle is to select the best low-rank approximation by means of a direct measure on the intrinsic ability of the resulting model to capture the desired regularity from the data. To obtain the family of models M corresponding to all possible low-rank approximation of A, Ramirez and Sapiro [253][254] applied the RPCA decomposition for a decreasing sequence of values of λ, {λ t : t = 1, 2, 3, ...} obtaining a corresponding sequence of decomposition {(L t , S t ), t = 1, 2, 3, ...}. This sequence is obtained via a simple modification of the ALM algorithm [188] to allow warm restarts, that is, where the initial ALM iterate for computing (L t , S t ) is (L t-1 , S t-1 ). Finally, Ramirez and Sapiro [253][254] select the pair (L t, S t), t = arg min t {M DL(L t ) + M DL(S t )} where M DL(L t ) + M DL(S t ) = M DL(A|M ) denoted the description length in bits of A under the description provided by a given model M ∈ M. Experimental results show that the best λ is not the one determined by the theory in Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF].

Another approach developed by Gao et al. [START_REF] Gao | Block-sparse RPCA for consistent foreground detection[END_REF] consist of two-pass RPCA process. The first-pass RPCA done on block resolution detect region with salient motion. Then, a saliency measure in each area is computed and permits to adapt the value of λ following the motion in the second-pass RPCA. Experimental results show that this block-sparse RPCA outperforms the original PCP [START_REF] Candes | Robust principal component analysis? International[END_REF] and the ReProCS [243]. In a similar way using a block-based RPCA, Biao and Lin [START_REF] Biao | Robust foreground detection using block based RPCA[END_REF] determined λ with the affiliation of block to the class of "moving objects". Experimental results show that this approach gives better robustness on the I2R dataset than the single Gaussian, the MOG and the KDE.

Modified-PCP

In the literature, there are several modifications which concern the improvements of the original PCP and they can be classified as follows:

1. Fixed rank: Leow et al. [174] proposed a fixed rank algorithm for solving background recovering problems because low-rank matrices have known ranks.

The decomposition involves the same model than PCP in Equation 20but the minimization problem differs by the constraint as follows:

min

L,S ||L|| * + λ||S|| l1 subj rank(L) = r ( 26
)
with r is the rank of the matrix L and r is known. Lai et al. [170] developed an incremental fixed rank algorithm.

2. Nuclear norm free: Another variant of PCP was formulated by Yuan et al.

[367] who proposed a nuclear-norm-free algorithm to avoid the SVD computation. The low-rank matrix is thus represented as u1 T where u ∈ R m and 1 denotes the vector R n . Accordingly, a noiseless decomposition is formulated as follows:

A = u1 T + S (27)
Then, the corresponding minimization problem is the following one:

min u ||A -u1 T || l1 (28)
Note the closed-form solution of Equation 28 is given by the median of the entries of the matrix A. In other words, the background is extracted as the median at each pixel location of all frames of a surveillance video. As no iteration is required at all to obtain the solution of 28, its computation for solving should be significantly cheaper than any iterative schemes for solving Equation 22numerically. Furthermore, this model extracts more accurately the background than the original PCP. Moreover, Yuan et al. [367] developed a noise and a blur and noise nuclear-norm-free models SPCP which are detailed in Section 2.2.4. In a similar way, Yang et al. [354] proposed a nonconvex model for background/foreground separation, that can incorporate both the nuclearnorm-free model and the use of nonconvex regularizers.

3. Capped norms: In another way, Sun et al. [299] presented a nonconvex formulation using the capped norms for matrices and vectors, which are the surrogates of the rank function and the l 0 -norm, and called capped nuclear norm and the capped l 1 -norm, respectively. The minimization problem is formulated as follows:

min

L,S rank(L) + λ||S|| l0 subj ||A -L -S|| 2 F ≤ σ 2 (29)
where σ 2 is the level of Gaussian noise. The capped nuclear norm is then:

1 θ 1 ||L|| * + p i=1 max(σ i (L) -θ 1 ), 0) (30)
and the capped l 1 -norm is formulated as follows:

1

θ 2 ||S|| l1 + p i=1 max(S ij ) -θ 2 ), 0) (31)
for some small parameters θ 1 and θ 2 . If all the singular values of L are greater than θ 1 and all the absolute values of elements in S are greater than θ 2 , then the approximation will become equality. The smaller θ 1 and θ 2 are, the more accurate the capped norm approximation would be. The recovery precision is controled via θ 1 and θ 2 . By carefully choosing θ 1 and θ 2 , L and S are more accurately determined than with the nuclear norm and the l 1 -norm approximation. This capped formulation can be solved via two algorithms. One is based on the Difference of Convex functions (DC) framework and the other tries to solve the sub-problems via a greedy approach. Experimental results [299] show better performance for the capped formulation of PCP than the original PCP [START_REF] Candes | Robust principal component analysis? International[END_REF] and SPCP [393] on the I2R dataset [177].

4. Inductive approach: Bao et al. [18] proposed the following decomposition:

A = P A + S (32)
where P ∈ R n×n is the low-rank projection matrix. The related optimization problem is formulated as follows:

min

P,S ||P || * + λ||S|| l1 subj A -P A -S = 0 ( 33
)
This is solved by IALM [188]. Furthermore, Bao et al. [18] developed an inductive version which requires less computational cost in processing new samples.

5. Partial Subspace Knowledge: Zhan and Vaswani [369] proposed a modified-PCP with partial subspace knowledge. They supposed that a partial estimate of the column subspace of the low-rank matrix L is available. This information is used to improve the PCP solution, i.e. allow recovery under weaker assumptions. So, the modified-PCP requires significantly weaker incoherence assumptions than PCP, when the available subspace knowledge is accurate.

The corresponding optimization problem is written as follows:

min

L,S ||L|| * + λ||S|| l1 subj L + P Γ ⊥ S = P Γ ⊥ A (34)
where P Γ ⊥ is a projection matrix, Γ is a linear space of matrices with column span equal to that of the columns of S, and Γ ⊥ is the orthogonal complement. Zhan and Vaswani [369] applied with success their modified-PCP to the background-foreground separation problem, in which the subspace spanned by the background images is not fixed but changes over time and the changes are gradual.

6. Schatten-p,lq-PCP (p,q-PCP): The introduced norms by Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] are not tight approximations, which may deviate the solution from the authentic one. Thus, Wang et al. [320] considered a non-convex relaxation which consists of a Schatten-p norm and a lq-norm with 0 < p, q ≤ 1 that strengthen the low-rank and sparsity, respectively. The Schatten-p norm (||.|| Sp ) is a popular non-convex surrogate of the rank function. Thus, the miminization problem is the following one: min

L,S ||L|| p Sp + λ||S|| lq subj A -L -S = 0 (35)
By replacing the Schatten-p norm and a lq-norm by their expression, the miminization problem can be written as follows:

min

L,S λ 1 min(m,n) i=1 (σ i (L)) p + λ 2 m i=1 n j=1 |S ij | q (36)
where σ i (L) denotes the i th singular values of L. When p = q = 1, p,q-PCP degenerates into the original convex PCP. Smaller values of p and q help p,q-PCP to well approximate the original formulation of RPCA. The solver used is a Proximal Iteratively Reweighted Algorithm (PIRA) based on alternating direction method of multipliers, where in each iteration the underlying objective function is linearized to have a closed form solution. Experimental results [320] on the I2R dataset [177] show better performance for p,q-PCP (in its stable formulation) in comparison to the original SPCP [393] and SPCP solved via NSA [12].

7. Modified Schatten-p,lq-PCP: Shao et al. [276] proposed a similar approach than p,q-PCP [320] but they used the Lq-seminorm as a surrogate to the l 1norm instead of the lq-norm. Thus, the miminization problem is the following one: min

L,S ||L|| p Sp + λ||S|| q Lq subj A -L -S = 0 (37)
Furthermore, Shao et al. [276] used two different solvers based on the ALM and the APG methods as well as efficient root-finder strategies.

Robust 2D-PCA:

To take into account the two-dimensional spatial information, Sun et al. [301] extracted a distinguished feature matrix for image representation, instead of matrix-to-vector conversion. Thus, the miminization problem is the following one:

min U,V,S λ||S|| l0 + 1 2 ||A -U ΣV T -S|| 2 F subj U U T = I, V V T = I ( 38
)
where

U ΣV T = L Different from l 1 -
L,S min(m,n) i=N +1 ||σ i (L)|| + λ||S|| l1 subj A -L -S = 0 (39)
where σ i (L) represents the i th singular value of the low-rank matrix L, and N is a constraint parameter for rank-N. Minimizing partial sum of singular values can minimize the rank of the matrix L and satisfy rank-N constraint. Then, Oh [231] applied the RPCA with Rank-1 Soft Constraint on the edge images for moving objects detection under global illumination changes.In the case of moving camera, Ebadi and Izquierdo [75][77] proposed a SVD-free algorithm to solve Rank-1 RPCA that achieved more than double the amount of speed-up in computation time for the same performance target compared to its counterpart with SVD. This approach [START_REF] Ebadi | Approximated robust principal component analysis for improved general scene background subtraction[END_REF] can handle camera movement, various foreground object sizes, and slow-moving foreground pixels as well as sudden and gradual illumination changes.

10. JVFSD-RPCA: Wen et al. [332] reconstructed the input video data and aimed to make the foreground pixels not only sparse in space but also sparse in "time" by using a Joint Video Frame Set Division and RPCA-based (JVFSD-RPCA) method. In addition, they used the motion as a priori knowledge. The proposed method consists of two phases. In the first phase, wa Lower Boundbased Within-Class Maximum Division (LBWCMD) method divided the video frame set into several subsets. In this way, the successive frames are assigned to different subsets in which the foregrounds are located at the scene randomly. In the second phase, aach subset with the frames are increased with a small quantity of motion. This method show robustness in the case of dynamic backgrounds.

11. NSMP/WNSMP: Wang and Feng [322] improved the RPCA method to find a new model to separate the background and foreground and reflect the correlation between them as well. For this, they proposed a "low-rank + dual" model and they use the reweighted dual function norm instead of the normal norms so as to get a better and faster model. So, the original minimization problem is improved by a nuclear norm and spectral norm minimization problem (NSMP). Thus, the minimization problem with the dual norm is the following one: min

L,S λ||L|| * + µ||S|| 2 subj A -L -S = 0 ( 40
)
where the spectral norm ||.|| 2 is the dual norm of the nuclear norm, and it corresponds to the largest singular value of the matrix [322]. As the nuclear norm regularized is not a perfect approximation of the rank function, Wang and Feng [322] proposed a weighted function nuclear norm and spectral norm minimization problem (WNSMP) with the corresponding minimization problem: min

L,S λ||ω(L)|| * + µ||ω -1 (S)|| 2 subj A -L -S = 0 ( 41
)
where ω() denotes the weighted function which directly adds the weights onto the singular values of the matrix, and, for any matrix X, weighted function norm is defined as follows:

||ω(X)|| * = min(m,n) i=1 ω i σ i (X) and ||ω -1 (X)|| 2 = max i 1 ωi σ i (X)
w where σ i (X) represents the i th singular value of the matrix X.

Although this minimization problem with the weighted function nuclear norm is nonconvex, fortunately it has a closed form solution due to the special choice of the value of weights, and it is also a better approximation to the rank function. NSMP and WNSMP show more robustness on the I2R dataset [177] than RPCA solved IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and GoDec [389].

12. Implicit Regularizers: He et al. [126] proposed a robust framework for lowrank matrix recovery via implicit regularizers of robust M -estimators (Huber, Welsch, l 1 -l 2 ) and their minimizer functions. Based on the additive form of halfquadratic optimization, proximity operators of implicit regularizers are developed such that both low-rank structure and corrupted errors can be alternately recovered. Thus, the minimization problem with implicit regularizers is formulated as follows:

min

L,S λ||L|| * + ϕ(S) subj A -L -S = 0 ( 42
)
where the implicit regularizer ϕ(y) is defined as the dual potential function of a robust loss function φ(x) where φ(x) = min

y 1 2 ||x -y|| 2 2 + ϕ(x). If φ(x)
is Huber M -estimator, the implicit regularizer ϕ(y) becomes µλ||(L)|| l1 . When the M -estimator φ(x) is Welsch M -estimator, the minimization problem becomes the sample based maximum correntropy problem. Compared with the mean square error, the model in Equation ( 40) is more robust to outliers due to Mestimation. Experimental results [126] on the I2R dataset [177] show that the Welsch M -estimator outperforms the Huber-estimator and the l 1 -l 2 -estimator.

RPCA via Stable Principal Component Pursuit

PCP is limited to the low-rank component being exactly low-rank and the sparse component being exactly sparse but the observations in real applications are often corrupted by noise affecting every entry of the data matrix. Therefore, Zhou et al.

[393] proposed a stable PCP (SPCP) that guarantee stable and accurate recovery in the presence of entry-wise noise. In the following sub-sections, we reviewed this method and all these modifications in terms of decomposition, solvers, incremental algorithms and real time implementations. Table 16 shows an overview of the Stable Principal Component Pursuit methods and their key characteristics.

Stable Principal Component Pursuit

Zhou et al. [393] proposed a stable PCP (SPCP) which assumes that the observation matrix A is represented as follows:

A = L + S + E ( 43
)
where E is a noise term (say i.i.d. noise on each entry of the matrix) and ||E|| F < δ for some δ > 0. To recover L and S, Zhou et al. [393] proposed to solve the following optimization problem, as a relaxed version to PCP: min

L,S ||L|| * + λ||S|| l1 subj ||A -L -S|| F < δ ( 44
)
where ||.|| F is the Frobenius norm and λ = 1 √ n .

Algorithms for solving SPCP

Just as in Equation ( 21) for PCP, Tao and Yuan [307] showed that an easy reformulation of the constrained convex programming for Equation ( 44) falls perfectly in the applicable scope of the classical ALM. Moreover, the favorable separable structure emerging in both the objective function and the constraints entails the idea of splitting the corresponding augmented Lagrangian function to derive efficient numerical algorithms. So, Tao and Yuan [307] developed the alternating splitting augmented Lagrangian method (ASALM) and its variant (VASALM), and the parallel splitting augmented Lagrangian method (PSALM) for solving Equation [START_REF] Cevher | Sparse signal recovery using markov random fields[END_REF]. ASALM and its variants converge to an optimal solution. However, ASALM iterations are too slow for real time application and its complexity is not known.

To address this problem, Aybat et al. [12] studied how fast first-order methods can be applied to SPCP with low complexity iterations and showed that the subproblems that arise when applying optimal gradient methods of Nesterov, alternating linearization methods and alternating direction augmented Lagrangian methods to the SPCP problem either have closed-form solutions or have solutions that can be obtained with very modest effort. Furthermore, Aybat et al. [12] developed a new first order algorithm called Non Smooth augmented Lagrangian Algorithm (NSA), based on partial variable splitting. All but one of the methods analyzed require at least one of the non-smooth terms in the objective function to be smoothed and obtain an ǫ-optimal solution to the SPCP problem in O( √ ǫ) iterations. NSA, which works directly with the fully non-smooth objective function, is proved to be convergent under mild conditions on the sequence of parameters it uses. NSA, although its complexity is not known, is the fastest among the optimal gradient methods, alternating linearization methods and alternating direction augmented Lagrangian methods algorithms and substantially outperforms ASALM. In a similar way, Aybat et al. [START_REF] Aybat | Efficient algorithms for robust and stable principal component pursuit[END_REF] proposed a proximal gradient algorithm called Partially Smooth Proximal Gradient (PSPG). Experimental results show that both the number of partial SVDs and the CPU time of PSPG are significantly less than those for NSA and ASALM. An overview of these algorithms as well as their complexity can be seen in Table 11. et al. [209] proposed a real time implementation framework, entitled Divide-Factor-Combine (DFC). DFC randomly divides the original matrix factorization task into cheaper subproblems, solves those subproblems in parallel using any base matrix factorization (MF) algorithm, and combines the solutions to the subproblem using an efficient technique from randomized matrix approximation. The inherent parallelism of DFC allows for near-linear to superlinear speedups in practice, while the theory provides high-probability recovery guarantees for DFC comparable to those provided by its base algorithm. So, Mackey et al. [209] proposed two algorithms called DFC-PROJ and DFC-NYS, that differ from the method used to divide the original matrix. DFC-PROJ randomly partitions the orthogonal projection of the matrix A onto the t l-column submatrices C 1 , ..., C t by using a column sampling method, while DFC-NYS selects an l-column submatrix C and an d-row submatrix R using the generalized Nyström method. DFC significantly reduces the per-iteration complexity to O(mlr C1 ) where r C1 is the rank of the matrix C 1 for the DFC-PROJ. The cost of combining the submatrix estimates is even smaller, since the outputs of standard MF algorithms are returned in factored form. Indeed, the column projection step of DFC-PROJ requires only [209] improved these real time implementations by using ensemble methods that improve performance of matrix approximation algorithms, while straightforwardly leveraging the parallelism of modern many-core and distributed architectures [165].

Methods for real time implementation of SPCP

Mackey

As such, an ensemble variants of the DFC algorithms have been developed reducing recovery error while introducing a negligible cost to the parallel running time.

For DFC-PROJ-ENS, rather than projecting only onto the column space of C 1 , the projection of C 1 , ..., C t is done onto the column space of each C i in parallel and then average the t resulting low-rank approximations.For DFC-NYS-ENS, a random d-row submatrix is chosen as in DFC-NYS and independently partition the columns of the matrix in l as in DFC-PROJ. After running the base MF algorithm on each submatrix, the generalized Nyström method is applied to each pair of matrices in parallel and then the t resulting low-rank approximations is obtained by average.

Modified-SPCP

In the literature, There are several modifications which concern the original SPCP and they can be classified as follows:

1. Bilateral factorization: Zhou and Tao [390] proposed a noisy robust PCA by replacing L with its bilateral factorization L = U V and regularizing the l 1 norm of S's entries. The corresponding minimization problem is then formulated as follows: min

U,V,S λ||S|| l1 + ||A -U V -S|| 2 F subj rank(U ) = rank(V ) ≤ r (45)
The l 1 regularization permits soft-thresholding in updating S. Zhou and Tao [390] solved this minimization problem using a Greedy Bilateral Smoothing algorithm (GreBsmo). GreBsmo considerably speed up the decomposition and performs 30-100 times faster than most existing algorithms such as IALM [188].

2. Nuclear norm free: Other variants of SPCP were given by Yuan et al. [367] who developed a nuclear-norm-free algorithm to avoid the SVD computation in SPCP. The low-rank matrix is thus represented as u1 T where u ∈ R m and 1 denotes the vector R n . Accordingly, a noiseless decomposition is formulated as follows:

A = u1 T + S + E (46)
Then, the corresponding minimization problem is the following one:

min S∈R m×n ,u∈R m ||S|| l1 + µ 2 ||A -u1 T -S|| 2 F subj rank(u1 T ) = 1 ( 47
)
where µ is a penalty parameter. This model has no closed-form solution and need to be solved iteratively. This model extracts more accurately the background than the original SPCP. Another advantage of this model against the original SPCP is that it has only one parameter in the minimization function to be tuned in the implementation. Considering, that there might be a blur in a noisy video surveillance video, Yuan et al. [367] extended the model developed in Equation 47to:

min S∈R m×n ,u∈R m min ||S|| l1 + µ 2 ||A + H(U 1 T + S)|| 2 F subj rank(u1 T ) = 1 ( 48
)
where H is the matrix representation of a regular blurring operator. The blur is assumed to appear in a frame-wise way. In this blur configuration, this nuclear-norm-free model extracts more robustly the background than the original SPCP with the blur option.

3. Under-complete dictionary: Considering a dictionary approach, Sprechman et al. [295] proposed the following decomposition:

A = U T + S + E ( 49
)
where U ∈ R m×r , T ∈ R r×n , and S ∈ R m×n . Then, the corresponding minimization problem is as follows:

min

U,T,S λ||S|| l1 + λ 1 2 (||U|| 2 F + ||T || 2 F ) + 1 2 ||A -U T -S)|| 2 F subj rank(U T ) ≤ r (50)
The low-rank component can be considered as an under-complete dictionary U , with r atoms, multiplied by a matrix T containing in its columns the corresponding coefficients for each data vector in A. This interpretation brings the SPCP problem close to that of dictionary learning in the sparse modeling domain. This is solved via an alternating minimization problem. Furthermore, Sprechman et al. [295] developed an online version of their SPCP via fast trainable encoders.

4.

Variational formulation: Aravkin et al. [START_REF] Aravkin | A variational approach to stable principal component pursuit[END_REF] proposed a convex variational framework which is accelerated with quasi-Newton methods. The corresponding minimization problem is then formulated as follows:

min Φ(L, S) subj ρ(L + S -A) ≤ ǫ
where in classical formulation Φ(L, S) = ||L|| * + λ||S|| l1 and ρ is assumed to be the Frobenius norm. As this restriction is not necessary, Aravkin et al. [START_REF] Aravkin | A variational approach to stable principal component pursuit[END_REF] considered ρ to be smooth and convex. Then, ρ is taken to be the robust Huber penalty. This approach offers advantages over the original SPCP formulation in terms of scalability and practical parameter selection.

5. Three Term Low-rank Optimization: Oreifej et al. [236] proposed a three term decomposition for video stabilization and moving object detection in turbulence as follows:

min

A,L,S λ||L|| * + λ 1 ||fπ(S)|| + λ 2 ||E|| 2 F subj A = L + S + E (51)
where the frames of the sequence are stacked in the matrix A. Thus, low-rank matrix L corresponds to the background, the sparse matrix E corresponds to the moving object and the dense error matrix E corresponds to the turbulence. The turbulence causes dense and Gaussian noise, and therefore can be captured by Frobenius norm. Therefore, Oreifej et al. [236] enforced an additional constraint on the objects with fπ(.) which is the object confidence map, which is a linear operator that weights the entries of S according to their confidence of corresponding to a moving object such that the most probable elements are unchanged and the least are set to zero. Oreifej et al. [236] used the solver IALM [188] to solve Equation 51. Experimental results [236] show that this decomposition outperforms Mixture of Gaussians (MoG) [298], Kernel Density Estimation [START_REF] Elgammal | Non-parametric model for background subtraction[END_REF] and PCA [233]. Furthermore, the code called ThreeWayDec32 is provided.

RPCA via Quantization based Principal Component Pursuit

Becker et al. [20] proposed a inequality constrained version of RPCA proposed by Candes et al. [START_REF] Candes | Robust principal component analysis? International[END_REF] to take into account the quantization error of the pixel's value. Indeed, each pixel has a value between 0, 1, 2, . . . , 255. This value is the quantized version of the real value which is between [0, 255]. So, the idea is to apply RPCA to the real observations instead of applying it to the quantized ones. Indeed, it is unlikely that the quantized observation can be split nicely into a low-rank and sparse component. So, Becker et al. [20] supposed that L + S is not exactly equal to A, but rather that L + S agrees with A up to the precision of the quantization. The quantization can induce at most an error of 0.5 in the pixel value. This measurement model assumes that the observation matrix A is represented as follows:

A = L + S + Q ( 52
)
where Q is the error of the quantization. Then, the objective function is the same than the equality version in Equation ( 22), but instead of the constraints L+S = A, the constraints are A -L -S l∞ ≤ 0.5. So, the quantization based PCP (QPCP) is formulated as follows:

min

L,S ||L|| * + λ||S|| l1 subj A -L -S l∞ ≤ 0.5 (53)
The l∞-norm allows to capture the quantization error of the observed value of the pixel.

Algorithms for solving QPCP: Becker et al. [20] used a general framework for solving this convex cone problem called Templates for First-Order Conic Solvers (TFOCS). First, this approach determines a conic formulation of the problem and then its dual. Then, Becker et al. [20] applied smoothing and solved the problem using an optimal first-order method. This approach allows to solve the problem in compressed sensing.

L,S ||L|| * + κ(1 -λ)||L|| l2,1 + κλ||S|| l2,1 subj A -L -S = 0 ( 54
)
where ||.|| * and ||.|| l2,1 are the nuclear norm and the l 2,1 -norm, respectively. The l 2,1 -norm corresponds to the l 1 -norm of the vector formed by taking the l 2 -norms of the columns of the underlying matrix. The term κ(1λ)||L|| l2,1 ensures that the recovered matrix L has exact zero columns corresponding to the outliers. In order to eliminate ambiguity, the columns of the low-rank matrix L corresponding to the outlier columns are assumed to be zeros.

Algorithm for solving BPCP: Tang and Nehorai [326] designed an efficient algorithm to solve the convex problem in Equation ( 54) based on the ALM method. This algorithm decomposed the matrix A in a low-rank and block-sparse matrices in respect to the ||.|| l2,1 and the extra term κ(1λ)||L|| l2,1 .

RPCA via Local Principal Component Pursuit

PCP is highly effective but the underlying model is not appropriate when the data are not modeled well by a single low-dimensional subspace. Wohlberg et al.

[333] proposed a decomposition corresponding to a more general underlying model consisting of a union of low-dimensional subspaces.

A = AU + S (55)
This idea can be implemented as the following problem:

min U,S α||U || l1 + β||U || l2,1 + β||S|| l1 subj A -AU -S = 0 (56)
The explicit notion of low-rank, and its nuclear norm proxy, is replaced by representability of a matrix as a sparse representation on itself. The l 2,1 -norm encourages rows of U to be zero, but does not discourage nonzero values among the entries of a nonzero row. The l 1 -norm encourages zero values within each nonzero row of S.

To better handle noisy data, Wohlberg et al. [333] modified the Equation ( 56) with a penalized form and add a Total Variation penalty on the sparse deviations for contiguous regions as follows:

min

U,S 1 2 ||A -DU -S|| 2 l2 + α||U || l1 + β||U || l2,1 + β||S|| l1 + δ||grad(S)|| l 1 subj A -DU -S = 0 ( 57
)
where the dictionary D is derived from the data A by mean subtraction and scaling, and grad(S) is a vector valued discretization of the 3D gradient of S. An appropriate sparse U can be viewed as generating a locally low-dimensional approximation DU of A -S. When the dictionary is simply the data (i.e., D = A), the sparse deviations (or outliers) S are also the deviations of the dictionary D, so constructing the locally low-dimensional approximation as (D -S)U , implying an adaptive dictionary D -S, should allow U to be even sparser.

Algorithm for solving LPCP: Wohlberg et al. [333] proposed to solve the Equation ( 56) using the Split Bregman Algorithm (SBA) [START_REF] Goldstein | The split Bregman method for l 1 -regularized problems[END_REF]. Adding terms relaxing the equality constraints of each quantity and its auxiliary variable, Wohlberg et al. [333] introduced Bregman variables in the Equation ( 56). So, the problem is split into an alternating minimization of five subproblems. Two subproblems are l 2 problems that are solved by techniques for solving linear systems such as conjugate gradient. The other three subproblems are solved very cheaply using shrinkage, i.e generalized shrinkage and soft shrinkage.

2.6 RPCA via Outlier Pursuit Xu et al. [342] proposed a robust PCA via Outlier Pursuit (OP) to obtain a robust decomposition when the outliers corrupted entire columns, that is every entry is corrupted in some columns. This method involves the nuclear norm minimization and recover the correct column space of the uncorrupted matrix, rather than the exact matrix itself. The decomposition involves the same model than PCP in Equation (20), that is A = L + S. A straightforward formulation to minimize the energy function can be written as follows:

min

L,S rank(L) + λ||S|| 0,c subj A -L -S = 0 ( 58
)
where ||S|| 0,c stands for the number of nonzero columns of a matrix, and it is equivalent to ||S|| l2,0 which corresponds to the number of non-zero columns too [375]. λ > 0 is an arbitrary balanced parameter. But this problem is N P -hard, typical solution might involve a search with combinatorial complexity. This research seeks to solve for L with the following optimization problem: min

L,S ||L|| * + λ||S|| l1,2 subj A -L -S = 0 ( 59
)
where ||.|| * and ||.|| l1,2 are the nuclear norm and the l 1,2 -norm, respectively. The l 1,2 -norm corresponds to the l 2 -norm of the vector formed by taking the l 1 -norms of the columns of the underlying matrix. λ > 0 is an arbitrary balanced parameter.

Adapting the OP algorithm to the noisy case,that is A = L + S + E, Xu et al. [342] proposed a robust PCA via Stable Outlier Pursuit (SOP): min

L,S ||L|| * + λ||S|| l1,2 subj ||A -L -S|| F < δ ( 60
)
where S is supported on at most γn columns and λ = 3 7 √ γn .

Algorithm for solving OP and SOP: Xu et al. [342] used the Singular Value Threshold (SVT) algorithm to solve these two minimization problems.

RPCA with Sparsity Control

Mateos and Giannakis [213][214] proposed a robust PCA using a bilinear decomposition with Sparsity Control (RPCA-SpaCtrl). The decomposition involves the following model:

A = M + P U T + S + E ( 61
)
where M is the mean matrix, the matrix U has orthogonal columns, P are the principal components matrix, S is the outliers matrix and E is a zero-mean matrix.

The percentage of outliers determines the degree of sparsity in S. The criterion for controlling outlier sparsity is seek to the relaxed estimation: min

U,S ||X + 1 N M T -P U T -S|| 2 F + λ||S|| l 2(r) subj U U T = Iq (62)
where ||S|| l 2(r) = m×n i=1 ||S i || l 2 is the row-wise l 2 -norm. The non-differentiable l 2norm regularization term controls rows-wise sparsity on the estimator of S. The sparsity is then also controlled by the parameter λ. To optimize Equation ( 62), Mateos and Giannakis [213][214] used an alternating minimization algorithm [396].

Algorithm for incremental SpaCtrl: An incremental version of Equation ( 62) is obtained using the Exponentially Weighted Least Squares (EWLS) estimator as follows:

min U,S m×n i=1 β (m×n)-i ||Xn + m -U T T i -S i || 2 l2 + λ||S i || l2 ( 63
)
where β is a learning rate between 0 and 1. So, the entire history of data is incorporated in the online estimation process. Whenever β < 1, past data are exponentially discarded thus enabling operation in nonstationary backgrounds. Towards deriving a real-time, computationally efficient, and recursive solver of Equation ( 63), an AM scheme is adopted in which iteration k coincides with the time scale i = 1, 2, ... of data acquisition. Experimental results [214] show that RPCA-SpaCtrl with λ = 9.69 × 10 -4 presents better performance than the naive PCA [233] and RSL [310] with less time computation.

RPCA with Sparse Corruption

Even if the matrix A is exactly the sum of a sparse matrix S and a low-rank matrix L, it may be impossible to identify these components from the sum. For example, the sparse matrix S may be low-rank, or the low-rank matrix L may be sparse. So, Hsu et al. [130] imposed conditions on the sparse and low-rank components in order to guarantee their identifiability from A . This method requires that S not be too dense in any single row or column, and that the singular vectors of L not be too sparse. The level of denseness and sparseness are considered jointly in the conditions in order to obtain the weakest possible conditions. This decomposition RPCA with Sparse Corruption (RPCA-SpaCorr) involves the same model than PCP in Equation ( 20), that is A = L + S. Then, Hsu et al. [130] proposed two convex formulations. The first is the following constrained formulation: min

L,S ||L|| * + λ||S|| l1 subj ||A -L -S|| l1 ≤ ǫ 1 and ||A -L -S|| * ≤ ǫ * (64)
where λ > 0, ǫ 1 ≥ 0 and ǫ * ≥ 0. The second is the regularized formulation: min

L,S ||L|| * + λ||S|| l1 + 1 2µ ||A -L -S|| 2 F subj ||A -L -S|| l1 ≤ ǫ 1 and ||A -L -S|| * ≤ ǫ * ( 65
)
where µ > 0 is the regularization parameter. Hsu et al. [130] added a constraint to control the entry-wise 64) and ||A -S|| l∞ ≤ b is added in Equation ( 65). The parameter b is a natural bound for L and is typically 510 for image processing. Hsu et al. [130] determined two identifiability conditions that guarantee the recovery. The first measures the maximum number of non-zero entries in any row or column of S. The second one measures the sparseness of the singular vectors L. Then, a mild strengthening of these measures is achieved for the recovery guarantees.

∞-norm of L, that is ||L|| l∞ . That is ||L|| l∞ ≤ b is added in Equation(

RPCA via Log-sum Heuristic Recovery

When the matrix has high intrinsic rank structure or the corrupted errors become dense, the convex approaches may not achieve good performances. Then, Deng et al. [66] used the Log-sum Heuristic Recovery (LHR) to learn the low-rank structure. The decomposition involves the same model than PCP in Equation ( 20), that is A = L + S. Although the objective in Equation ( 22) involves the nuclear norm and the l 1 -norm, it is based on the l 1 heuristic since nuclear norm can be regarded as a specific case of l 1 -norm [START_REF] Deng | Low-rank structure learning via log-sum heuristic recovery[END_REF]. Replacing the nuclear norm by its l 1 -norm formulation, the problem can be solved as follows: [START_REF] Deng | Low-rank structure learning via log-sum heuristic recovery[END_REF] where X = {Y, Z, L, S} and D = (Y, Z, L, S) :

min X∈ D 1 2 (||diag(Y )|| l1 + ||diag(Z)|| l1 ) + λ||E|| l1 (
Y L L T Z ≥ 0, (L, S) ∈ C
(L, S) ∈ C stands for convex constraint. Y and Z are both symmetric and positive definite. ≥ represents semi-positive definite. The convex problem with two norms in Equation ( 22) has been converted to an optimization only with l 1 -norm and therefore it is called l 1 -heuristic. Next, Deng et al. [66] used the logsum term to represent the sparsity of signals and obtained the Log-sum Heuristic Recovery (LHR) model:

min X∈ D 1 2 (||diag(Y )|| L + ||diag(Z)|| L ) + λ||E|| L ( 67
)
where

||X|| L = ij log (|X i j| + δ) with δ > 0 is a small regularization constant.
This model is non convex but the convex upper bound can be easily defined. LHR can remove much denser errors from the corrupted matrix rather compared to PCP.

Algorithm for solving LHR: Deng et al. [START_REF] Deng | Low-rank structure learning via log-sum heuristic recovery[END_REF] used the majorization-minimization (MM) [START_REF] Fazel | Matrix rank minimization with applications[END_REF][171] algorithm that replaces the hard problem by a sequence of easier ones. It proceeds in an Expectation Maximization (EM)-like fashion by repeating two steps of majorization and minimization in an iterative way. During the majorization step, it constructs the convex upper bound of the non-convex objective.

In the minimization step, it minimizes the upper bound.

2.10 RPCA via Iteratively Reweighted Least Squares Minimization Guyon et al. [113] proposed the decomposition solved via IRLS with the following model:

A = L + S = U V + S ( 68
)
where U is a low-rank matrix corresponding to the background model plus noise, and V reconstructs L by linear combination. S corresponds to the moving objects. The model involves the error reconstruction determined by the following constraints:

min U∈R n×p ,V ∈R p×m µ||U V || * + ||(A -U V ) • W 1 || l α,β (69)
where ||.|| * denote the nuclear norm and ||.|| l α,β is a mixed norm. W 1 which is a weighted matrix is iteratively computed and aims to enforce the fit exclusively on guessed background region. A function Φ(.) smoothes the error like spatial median filtering and transforms the error for obtain a suitable weighted mask for regression:

W = Φ(A -U V ), Φ(x) = e -γT V (A-U V ) (70)
By including local penalty as a constraint in RPCA, it explicitly increases local coherence of the sparse component as filled/plain shapes for moving objects. Furthermore, the decomposition is split into two parts. The first part tracks 1-Rank decomposition since the first eigen-vector is strongly dominant in video surveillance. For the mixed norm, Guyon et al. [113] used ||.|| l2,1 instead of the usual ||.|| l1,1 because it forces spatial homogeneous fitting. Thus, the SVD algorithm can be seen as an iterative regression and then IRLS algorithm is then used. So, Guyon et al. [113] increased local coherence of the error for moving objects by including local penalty as a constraint in the decomposition. Using the same approach, Guyon et al. [112] added spatial constraint in the minimization based on the gradient and Guyon et al. [111] proposed a spatio-temporal version. Another variant of RPCA via IRLS have been developed by Lu et al [202]. 2.11 RPCA via Stochastic Optimization Feng et al. [83] proposed an Online Robust PCA (OR-PCA) algorithm. The main idea is to develop a stochastic optimization algorithm to minimize the empirical cost function, which processes one sample per time instance in an online manner.

The coefficients which correspond to noise and the basis are optimized in an alternative manner. The low dimensional subspace called low-rank matric basis are first initialized randomly and then updated after every frame per time instance. Moreover, OR-PCA decomposes the nuclear norm of the objective function of the traditional PCP algorithms into an explicit product of two low-rank matrices, i.e. basis and coefficients. The main function in OR-PCA is formulated as: min

L∈R n×p ,R∈R n×r 1 2 ||A -LR T -S|| 2 F + λ 1 2 (||L|| 2 F + ||R|| 2 F ) + λ 2 ||S|| l1 ( 71
)
where R is a coefficient matrix. λ 1 controls the basis and coefficients for low-rank matrix, whereas λ 2 controls the sparsity pattern, which can be tunned according to video analysis. In addition, basis and coefficient depend on the value of rank.

In case of video background modeling, no visual results [START_REF] Feng | Online robust PCA via stochastic optimization[END_REF] have been found using this technique. Therefore, Javed et. al [144] modified OR-PCA via stochastic optimization method for background subtraction applications. An initialization scheme is adopted which converges the algorithm very fastly as compared to original OR-PCA.

In order to perform OR-PCA, a number of video frames are first initialized as a low dimensional basis then stochastic optimization is performed on each input frame to separate the low-rank and sparse component. As compare to conventional RPCA via PCP based schemes, no batch optimizations are needed therefore OR-PCA is applicable for real time processing. In addition, a global pre-processing steps such as Laplacian and Gaussian images are introduced in modified OR-PCA which increase the detection rate. Using these modifications in original scheme, both memory cost and computational time is decreased, since the idea is based on to process one single frame per time instance, but the method shows some weak performance when large variations in the background scenes occurs such as waving trees and water surface. Therefore, Javed et al. [147][143] further improved the foreground segmentation using the continuous constraints such as Markov Random Field (MRF). OR-PCA via image decomposition using initialization scheme including continuous MRF with tuned parameters shows a drastic improvements in the experimental results specially in case highly dynamic backgrounds. In their work a good parameters range is provided according to different background scenarios. A huge amount of experimental results are provided which shows a very nice potential for its

Bayesian Robust Principal Component Analysis

Bayesian Robust Principal Component Analysis approaches have also been investigated for RPCA and used a Bayesian framework in the decomposition into low-rank plus sparse matrices. Ding et al. [71] modeled the singular values of L and theentries of S with beta-Bernoulli priors, and used a Markov chain Monte Carlo (MCMC) sampling scheme to perform inference. This method called Bayesian RPCA (BRPCA) needs many sampling iterations, always hampering its practical use. In a similar approach, Babacan et al. [17] adopted the automatic relevance determination (ARD) approach to model both L and S, and utilized the variational Bayes (VB) method to do inference. This method called Variational Bayesian RPCA (VBRPCA) is more computationally efficient. However, these three methods assume a certain noise prior (a sparse noise plus a dense noise), which cannot always effectively model the diverse types of noises occurring in practice. To address this problem, Zhao et al. [385] proposed a generative RPCA model under the Bayesian framework by modeling data noise as a mixture of Gaussians (MoG).

Table 17 shows an overview of the Bayesian Robust Principal Component Analysis methods.

2.12.1 Bayesian Robust PCA Ding et al.[71] proposed a Bayesian Robust PCA (BRPCA). Assuming that the observed data matrix A can be decomposed in three matrix as in SPCP [393], the Bayesian model is then as follows:

A = D(ZG)W 2 + BX + E ( 72
)
where D ∈ R n×r , W ∈ R r×m and G ∈ R r×r are diagonal matrices and X ∈ R n×m . The diagonal matrix Z has binaries entries along the diagonal, and the binary matrix B ∈ 0, 1 n×m is sparse. r defines the largest possible rank that may inferred for L, and r is set to a large value. The low-rank, sparse and noise component are obtained as follows.

-Low-rank Component: The low-rank component is modeled as L = D(ZG)W 2 . This is similar to SVD excepted for the extra diagonal matrix Z with diagonal elements z k,k ∈ 0, 1 for k = 1, ..., r. The product ZG is a diagonal matrix too. The use of Z decouples the rank learning and the singular value learning. r is chosen large and then the diagonal entries of ZG are sparse. The binary diagonal matrix Z is modeled as follows: [START_REF] Ebadi | Approximated RPCA for fast and efficient recovery of corrupted and linearly correlated images and video frames[END_REF] with α 0 > 0 and β 0 > 0. The parameters α 0 and β 0 imposed the sparseness of the diagonal of Z. α 0 and β 0 are set respectively to 1/K and (K -1)/K. Each diagonal entry in G, denoted as g k,k for k = 1, ..., r, is obtained from a normal-gamma distribution:

z k,k ∼ Bernoulli(p k ) ( 73
)
p k ∼ Beta(α 0 , β 0 ), k = 1, ...r (
g k,k ∼ N (0, τ -1 )k = 1, ...r. ( 75
) τ ∼ Gamma(a 0 , b 0 ) ( 76
)
with a 0 > 0 and b 0 > 0. a 0 and b 0 are set to 10 -7 . The column of matrices D and W 2 are obtained from normal distribution:

d k ∼ N (0, (1/N)I N )k = 1, ...K. ( 77
)
w 2,m ∼ N (0, (1/K)I K )m = 1, ...M. ( 78
)
with I N is the N × N identity matrix. The decomposition can be rewritten as follows:

lm = D(ZG)w 2,m = K k=1 z k,k g k,k w 2,(k,m) d k , m = 1, ..., M (79)
So, each column of L is the weighted sum of the dictionary elements in D, and K is the size of the dictionary. The weights z k,k g k,k k=1:K determined the dictionary elements that are active to construct L. The weights w 2,(k,m) k=1:K determined the importance of the selected dictionary elements for the representation of the m th column of L. -Sparse Component: The sparse component is modeled as S = BX, where B is a binary matrix. This decomposition separates the learning sparseness from the learning of values. Each column of B is modeled as follows:

bm

∼ N n=1 Bernoulli(πn), m = 1, ..., M (80)
πn ∼ Beta(α 1 , β 1 ), n = 1, ...N. ( 81
)
The sparseness prior is made with the parameters α 1 and β 1 . α 1 and β 1 are set respectively to 1/N and (N -1)/N. The columns of X are obtained from a normal-gamma distribution:

xm ∼ N (0, v -1 I N ), m = 1, ..., M (82)
v ∼ Gamma(c 0 , d 0 ). [START_REF] Feng | Online robust PCA via stochastic optimization[END_REF] with c 0 > 0 and d 0 > 0. c 0 and d 0 are set to 10 -6 . Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] addressed the dependency of the sparse component in time and space with a Markov structure. If the parent node of I t (i, j) noted I t-1 (i, j) is non-zero, its child node is also non-zero with a high probability. To introduce spatial dependence, Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] defined the state of F t (i, j) as follows:

S(F t (i, j)) = active if ||N (F t (i, j)|| 0 ≥ ρ (84)
S(F t (i, j)) = inactive otherwise [START_REF] Gan | Multi-feature robust principal component analysis for video moving object segmentation[END_REF] where ρ = 5 which imposes that a node is active if the sparse component contains at least 5 non-zero members in its neighborhood defined by

N (F t (i, j)) = F (k,l) : |k -i| ≤ 1, |l -j| ≤ 1 .
Then, a child node depends on its parent node in time and on its neighbors in space. Markov dependency is then imposed by modifying the equations ( 80) and ( 81) as follows. [START_REF] Gerhart | Convex optimization techniques and their application in hyperspectral video processing[END_REF] with n = 1, ...N, t = 2, ..., M.

b t ∼ N n=1 Bernoulli(π nt ), t = 1, ..., M (86)
π nt ∼ Beta(α H , β H ) if S(b n,t-1 ) = active with n = 1, ...N, t = 2, ..., M. ( 87
)
π nt ∼ Beta(α L , β L ) if S(b n,t-1 ) = inactive (
where H and L indicate the high and low states in the Markov model and α H ,α H , β L and β H are set to assume that the sparseness will be propagated along time with high probability. For t = 1, the equation ( 80) and ( 81) are used since there are no parent nodes for the first frame. -Noise Component: The noise is modeled by a Gaussian distribution as follows:

en,m ∼ N (0, γ -1 m ), with n = 1, ..., N (89)
γm ∼ Gamma(e 0 , f 0 ) for m = 1, ..., M, [START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex function[END_REF] with en,m is the entry at row n and column m of E. c 0 and d 0 are set to 10 -6 .

Then, the posterior density function of the BRPCA is as follows:

-log (p(Θ|A, H)) = τ 2 ||G|| 2 F -log ([f BB (Z; H))] + N 2 r k=1 ||d k || 2 l 2 + 1 2 M m=1 ||wm|| 2 l 2 + v 2 ||X|| 2 F -log ([f BB (B; H))] + 1 2 ||Y -L -S|| 2 F -log [Gamma(τ |H)Gamma(v|H)Gamma(γ|H)]
+ constant [START_REF] Goldstein | The split Bregman method for l 1 -regularized problems[END_REF] where Θ represents all model parameters, f BB (.|H) represents the beta-Bernoulli prior, and H = {α 0 , α 1 , β 0 , β 1 , a 0 , b 0 , c 0 , d 0 , e 0 , f 0 } are model hyper parameters.

Algorithms for solving BRPCA: Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] proposed to approximate the posterior density function in Equation ( 91) with two algorithms:

-Markov chain Monte Carlo (MCMC) analysis implemented with Gibbs sampler [256]: The posterior distribution is approximated by a set of samples, collected by iteratively drawing each random variable from its conditional posterior distribution given the most recent values of all the other parameters. -Variational Bayesian inference (VB) [19]: A set of distributions q(Θ) allow to approximate the true posterior distributions p(Θ|A), and uses a lower bound to approximate the true log-likelihood of the model log (p(A|Θ). The algorithm iteratively updates q(Θ) so that the lower bound approaches to log (p(A|Θ).

The computational complexity of MCMC and VB iteration is approximatively the same. The VB solution may find a local-optimal solution which may be not be the global-optimal best solution. Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] found that MCMC work quite effectively in practice.

Relation to PCP and SPCP: For the low-rank component, Ding et al. [START_REF] Ding | Bayesian robust principal component analysis[END_REF] employed a Gaussian prior to obtain a constraint on Frobenius norm ||G|| 2 F with a beta-Bernoulli distribution to address the sparseness of singular value and to obtain a small number of non-zero singular values, while PCP employs the rank function that is relaxed to the nuclear norm when solving the problem in a convex way. For the sparse component, the constraint on Frobenius matrix norm ||X|| 2 F and the beta-Bernoulli distribution are used to impose sparseness while PCP use the l 0 -norm that is relaxed to the l 1 -norm. The error term (2µ SPCP [335] corresponds to the Gaussian prior placed on the measurement noise in Equation [START_REF] Goes | Robust stochastic principal component analysis[END_REF]. For solving the problem, the main difference is that BRPCA uses numerical methods to estimate a distribution for the unknown parameters, whereas optimization based methods effectively search a single solution that minimizes a analogous functional to -log (p(Θ|A, H)).

) -1 ||A -L -S|| 2 F in

Variational Bayesian Robust Principal Component Analysis

Babacan et al. [17] proposed a Variational Bayesian Robust PCA (VBRPCA). Assuming that the observed data matrix A can be decomposed in three matrix as in SPCP [393], the variational Bayesian model is then as follows:

A = DB T + S + E ( 92
)
where DB T is the low-rank component with D ∈ R m×r and B ∈ R r×n , S is the sparse component with arbitrarily large coefficients and E is the dense error matrix with relatively smaller coefficients. The low-rank, sparse and noise component are obtained as follows. The low-rank component L is then given by DB T . So, L is the sum of outer-products of the columns of D and B, that is, [START_REF] Goyette | changedetection.net: A new change detection benchmark dataset[END_REF] where k ≥ r. d .i and d i. denote the i th column and row of D, respectively. To impose column sparsity in D and B, such that most columns in D and B are set equal to zero, the columns are defined with Gaussians priors as follows:

L = k i=1 d .i b T .i (
p(D|γ) = k i=1 N (d .i |0, σ i I) (94) p(B|γ) = k i=1 N (b .i |0, σ i I) (95)
where σ i is the variance. Most of the variances are very small values during inference to reduce the rank of the estimate. Then, the following conditional distribution for the observations are obtained:

p(A|D, B, S, β) = N (A|DB T + S, γ -1 I) (96) = exp[ β 2 ||A -DB T -S|| 2 F ( 97
)
where β is a uniform hyperprior. The modeling of the sparse component S is done by using independent Gaussian priors on its coefficients S ij as follows:

p(S|α) = m i=1 n j=1 N (S ij |0, α -1 ij ) ( 98
)
where α = α ij , α ij is the precision of the Gaussian on the (i, j) th coefficient and p(α ij )=const ∀i, j. Finally, the joint distribution is expressed as follows:

p(A, D, B, S, γ, α, β) = p(A|D, B, S, β)p(D|γ)p(B|γ)p(S|α)p(γ)p(α)p(β) (99)
where p(γ i ) = 1 γi a+1 exp( -b γi ) and p(β) is a constant assuming that the noise precision have a uniform prior.

Algorithm for solving VBRPCA: The exact full-Bayesian inference using joint distributions in Equation ( 99) is intractable because p(y) can't be computed by marginalizing all variables. Therefore, Babacan et al. [17] used an inference procedure based on mean field variational Bayes. The aim is to compute posterior distribution approximations by minimizing the Kullback-Leibler divergence in an alternating way for each variable. Let z = (D, D, S, γ, α, β), the posterior approximation q(z k ) of each variable z k ∈ z is then determined as follows:

log (q(z k )) = log (p(A, z)) z z k + const (100)
where z z k is the set z without z k . The distribution p(A, z) is the joint probability distribution given in Equation ( 99). The posterior factorization q(z) = q(z k ) is used such that the posterior distribution of each unknown is estimated by holding the others fixed using their most recent distributions. Thus, the expectations of all parameters in the joint distribution are taken with respect to their most recent distributions, and the result is normalized to find the approximate posterior distribution. Since all distributions are in the conjugate exponential family, the form of each posterior approximation is easily determined.

Factorized Variational Bayesian RPCA (FVBRPCA)

Aicher [5] proposed a Factorized Variational Bayesian RPCA. This model is slightly different from BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF] and VBRPCA [17] in how sparse noise is modeled and incorporated as well as the use of variational Bayes instead of MCMC.

A = U V T + Z * • B + E ( 101
)
where • denotes the Hadamard element-wise multiplication. The low-rank matrix is L = U V T and U is restrited to be an n × r matrix and V to be an r × m matrix so that the rank of L less than or equal to r. The sparse matrix is S = Z * • B and B is set to be a sparse binary matrix and Z * is without constraint. For numerical reasons, Z * is treated as a very diffuse Gaussian matrix. To induce sparsity in S, a prior on B is selected such that it is sparse. E is a small Gaussian noise term and the prior on its variance small compared with the variance of Z * . Instead of solving Equation 101, it is more numerical convenient to solve the following problem:

A = U V T + Z * • B + E • (1 -B) (102)
To infer U ,V ,B,Z and E, Aicher [5] approximated the posterior distribution with a factorizable distribution. This a variational approach selects the distribution q closest to the posterior in the sense of Kullback-Leibler (KL) divergence. By parameterizing q, Aicher [5] converted the inference scheme back into an objective maximization problem. After selecting a distribution to approximate the posterior, the expectations of U ,V ,B,Z and E are taken to estimate them. Experimental results [5] show that FVBRPCA performs slightly better than RPCA solved via IALM [188], VBRPCA [17] and GoDec [389].

2.12.4 Bayesian RPCA with MoG noise(MoG-BRPCA) Zhao et al. [385] developed a generative RPCA model under the Bayesian framework by modeling data noise as a mixture of Gaussians (MoG). The MoG is a universal estimator to continuous distributions and thus MoG-BRPCA is able to fit a wide range of noises such as Laplacian, Gaussian, sparse noises and any combinations of them.

2.13 Approximated RPCA 2.13.1 "Go Decomposition" (GoDec) Zhou and Tao [389] proposed a randomized low-rank and sparse matrix decomposition called "Go Decomposition" (GoDec). GoDec estimates the low-rank part L and the sparse part S by using the same decomposition than SPCP [393]:

A = L + S + E ( 103
)
To solve the problem in Equation ( 103), GoDec alternatively assigns the lowrank approximation to A -S to L and the sparse approximation to A -L to S.

This approximated decomposition problem seeks to solve the minimization of the following decomposition error: min

L,S ||A -L -S|| 2 F subj rank(L) ≤ r, card(S) ≤ k. ( 104
)
Algorithm for solving GoDec: The optimization problem in Equation ( 104) is solved by alternatively solving the two following subproblems:

L t = arg min rank(L)≤e ||A -L -S t-1 || 2 F ( 105
) S t = arg min card(S)≤e ||A -L t -S|| 2 F ( 106
)
Although both subproblems have nonconvex constraints, their global solutions L t and S t exist. Indeed, these subproblems can be solved by updating L t via singular value hard thresholding of A -S t-1 and updating S t via entry-wise hard thresholding of A -L t , respectively as follows:

L t = r i=1 λ i U i V T i with SV D(A -S t-1 ) = U GV T ( 107
)
S t = P Ω (A -L t ) with Ω : (A -L t ) i,j∈Ω = 0 and ≥ (A -L t ) i,j∈ Ω , |Ω| ≥ k ( 108
)
where P Ω (.) is defined as the projection of the matrix on the observed entries following the sampling set Ω. The main computation time is due to the computation of the SVD for A -S t-1 in the updating L t sequence. To significantly reduce the time cost, Zhou and Tao [389] replaced the SVD by a Bilateral Random Projection(BRP) based low-rank approximation.

Semi-Soft GoDec

Zhou and Tao [389] proposed a Semi-Soft GoDec which adopts soft thresholding to the entries of S, instead of GoDec which imposes hard thresholding to both the singular values of the low-rank part L and the entries of the sparse part S. This improvement has two two main advantages: 1) the parameter k in constraint card(S) ≤ k is automatically determined by a soft-threshold τ , thus avoids the situation when k is chosen too large and some part of noise E is leaked into S;

2) the time cost is substantially smaller than the ordinary GoDec. For example, the background modeling experiments can be accomplished with a speed 4 times faster then ordinary GoDec, while the error is kept the same or even smaller.

The approximated decomposition problem seeks to solve the minimization of the following decomposition error: min

L,S ||A -L -S|| 2 F subj rank(L) ≤ r, card(S) ≤ τ. ( 109
)
where τ is the soft threshold. Chen et al. [START_REF] Chen | Incremental low-rank and sparse decomposition for compressing videos captured by fixed cameras[END_REF] proposed to use Semi-Soft GoDec for video coding in the existing standard codecs H.264/AVC and HEVC via background/foreground separation. For this, Chen et al. [START_REF] Chen | Incremental low-rank and sparse decomposition for compressing videos captured by fixed cameras[END_REF] developed an extension of the Semi-Soft GoDec that is able to perform LRSD on new matrix columns with a given low-rank structure, which is called incremental low-rank and sparse decomposition (ILRSD).

2.14 Sparse Additive Matrix Factorization (SAMF) Nakajima et al. [222][223] extented the original robust PCA [START_REF] Candes | Robust principal component analysis? International[END_REF] by proposing a unified view called Sparse Additive Matrix Factorization (SAMF). Instead of RPCA which only copes with element-wise sparsity (spiky noise) and low-rank sparsity (low-dimensional matrix), SAMF handles various types of sparse noise such as row-wise and column-wise sparsity. Thus, the decomposition is written as follows:

A = K k=0 S + E ( 110
)
where K is the number of sparse matrices. K = 2 in the original RPCA [START_REF] Candes | Robust principal component analysis? International[END_REF] in which the element-wise sparse term is added to the low-rank term. For background/foreground separation, the low-rank term and the element-wise sparse term capture the static background and the moving foreground, respectively. Nakajima et al. [222][223] relied on the natural assumption that a pixel segment which has similar intensity values in an image tends to belong to the same object. Thus, Nakajima et al. [222][223] adopted a segment-wise sparse term, where the matrix is constructed using a precomputed over-segmented. Experimental results [222][223] on the CAVIAR dataset [START_REF] Fisher | CAVIAR: context aware vision using image-based active recognition[END_REF] show that SAMF based on image segmentation (sSAMF) outperforms PCP via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] which correponds to 'LE .

Algorithm for solving SAMF: First, Nakajima et al. [222][223] reduced the partial SAMF problem to the standard MF problem, which can be solved analytically. Then, Nakajima et al. [222][223] derived an iterative algorithm called the mean update (MU) for the variational Bayesian approximation to SAMF, which gives the global optimal solution for a large subset of parameters in each step.

Variational Bayesian Sparse Estimator (VBSE)

Chen et al. [57][58][56] proposed a generalization of the original RPCA [START_REF] Candes | Robust principal component analysis? International[END_REF], where a linear transformation through the use of a known measurement matrix, is applied to the outlier corrupted data. The aim is to estimate the outlier amplitudes given the transformed observation. This approach called variational Bayesian Sparse Estimator (VBSE) can achieved background/foreground separation in blurred and noisy video sequences.Thus, the decomposition is written as follows:

A = L + RS + E (111)
where R models the linear transformation performed on the data. The aim is to obtain accurate estimates for the sparse term S and the low-rank term L, given the noise corrupted observation A. Although S is sparse, the multiplication with a wide matrix R has an effect of compression, and hence the product RE is not necessarily sparse. Then, Chen et al. [57][58] modeled the lowk-rank part as follows:

||L|| * = min U,V 1 2 ||U || 2 F + ||V || 2 F subj L = U V T (112)
With these relaxation and parametrization, Chen et al. [57][58] obtained the following optimization problem:

min U,V,S 1 2 ||A -U V T -RS|| 2 F + λ * (||U || 2 F + ||V || 2 F ) + λ 1 ||E|| l1 subj L = U V T (113)
where λ * and λ 1 are regularization parameters. To enforce column sparsity in U and V , the columns of U and V are modeled with Gaussian priors of precision. Then, Chen et al. [57][58] incorporated conjugate Gamma hyperprior on the precisions. The sparse part S is modeled by setting the entries be independent of each other, and their amplitudes are modeled by zero-mean Gaussian distributions with independent precisions. For the noise part E, Gaussian priors with zero mean are used to model the dense observation noise. By combining these different stages in a hierarchical Bayesian model, a joint distribution of the observation and all the unknowns is expressed as follows:

ρ(A, U, V, S, γ, α, β) (114)
where γ and α are hyperparameters and β is the noise precision. To solve VBSE, Chen et al. [57][58] used an an approximate Bayesian inference. Experimental results [57][58] on the CAVIAR dataset [START_REF] Fisher | CAVIAR: context aware vision using image-based active recognition[END_REF] show that VBSE outperforms PCP solved via APG [189] and PCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF].

3 Robust Non-negative Matrix Factorization

Non-negative matrix factorization (NMF) approximates a non-negative matrix A by a product of two non-negative low-rank factor matrices W and H. Table 5 shows an overview of the different robust NMF decompositions. Their corresponding solvers as well as their complexity can be seen in Table 12.

3.1 Manhattan Non-negative Matrix Factorization (MahNMF)

Guan et al. [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF] proposed a robust non-negative matrix factorization when the noise distribution is heavy tailed. The method called Manhattan NMF (MahNMF) minimizes the Manhattan distance between A and W T H for modeling the heavy tailed Laplacian noise. Thus, Guan et al. [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF] minimized the Manhattan distance between an m × n-dimensional non-negative matrix A and W T H as follows:

min W ≥0,H≥0 f (W, H) = ||A -W T H|| M (115)
where ||.|| M is the Manhattan distance which is equals the summation of the l 1 norm, and the reduced dimensionality r satisfies that r ≪ min(m, n). Since W and H are low-rank matrices, MahNMF actually estimates the nonnegative low-rank part, i.e., W T H, and the sparse part, i.e., A -W T H, of a non-negative matrix A. MahNMF performs effectively and robustly when data are contaminated by outliers because it benefits from both the modeling ability of Laplace distribution to the heavy tailed behavior of noise and the robust recovery capability of the sparse and low-rank decomposition. Experimental results [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF] on the I2R dataset [177] show that MahNMF gives similar visual results than PCP via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and GoDec [389].

Algorithms for solving MahNMF: Two fast optimization algorithms for Mah-NMF were developed by Guan et al. [99] and called the rank-one residual iteration (RRI) method and Nesterov's smoothing method, respectively. By approximating the residual matrix by the outer product of one row of W and one row of H, the RRI method iteratively updates each variable of W and H in a closed form solution. RRI is neither scalable to large scale matrices nor flexible enough to optimize all MahNMF extensions. As the objective functions of MahNMF are neither convex nor smooth, Guan et al. [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF] proposed a Nesterov's smoothing method to recursively optimize one factor matrix with another matrix fixed. By setting the smoothing parameter inversely proportional to the iteration number, the approximation accuracy is improved iteratively.

Near-separable Non-negative Matrix Factorization (NS-NMF)

Promising robust NMF approaches have emerged recently based under the assumption that the data matrix satisfies a separability condition which enables the NMF problem to be solved efficiently and exactly. Under this assumption, the data matrix A is said to be r-separable if all columns of X are contained in the conical hull generated by a subset of r columns of A. In other words, if A has a factorization W H then the separability assumption states that the columns of W are present in A at positions given by an unknown index set B of size r. Equivalently, the corresponding columns of the right factor matrix H constitute the r×r identity matrix, i.e., H B = I. These columns indexed by B are called anchor columns. In this framework, Kumar et al. [164] proposed a family of conical hull finding procedures called Xray for near-separable NMF (NS-NMF) problems with Frobenius norm loss. The minimization problem can be formulated as follows:

min AB ≥0,H≥0 ||A -A B H|| 2 F subj A B ≥ 0, H ≥ 0 (116)
Geometrically, Xray finds anchor columns one after the other, incrementally expanding the cone and using exterior columns to locate the next anchor. Xray present several advantages for background/foreground separation: (1) it requires no more than r iterations each of which is parallelizable, (2) it empirically robust to noise, (3) it admits efficient model selection, and (4) it does not require normalizations or preprocessing needed in other methods. However, in the presence of outliers or different noise characteristics, the use of Frobenius norm approximations is not optimal. In this context, Kumar and Sindhwani [163] improved Xray to provide robust factorizations with respect to 1 1 loss, and approximations with respect to the family of Bregman divergences. In the case of background/foreground separation, it is natural to seek a low-rank background matrix L that minimizes ||A -L|| where A is the frame-by-pixel video matrix, and the l 1 loss imposes a sparsity prior on the residual foreground. Instead of low-rank approximations in Frobenius or spectral norms, there is not a SVD-like tractable solution. For this, Kumar et al. [163] recovered tractability by imposing the separable NMF assumption on the background matrix. This implies that the variability of pixels across the frames can be considered as observed variability in a small set of anchor pixels. Under a more restrictive setting, it is equivalent to median filtering on the video frames, while a full near-separable NMF model imparts more degrees of freedom to model the background. The minimization problem can be formulated as follows:

min AB ≥0,H≥0 ||A -A B H|| l1 subj A B ≥ 0, H ≥ 0 ( 117
)
where A B are the columns of A indexed by set B ⊂ 1, 2, ..., n. Experimental results [163] on the I2R dataset [177] show that RobustXray outperforms the robust NMF (local search) which minimizes min and PCP via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF].

Algorithms for solving XRay-l 2 and RobustXray: Algorithms of Kumar et al. [164] are not suitable for noise distributions other than Gaussian. The algorithm for RobustXray proceeds by identifying one anchor column in each iteration and adding it to the current set of anchors, thus expanding the cone generated by anchors. Each iteration consists of two steps: [START_REF] Abdel-Hakim | FRPCA: fast robust principal component analysis[END_REF] anchor selection step that finds the column of A to be added as an anchor, and (2) a projection step where all data points (columns of A ) are projected to the current cone in terms of minimizing the l 1 norm.

3.3 Robust Asymmetric Non-negative Matrix Factorization (RANMF) Woo and Park [334] proposed a formulation called l∞-norm based robust asymmetric nonnegative matrix factorization (RANMF) for the grouped outliers and low nonnegative rank separation problems. The main advantage of RANMF is that the denseness of the low nonnegative rank factor matrices can be controlled. To control distinguishability of the column vectors in the low nonnegative rank factor matrices for stable basis, Woo and Park [334] imposed asymmetric constraints, i.e., denseness condition on the coefficient factor matrix only. As a byproduct, a wellconditioned basis factor matrix is obtained. Compared to the nuclear norm based low-rank enforcing models, RANMF is not sensitive to the nonnegative rank constraint parameter due to the soft regularization method. Thus, the decomposition is made as follows:

A = L + S = W ΛH + S (118)
where L = W ΛH is low nonnegative rank matrix and S contains the grouped outliers. Λ is a diagonal matrix with Λ ii = λ i which is be considered as an asymmetric singular value of L. Since Λ is subsumed into W or H, Equation 118 is a typical nonnegative matrix factorization (NMF) formulation and its corresponding minimization problem can be written as follows:

min L,S,Φ α 2 ||A-L-S|| 2 F +Φ(S)+βΨ(S, Φ)+γT V (Q(Φ)) subj Rank(L) ≤ τ, 0 < L < b L (119
) where Φ(.) is a sparsity enforcing function, such as lp-norm (0 < p < 1) or logfunction. Q is the reshaping operator from 2D to 3D, T V (.) is the 3D Total Variation and b L = 255 for image data. Rank(L) is a low nonnegative rank enforcing function such as nonnegative nuclear norm, and τ ≥ 0 is a rank nonnegative constraint parameter. As grouped outliers is not foreground mask and can be very noisy, the TV appears as an additional denoising/segmentation process to detect foreground mask. To solve Equation 119, Woo and Park [334] developed a Soft Regularized Asymmetric alternating Minimization (SRAM) algorithm.

Experimental results [334] on the I2R dataset [177] show that RANMF outperforms PCP via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and DECOLOR [391].

Robust Matrix Completion

Robust matrix completion RMC, also called RPCA plus matrix completion problem can also be used for background/foreground separation. Althought RPCA via principal component pursuit [START_REF] Candes | Robust principal component analysis? International[END_REF] can be considered as RMC using l 1 -norm loss function, the main difference lies in that in RMC problems the support of missing entries is given, whereas in RPCA corrupted entries are never known [335][360]. From a statistical learning viewpoint, RPCA is a typical unsupervised learning problem while the RMC can be interpreted as a supervised learning problem [360]. Table 5 shows an overview of the different RMC decompositions. Their corresponding solvers as well as their complexity can be seen in Table 12.

lσ-norm loss function (RMC-lσ)

RPCA via principal component pursuit can be considered as RMC using l 1 -norm loss function. Following this idea, Yang et al. [360] proposed a nonconvex relaxation approach to the matrix completion problems when the entries are contaminated by non-Gaussian noise or outliers. Based on a nonconvex lσ loss function, Yang et al. [360] developed a rank constrained as well as a nuclear norm regularized model. The nuclear norm heuristic model is formulated in the following form: min

S∈R m×n λ||L|| * + lσ(L) ( 120
)
where λ is a regularization parameter and the data fitting risk lσ(L) is given by:

lσ(L) = σ 2 2 (i,j)∈Ω (1 -exp(-(L ij -A ij ) 2 /σ 2 )) (121)
Experimental results [360] show that RMC-lσ performs slightly better than PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] because details of the background images are recovered well, whereas PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] does not seem to perform as well as RMC-lσ where some details of the background are added to the foreground. It can be also observed that none of the two methods can recover the missing entries in the foreground. Furthermore, RMC-lσ is more than 3 times faster than PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF].

Algorithms for solving RMC-lσ-IHT: lσ-IHT leads to computational difficulty due to its nonconvexity. To solve this problem, Yang et al. [360] developed two algorithms based on iterative soft thresholding (IST) and iterative hard thresholding (IHT). These two algorithms are called lσ-IST and lσ-IHT. By verifying the Lipschitz continuity of the gradient of the datafitting risk, lσ-IST and lσ-IHT converge. Under proper conditions, the recoverability as well as the linear convergence rate results are obtained. Only RMC-lσ-IHT was tested on the I2R dataset [177].

Robust Bilateral Factorization (RMC-RBF)

For RMC, Shang et al. [275][273] proposed a scalable and provable structured low-rank matrix factorization method to recover the low-rank plus sparse matrices from missing and grossly corrupted data. Thus, a Robust Bilinear Factorization (RBF) method recovered the low-rank plus sparse matrices from incomplete and/or corrupted data, or a small set of linear measurements. The decomposition is the followong one:

A = L + S = U V T + S (122)
where U and V are two matrices of compatible dimensions, where U has orthogonal columns. Then, the corresponding minimization problemn is formulated as follows:

min U,V,S λ||V || * + ||P Ω (S)|| l1 + subj P Ω (A) = P Ω (UV T + S), U T U = I (123)
where λ ≥ 0 is the regularization parameter, ||V || * is the nuclear norm of the lowrank matrix V ∈ R m×n , S ∈ R m×n is the sparse error matrix. Ω is the index set of observed entries and P Ω (.) is the projection operator onto that subspace. RBF not only takes into account the fact that the observation is contaminated by additive outliers or missing data, but can also identify both low-rank and sparse noisy components from incomplete and grossly corrupted measurements. So, Shang et al.

[275] developed two small-scale matrix nuclear norm regularized bilinear structured factorization models for RMC problems, in which repetitively calculating SVD of a large-scale matrix is replaced by updating two much smaller factor matrices. Then, Shang et al. [275][273] applied the alternating direction method of multipliers (ADMM) to efficiently solve the RMC problems. Experimental results show that RBF gives similar visual results than RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and GRASTA [119] but RBF is more than 3 times faster than GRASTA [119] and more than 2 times faster than RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF].

Matrix Factorization (RMC-MF)

The general RMC problem aims to simultaneously recover both low-rank and sparse components from incomplete and grossly corrupted observations via the following convex optimization problem: min

L,S ||L|| * + λ||S|| l1 subj P Ω (L + S) = P Ω (A) (124)
where P Ω (A) is defined as the projection of the matrix A on the observed entries Ω: P Ω (A ij ) = A ij if (i, j) ∈ Ω and P Ω (A ij ) = 0 otherwise. From Equation 124, Shang et al. [274] find that the optimal solution E Ω C = 0 where Ω C is the complement of Ω and corresponds to the index set of the unobserved entries. So, the RMC problem is equivalent to the following convex optimization problem: min

L,S ||L|| * + λ||P Ω (S)|| l1 subj P Ω (L + S) = P Ω (A), E Ω C = 0 (125)
To efficiently solve the RMC problem and avoid introducing some auxiliary variables, Shang et al. [274] assumed that the constraint with a linear projection operator P Ω can be simplified into A = L + S. To further improve the efficiency of this convex model and the scalability of handling large data sets, Shang et al. [274] proposed a scalable non-convex model in which the desired low-rank matrix L is factorized into two much smaller matrices U ∈ R m×d and V ∈ R d×n where d is an upper bound for the rank of the matrix L, i.e., d ≥ r = rank(L). Thus, the decomposition problem is formulated as follows:

A = U V T + S (126)
Finally, the RMC with Matrix Factorization (RMC-MF) problem is equivalent to the following convex optimization problem:

min U,V,S ||V || * + λ||P Ω (S)|| l1 subj A = U V T + S, U T U = I (127)
Experimental results show that RMC (convex formulation) and RMC-MF (nonconvex formulation) are slightly better than that of GRASTA [119] and UNN-BF [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF]. The theoretical reason for the unsatisfactory performance of the l 1 -penalty is that the irrepresentable condition is not met. Hence, RMC-MF incorporating with matrix factorization is more accurate in recovering the low-rank matrix than RMC (convex formulation). Furthermore, RMC-MF is more than 7 times faster than RMC (convex formulation), more than 4 times faster than GRASTA [119], and more than 2 times faster than UNN-BF [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF].

Algorithms for solving RMC and RMC-MF: Shang et al. [274] developed two efficient alternating direction augmented Lagrangian (ADAL) solvers for solving the convex model and the non-convex model, respectively. For the convex problem, the running time of the corresponding algorithm is dominated by that of performing SVD on the matrix of size m × n. For the non-convex problem, the corresponding algorithm performs SVD on much smaller matrices of sizes m × d and d×n, and some matrix multiplications. Hence, the total time complexity of the algorithm for the convex RMC and the algorithm for the non-convex are O(tmn 2 ) and O(t(d 2 m + mnd)) with (d ≪ n < m), respectively, where t is the number of iterations.

Factorized Robust Matrix Completion (FRMC)

Mansour and Vetro [212] developed a factorized robust matrix completion (FRMC) algorithm with global motion compensation to solve the background/foreground separation problem. Since the main drawbacks in Equation 22is that it requires the computation of full or partial singular value decompositions of L in every iteration of the algorithm, Mansour and Vetro [212] adopted a surrogate for the nuclear norm of a rank-r matrix L defined by the following factorization:

||L|| * = min LL∈R m,r ,LR∈R n,r 1 2 (||L L || 2 F + ||L R || 2 F ) subj L L L T R = L (128)
This nuclear norm surrogate can be used in standard nuclear norm minimization algorithms that scale to very large matrix completion problems. Moreover, it was shown that when the factors L L and L R have a rank greater than or equal to the true rank of L. Then, each subproblem in the FRMC algorithm is a Lasso problem that Mansour and Vetro [212] solved using spectral projected gradient iterations. FRMC was developed in batch mode and online mode. The FRMC algorithm in online mode completes the recovery 7 to 9 times faster than GRASTA [119] and results in a comparable separation quality.

4.5 Motion-Assisted Matrix Completion (MAMC) Yang et al. [351][365] proposed a motion-assisted matrix completion (MAMC) model for foreground-background separation. Thus, a dense motion field is estimated for each frame, and mapped into a weighting matrix W 3 which indicates the likelihood that each pixel belongs to the background as follows: min

L,S ||L|| * + λ||S|| l1 subj W 3 • A = W 3 • (L + S) (129)
where • is denotes the element-wise multiplication of two matrices and W 3 is constructed from motion information. By incorporating this information, areas dominated by slowly-moving objects are suppressed while the background that appears at only a few frames has more chances to be recovered in the foreground detection results. The influence of light conditions, camouflages, and dynamic backgrounds can also be decreased.

In addition, Yang et al. [351][365] extended MAMC to a robust MAMC model (RMAMC) which is robust to noise for practical applications as follows:

min L,S,E ||L|| * + λ||S|| l1 + γ||E|| 2 F subj W 3 • A = W 3 • (L + S + E) ( 130
)
where γ is a positive constant and E is the matrix which contains the noise. Yang et al. [351] adapted the ALM algorithm [188] to solve MAMC and RMAMC. Experimental results [351][365] on several datasets show that RMAMC outperforms RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF].

Robust Subspace Recovery

In this category are the robust decompositions other than RPCA and RNMF decompositions. Table 6 shows an overview of the different RSR decompositions.

Their corresponding solvers as well as their complexity can be seen in Table 12.

Robust Subspace Recovery via Bi-Sparsity (RoSuRe)

High dimensional data is distributed in a union of low dimensional subspaces in sparse models but the underlying structure may be affected by sparse errors and/or outliers. To solve this problem, Bian and Krim [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] proposed a bi-sparsity model and provided an algorithm to recover the union of subspaces in presence of sparse corruption. This, the proposed decomposition is the following one:

A = L + S = LW + S ( 131
)
where W 4 is sparse matrix. Then, the corresponding minimization problem is formulated as follows:

min

W,S ||W 4 || l1 + λ||S|| l1 subj A = L + S, L = LW 4 , W 4,(ii) = 0, ∀i (132)
Experimental results [24][23] on static and moving camera cases show the ability of RoSuRe to separate foreground from background in both cases. More interestingly, the sparse coefficient matrix W 4 gives information about the relations among data points, which potentially may be used to cluster data into individual clusters. Indeed, for each column of the coefficient matrix W 4 , the nonzero entries appear periodically. In considering the periodic motion of the camera, every frame is mainly represented by the frames when the camera is in a similar position, i.e. a similar background, with the foreground moving objects as outliers. After permuting the rows and columns of W 4 according to the position of cameras, a block-diagonal structure appears, where images with similar backgrounds are clustered as one subspace.

Algorithms for solving RoSuRe: Bian and Krim [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] developed an algorithm via Bi-Sparsity Pursuit based on linearized ADMM [190]. Practically, Bian and Krim [START_REF] Bian | Robust subspace recovery via bi-sparsity pursuit[END_REF] pursued the sparsity of S and W 4 alternatively until convergence. Besides the effectiveness of ADMM on l 1 minimization problems, the augmented Lagrange multiplier (ALM) method can address the non-convexity of Equation 132. It hence follows that with a sufficiently large augmented Lagrange multiplier, the global optimizer is approximated by solving the dual problem.

5.2 Robust Orthonomal Subspace Learning (ROSL) Shu et al. [286][287] presented a computationally efficient low-rank recovery method, called as Robust Orthonormal Subspace Learning (ROSL). ROSL speeds the rankminimization of a matrix L by imposing the group sparsity of its coefficients α under orthonormal subspace spanned by orthonormal bases D. Its underlying idea is that, given the subspace representation L = Dα, the rank of L is upper bounded by the number of non-zero rows of α, that is ||α|| row-0 . ROSL can be regarded as a non-convex relaxation of RPCA by replacing nuclear norm with this rank heuristic [286]. So, ROSL involved the following decomposition:

A = Dα + S (133)
Thus, ROSL recovers the low-rank matrix L from A by minimizing the number of non-zero rows of α, and the sparsity of S as follows: min S,D,α

||α|| row-0 + λ||S|| l0 subj Dα + S = A, D T D = I k , ∀i (134)
where the subspace bases D = U , the coefficients α = SV T and L = U SV T obtained by SVD. As the sparsity-inducing l 1 -norm is an acceptable surrogate for the sparsity measure with l 0 -norm, Shu et al. [287] reformulated ROSL as the following non-convex optimization problem: min S,D,α

||α|| row-1 + λ||S|| l1 subj Dα + S = A, D T D = I k , ∀i (135)
where the row -1-norm is defined as

||α|| row-1 = k i ||α i || l 2
Experimental results [286][287] on the I2R dataset [177] show that the recovery accuracy and efficiency of ROSL is slightly better than PCP solved via EALM [188], PCP solved via IALM [188], PCP solved by Random Projection [221] and PCP solved by LMaFit [281]. Furthermore, ROSL is more than 10 times faster than RPCA solved via IALM [188].

Algorithms for solving ROSL: Shu et al. [287] presented an efficient sparse coding algorithm to minimize this rank measure and recoverthe low-rank matrix at quadratic complexity of the matrix size. This, ROSL is solved using inexact ADM (Alternating Direction Method) at the higher scale and inexact BCD (Block Coordinate Descent) at the lower scale. This solver is called inexact ADM/BCD. Finally, Shu et al. [287] developed a random sampling algorithm to further speed up ROSL such that its accelerated version (ROSL+) has linear complexity with respect to the matrix size. ROSL+ is more than 92 times faster than RPCA solved via IALM [188].

Robust Orthogonal Complement Principal Component Analysis (ROC-PCA)

She et al. [277] proposed a robust orthogonal complement principal component analysis (ROC-PCA). The aim is to deal with orthogonal outliers that are not necessarily apparent in the original observation space but could affect the principal subspace estimation. For this, She et al. [277] introduced a projected mean-shift decomposition as follows:

AV ⊥ = L + S + E (136)
where V ⊥ is n × m matrix verifying V T ⊥ V ⊥ = I and characterizes the subspace orthorgonal to the rank-r principal component subspace. AV ⊥ gives the coordinates after projecting the data onto the orthogonal components subspace and it decomposed into three parts: mean L, outlier S and noise E. The corresponding minimization problem of Equation 136 is formulated as follows: min

V ⊥ ,µ,S 1 2 ||AV ⊥ -L -S|| 2 F + ij P (||s ij || l 2 ; λ i ) subj V T ⊥ V ⊥ = I ( 137
)
where I is the identity matrix and s ij is the i th row vector of Subspace tracking addresses the problem when new observations come progressively as in online streaming application. The algorithm cannot store all the input data in memory. Thus, the new incoming observations need to be processed and then discarded. The involved subspaces can have low-rank and/or sparse structures as in the previous decomposition frameworks. Table 6 shows an overview of the different RMC decompositions. Their corresponding solvers as well as their complexity can be seen in Table 13.

6.1 Grassmannian Subspace Tracking (GRASTA) He et al. [119][120] proposed a Grassmannian robust adaptive subspace tracking algorithm (GRASTA). This algorithm uses a robust l 1 -norm cost function in order to estimate and track non-stationary subspaces when the streaming data vectors are corrupted with outliers. This problem is solved via an efficient Grassmannian augmented Lagrangian Alternating Direction Method. Let denote the evolving subspace of R n×m as S t at time t with its dimension d that is supposed to be much smaller than m and n. Let the columns of an mn × d matrix U t be orthonormal and span S t . Tracking the subspace S t is equivalent to estimating U t at each time step t. The observed vector data A t is assumed to be generated at each time step t as follows:

A t = U t w t + S t + E t ( 138
)
where U t w t = L t posses a low-rank structure, and w t is a d × 1 weight vector. The orthonormal columns of Ut span the low-rank subspace of the images. The set of all subspaces of R n dimension d is called the Grassmannian, which is a compact Riemannian manifold and is denoted by G(d, n). S t is the n×1 sparse outlier vector whose nonzero entries may be arbitrarily large, and S t models foreground pixels in the background/foreground separation. E t is the n × 1 zero-mean Gaussian white noise with small variance. Then, He et al. [119][120] subsampled A t on the index set Ω t ⊂ {1, ..., n}. So, only a small subset of entries of A t are kept. U Ωt is the submatrix of U t consisting of the rows indexed by Ω t . For a vector A t ∈ R n , A Ωt is the vector in R |Ωt| whose entries are those of A t indexed by Ω t . To quantify the subspace error when the data are incomplete and corrupted, GRASTA uses the l 1 -norm to measure the subspace error from the subspace spanned by the column of U t to the observed vector A Ωt :

F (S; t) = min w ||U Ωt w -A Ωt || l1 (139)
If U Ωt is known or can be estimated, this l 1 -minimization problem can be solved by Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distribute optimization and statistical learning via the alternating direction method of multipliers[END_REF]. According to ADM, the Equation ( 139) is equivalent to the following problem by introducing a sparse outlier vector S ∈ R |Ωt| : min

U,w,S ||S|| l1 subj A Ωt = U Ωt w + S, , U ∈ G(d, n) ( 140
)
This problem is not convex but it offers much more computationally efficient solutions. Experimental results [119] on the I2R dataset [177] show that GRASTA show more robustness than RPCA solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and ReProCS [243] with less time requirement.

Algorithm for solving GRASTA: The problem in Equation ( 140) is solved via the augmented Lagrangian function:

L(S, w, y) = ||S|| 1 + y T (U Ωt w + S -A Ωt ) + ρ 2 ||U Ωt w + S -A Ωt || 2 l2 ( 141
)
where y is the dual vector. The unknowns are S, w, y and U . If U is fixed, the triple (S, w, y) is solved by ADMM [START_REF] Boyd | Distribute optimization and statistical learning via the alternating direction method of multipliers[END_REF] and if the triple (S, w, y) is fixed, U is estimate by Grassmannian geodesic gradient descent [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF]. GRASTA is then composed by this alternating approach. The total computational cost of GRASTA is

O(|Ω| d 3 + T d |Ω| + nd 2 )
where |Ω| is the number of samples per vector used, d is the dimension of the subspace, n is the ambient dimension, and T is the number of ADMM iterations.

||S|| l1 subj A Ωt • τ = U Ωt w + S , U ∈ G(d, n) (142)
where τ are the transformations. He et al. [122][123] developed batch mode and online mode algorithms. For batch mode, U is the iteratively learned aligned subspace in each iteration; while for online mode, U is a collection of subspaces which are used for approximating the nonlinear transform, and they are updated iteratively for each video frame. To solve t-GRASTA, He et al. [122][123] used a ADDM solver suitable for the locally linearized problem.Experimental results [122][123] on sequences with simulating camera jitters show that t-GRASTA outperforms RASL [239] and GRASTA [119]. Furthermore, t-GRASTA is four faster than state-of-theart algorithms and has half the memory requirement.

6.3 Grassmannian Adaptive Stochastic Gradient with l 2,1 -norm (GASG21)

In the presence of column outliers corruption, He and Zhang [124] formulated the Grassmannian Adaptive Stochastic Gradient for l 2,1 -norm minimization (GASG21). Moreover, the classical matrix l 2,1 -norm minimization problem is formulated in its stochastic programming coun-terpart.

min w ||U w -A|| l2,1 - m j=1 ||U j w j -A j || l2 subj U ∈ G(d, n) (143)
The l 2,1 -norm minimization is well suitable for column outliers corruption. For inliers which can be well represented by the subspace, the residues are small. For outliers which can not be fitted into the subspace, the residues are large. Then Equation 143 means that HHe and Zhang [124] are optimizing U which can best fit inliers to reduce the sum of l 2 fit residues. To solve GAS21, He and Zhang [124] solved the l 2,1 -norm minimization by stochastic gradient descent (SGD).

Experimental results [124] show that GAS21 outperfoms slightly OP [342] with less time computation.

6.4

Lp-norm Robust Online Subspace Tracking (pROST)

GRASTA performs a gradient descent on the Grassmannian and aims at optimizing an l 1 -cost function to mitigate the effects of heavy outliers in the subspace tracking stage. He et al. [119] overcame the nondifferentiability of the l 1 -norm by formulating an augmented Lagrangian optimization problem at the cost of doubling the number of unknown parameters. Although the l 1 -norm leads to favorably conditioned optimization problems it is well-known that penalizing with nonconvex l 0 -surrogates allows reconstruction even in the case when l 1 -based methods fail. Therefore, Hage and Kleinsteuber [115][268] proposed a framework that combines the advantages of Grassmannian optimization with non-convex sparsity measures. This method called pROST focuses primarily on reconstructing and tracking the underlying subspace and can operate on both fully and incompletely observed data sets. The involved decomposition is the same than in GRASTA (See Equation 138) but the minimization problem is formulated as follows:

min

Rank(L)≤k ||U W -A|| Lp (144)
where U W = L. ||.|| Lp is a smoothed and weighted Lp-quasi-norm cost function to achieve robustness against outliers, and it is defined as follows:

Lp(X) = m i=1 (x 2 i + µ) 1 p 0 < p < 1 (145)
Spatial and temporal proximity of pixels leads that corresponding pixels in consecutive frames are likely to have the same label. Thus, Hage and Kleinsteuber [115][268] used this knowledge to increase the robustness of the residual cost. The contribution of labeled foreground pixels to the overall penalty is reduced by introducing additional pixel weights w i . If the pixel was previously labeled a foreground pixel and is then likely to remain an outlier in the current frame, the weight is small to avoid foreground objects compromising the background. If the pixel was labeled a background pixel, the weight is equal to one for an adaptive model maintenance. Thus, the weighted smoothed Lp-quasi-norm cost function is defined as follows:

Lp(X) = m i=1 w i (x 2 i + µ) 1 p 0 < p < 1 (146)
Then, an alternating online optimization framework for estimating the subspace makes the algorithm suitable for online subspace tracking. In contrast to GRASTA, the method presented directly optimizes the cost function and thus operates with less than half the number of unknowns.

pROST can be applied in real-time background/foreground separation and makes use of the spatio-temporal dependencies between pixel labels. This leads to especially good performance in presence of bootstrapping, large foreground objects (which often arise in RPCA-based methods) and jittery cameras. Experimental results [268] on the ChangeDetection.net dataset [START_REF] Goyette | changedetection.net: A new change detection benchmark dataset[END_REF] confirm that the proposed method can cope with more outliers and with an underlying matrix of higher rank than GRASTA. Particularly, pROST outperforms GRASTA in the case of multimodal backgrounds.

Algorithm for solving pROST: Hage and Kleinsteuber [115] [268] used a Conjugate Gradient (CG) type algorithm on the Grassmannian for solving the individual minimization tasks. Like all optimization methods on the Grassmannian, the algorithm allows to upper-bound the dimension of the underlying subspace and easily extends to the problem of robustly tracking this subspace.

6.5 Grassmannian Online Subspace Updates with Structured-sparsity (GOSUS) Xu et al. [343] studied the problem of online subspace learning when sequential observations involves structured perturbations. As the observations are an unknown mixture of two components presented to the model sequentially, if no additional constraints is imposed on the residual, it often corresponds to noise terms in the signal which were unaccounted for by the main effect. To address this problem, Xu et al. [343] imposed structural contiguity, which has the effect of leveraging the secondary terms as a covariate that helps the estimation of the subspace itself, instead of merely serving as a noise residual. min

U T U=I d ,W,S l i=1 µ i ||D i S|| l2 + λ 2 ||U W + S -A|| 2 l2 ( 147
)
where U W = L posses a low-rank structure, S posses a sparse structure, and D i is diagonal matrix which denotes a group i. Each diagonal element of D i corresponds to the presence/absence of a pixel in the i t h group. So, D ij is equal to one if pixel j is in group i, and it is equal to zero otherwise. Thus, the term l i=1 µ i ||D i S|| l2 is penalty function where µ i is the weight for group i and l is the number of such groups. This group sparsity function has a mixed norm structure. The inner norm is either l 2 forcing pixels in the corresponding group to be similarly weighted, and the outer norm is l 1 which encourages sparsity, that is.only few groups are selected. The corresponding online estimation procedure for Equation 147 is written as an approximate optimization process on a Grassmannian, called Grassmannian Online Subspace Updates with Structured-sparsity (GOSUS). GOSUS is solvable via an alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distribute optimization and statistical learning via the alternating direction method of multipliers[END_REF] applied in a blockwise manner. Experimental results [343] on the Wallflower dataset and the I2R dataset [177] show that GOSUS outperforms slightly RPCA solved via IALM [188], GRASTA[119] and BRMF [366].

6.6 Fast Adaptive Robust Subspace Tracking (FARST)

In spite of these good properties, the global convergence of GRASTA is note proved as developed in [4]. Empirically, it is slowly adapted to background change or not adapted to dynamic background. On the contrary, FARST [4] has an optimal global convergence while sharing some favorable properties with GRASTA.

FARST shares the procedure of separating frames into background and foreground with GRASTA, but it uses a recursive least square algorithm for subspace tracking, which makes it fast adapted to background change and dynamics. Every time a video frame streams in, two alternating procedures are repeatedly done. First, basis images are updated by a recursive least square algorithm. Secondly, foreground images are extracted by solving the l 1 -minimization problem. Furthermore, FARST is an online algorithm fast adapted to background change. Results [4] show that FARST outperforms GRASTA [119] and PRMF [321] in the presence of dynamic backgrounds. FARST is solvable via an alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distribute optimization and statistical learning via the alternating direction method of multipliers[END_REF].

Robust Low Rank Minimization

Low-rank minimization is a minimization problem involving a cost function which measures the fit between a given data matrix A and an approximating low-rank matrix L. Table 18 shows an overview of the different Robust Low Rank Minimization (RLRM) methods and the corresponding solvers. Furthermore, the complexity of the solverts can be seen in Table 13.

7.1 LRM with contiguous outliers detection (DECOLOR) Zhou et al. [391] proposed a formulation of outlier detection in the low-rank representation, in which the outlier support and the low-rank matrix are estimated. This method is called Detecting Contiguous Outlier detection in the Low-rank Representation (DECOLOR). The decomposition involves the same model than SPCP in Equation [START_REF] Cao | Total variation regularized RPCA for irregularly moving object detection under dynamic background[END_REF], that is A = L + S + E. So, the following formulation is proposed to achieve the decomposition:

min L,S αrank(L) + β||S|| l0 + 1 2 ||A -L -S|| 2 F subj rank(L) ≤ r (148) Categories Authors -Dates subj rank(L) ≤ r (152)
The binary states of entries in foreground support F are modeled by a Markov Random Field because the foreground objects are contiguous pieces with relatively small size. Based on the first order MRFs, the following regularizer on F is used:

||Cvec(F )|| l 1 = (ij,kl)∈N F ij -F kl ( 153
)
where C is the node-edge incidence matrix of a graph G with m×n nodes, and N is the set of all pairs of adjacent nodes in G. ij and kl denote the nodes corresponding to S ij and S kl , respectively. Then, the following energy function is obtained by relaxing rank(L) with the nuclear norm and adding the continuity constraint on F : min

L,S α||L|| * + β||F || l1 + 1 2 ||P F (A -L)|| 2 F + γ||Cvec(F )|| l 1 (154)
where α > 0 is a parameter which controls the complexity of the background model. S is recovered by S = P F (L -S).

Experimental results [391] on the I2R dataset [177] show that DECOLOR outperforms PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and two conventional models which are the mean model and the MOG model [298] but the main drawback of DECOLOR is its prohibitive computation time. However, the code called DECOLOR33 is provided

Algorithm for solving DECOLOR: The objective function is non-convex and its include both continuous and discrete variable. Zhou et al. [307] adopted an alternating algorithm that separates the energy minimization over B and F into two steps. B-step is a convex optimization problem and F -step is a combinatorial optimization problem. The optimal B is computed efficiently by the SOFT-IMPUTE [215] algorithm and the outlier support S is estimated in polynomial time using graph-cut [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][161].

7.2 LRM with Direct Robust Matrix Factorization (DRMF) Xiong et al. [341] proposed a direct robust matrix factorization (DRMF) assuming that a small portion of the matrix A has been corrupted by some arbitrary outliers.

The aim is to get a reliable estimation of the true low-rank structure of this matrix and to identify the outliers. To achieve this, the outliers are excluded from the model estimation. The decomposition involves the same model than PCP in Equation (20), that is A = L + S. The direct formulation of DRMF is written as follows: min

L,S ||A -S -L|| F subj rank(L) ≤ r, ||S|| l0 ≤ p ( 155
)
where L is the low-rank approximation, r is the rank, S is the matrix of outliers, and p is the maximal number of entries that can be ignored as outliers. Comparing DRMF to the conventional LRM, the difference is that the outliers S can be excluded from the low-rank approximation, as long as the number of outliers is not too large, that is, S is sufficiently sparse. By excluding the outliers from the low-rank approximation, Xiong et al. [341] ensured the reliability of the estimated low-rank structure. PCP [START_REF] Candes | Robust principal component analysis? International[END_REF] is the convex relaxation of DRMF.

Experimental results [341] on the I2R dataset [177] show that DRMF outperforms SVD [START_REF] Golub | Matrix computation[END_REF], PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and SPCP [393]. Furthermore, the code called DRMF 34 is provided.

Algorithm for solving DRMF: Optimization problems involving the rank or the l 0 -norm that is set cardinality are difficult to solve. Nevertheless, the DRMF problem admits a simple solution due to its decomposable structure in L and S. This problem can be solved by a block coordinate descent strategy. First, the current outliers S are fixed. Secondly, they are excluded from A to get the clean data C. Then, Xiong et al. [341] fit L based on C. Then, the outliers S are updated based on the current errors E = A -L. The algorithm solved the factorization problem as follows:

arg min L L,S ||C -L|| F subj rank(L) ≤ r (156)
where C = A -S.Then, the outlier detection problem is solved as follows:

arg min S L,S ||E -L|| F subj ||S|| l0 ≤ p (157)
where E = X -L. The solution to the low-rank approximation problem ( 156) is directly given by SVD. Since only the first r singular vectors are required, the computation is accelerated using a partial SVD algorithm.

7.3 LRM with Direct Robust Matrix Factorization for Rows (DRMF-R) Xiong et al. [341] proposed an extension of DRMF to deal with the presence of outliers in entire columns. This method is called DRMF-Row (DRMF-R). Instead of counting the number of outlier entries, the number of of outliers patterns is counted using the structured l 2,0 -norm. The direct formulation of DRMF-R is written as follows:

min

L,S ||A -S -L|| F subj rank(L) ≤ r, ||S|| l2,0 ≤ p ( 158
)
where p is the maximal number of outlier rows allowed. OP [342] is the convex relaxation of DRMF-R. No experimental results in Xiong et al. [341] are provided for DRMF-R on background/foreground separation.

7.4 Probabilistic Robust Matrix Factorization (PRMF) Wang et al. [321] proposed a probabilistic method for robust matrix factorization (PRMF) based on the l 1 -norm loss and l 2 -regularizer, which bear duality with the Laplace error and Gaussian prior, respectively. For model learning, Wang et al.

[321] used an efficient expectation-maximization (EM) algorithm by exploiting a hierarchical representation of the Laplace distribution as a scaled mixture of Gaussians. So, Wang et al. [321] considered the following probabilistic model:

A = U V ′ + S subj U ij |λ U ∼ N (U ij |0, λ -1 U ), V ij |λ V ∼ N (V ij |0, λ -1 V ) (159)
where U B -1 BV ′ = U V ′ holds for any r × r non singular matrix B. By exploiting U and V as model parameters with λ U , λ V and λ as hyperparameters with fixed values, MAP estimation is used to find U and V . From the rule of Bayes, the following equivalence can be written:

p(U, V |A, λ, λ U , λ V ) ∝ p(A|U, V, λ)p(U, λ U )p(V, λ V ) (160) Thus, log (p(U, V |A, λ, λ U , λ V )) = -λ||A -U V ′ || l1 - λ U 2 ||U || 2 l2 - λ V 2 ||V || 2 l2 + C (161)
where C is a constant term independent of U and V . Maximizing log (p(U, V |A, λ, λ U , λ V )) w.r.t. U and V is equivalent to the following minimization problem:

min U,V ||A -U V ′ || l1 + λ ′ U 2 ||U || 2 l2 + λ ′ V 2 ||V || 2 l2 ( 162
)
where λ ′ U = λU λ and λ ′ V = λV λ . Experimental results [321] on the I2R dataset [177] show that PRMF gives similar visual results than PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF], GoDec [389] and BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF] with less computation time. The corresponding code for PRMF 35 are provided in batch mode and online mode.

Algorithm for solving PRMF: While the model formulation given in Equation 161is rather straightforward, solving the optimization problem directly would be computationally challenging due to the non-smooth nature of the Laplace distribution. To address this computational issue, Wang et al. [321] reformulated the model by exploiting a two-level hierarchical representation of the Laplace distribution and EM algorithm is then used to solve this hierarchical model.

Bayesian Robust Matrix Factorization (BRMF)

PRMF present the following drawbacks: 1) it assumed that the basis and coefficient matrices are generated from zero-mean fixed-variance Gaussian distributions. This assumption is too restrictive, limiting the model flexibility needed for many realworld applications, 2) PRMF treats each pixel independently with no clustering effect but the moving objects in the foreground usually form groups with high within-group spatial or temporal proximity, 3) the loss function is defined based on the l 1 norm and it results to be not robust enough when the number of outliers is large. To address these three limitations, Wang et al. [366] proposed a full Bayesian formulation called Bayesian Robust Matrix Factorization (BRMF). The framework is similar to PRMF but the differences are the following ones:

1. For the generative proces, BRMF assumes that the mean vectors and precision matrices of the rows of U and V have conjugate priors. Learning the mean vectors offers more flexibility for the generation of A and learning the precision matrices captures the correlation between different features.

2. BRMF used a Laplace mixture with the generalized inverse Gaussian distribution as the noise model to further enhance model robustness. 3. BRMF contained a Markov extension (MBRMF) which assumes that the outliers exhibit spatial or temporal proximity.

Experimental results [366] on the SABS dataset [START_REF] Brutzer | Evaluation of background subtraction techniques for video surveillance[END_REF] show that MBRMF outperforms PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF], BRPCA [START_REF] Ding | Bayesian robust principal component analysis[END_REF], VBRPCA [17], DECOLOR [391] and PRMF [321]. The corresponding codes for BRMF 36 and MBRMF 36 are provided.

Practical Low-Rank Matrix Factorization (PLRMF)

Several LRM methods usually fail to minimize the l 1 -based nonconvex objective function sufficiently. Zheng et al. [386] proposed to add a convex nuclear-norm regularization term to improve convergence, without introducing too much heterogenous information. Thus, the robust l 1 -norm is choosen as the measurement. Then, Zheng et al. [386] enforced U to be column orthogonal to shrink the solution space, and added a nuclear-norm regularization term so as to improve convergence. This method called Practical Low-Rank Matrix Factorization (PLRMF) is based on the following problem: min

U,V ||W 5 ⊙ (A -U V )|| l1 (163)
where Ais the data matrix with a priori rank of r. The unknowns are U ∈ R m×r and V ∈ R r×n , which account for the rank constraint implicitly. W 5 is an indicator matrix of the same size as A where the entry value of 1 means that the component at the same position in A is observed, and 0 otherwise. The operator ⊙ denotes Hadamard element-wise matrix multiplication. U is constrained to be an orthogonal matrix and the following minimization problem is then obtained: min

U,V ||W 5 ⊙ (A -U V )|| l1 subj U T U = Ir (
U,V ||W 5 ⊙ (A -U V )|| l1 + λ||V || * subj U T U = Ir ( 165
)
where λ is is a weighting factor and should be small enough to keep the regularized objective in Equation 165. Experimental results [386] are presented on Structure from Motion (SFM) but not on background/foreground separation.

Algorithm for solving PLRMF (RegL1-ALM): Zheng et al. [386] developed a scalable first-order optimization algorithm to solve the regularized formulation on the basis of the augmented Lagrange multiplier (ALM) method with Gauss-Seidel iteration. The corresponding code called RegL1-ALM 37 is provided.

Low Rank Matrix Factorization with MoG noise (LRMF-MOG)

The previous low-rank factorization used loss functions such as the l 2 -norm and l 1norm losses. l 2 -norm is optimal for Gaussian noise, while l 1 -norm is for Laplacian distributed noise. However, real data in video are often corrupted by an unknown noise distribution, which is unlikely to be purely Gaussian or Laplacian. To address this problem, Meng et al. [216] proposed a low-rank matrix factorization problem with a Mixture of Gaussians (MoG) noise model. Since the MoG model is a universal estimator for any continuous distribution, it is able to represent a wider range of noise distributions. The parameters of the MoG model are estimated with a maximum likelihood method, while the subspace is computed with standard approaches. Thus, the decomposition is made as follows:

A = U V T + S.
Given the likelihood related to the MoG distribution, the aim is to maximize the loglikelihood function w.r.t the MoG parameters (Π,Σ) and the LRMF parameters (U, V ) as follows:

max U,V,Π,Σ i,j∈Ω K k=1 π k N(x ij |(u i ) T v j , σ 2 k ) (166)
with the MoG distribution constraints on S.

Experimental results [216] on the I2R dataset [177] show that LRMF-MOG outperforms SVD [START_REF] Golub | Matrix computation[END_REF], RSL [310], PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF], CWM [168], PCAl 1 [168].

Unifying Nuclear Norm and Bilinear Factorization

Cabral et al. [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF] proposed a low-rank matrix decomposition which can be achieved with both bilinear factorization and nuclear norm regularization models. By analyzing the conditions under which these two decompositions are equivalent, Cabral et al. [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF] proposed a unified model that inherits the benefits of both which is formulated as follows:

min U,V f (A -U V T ) + λ 2 (||U|| 2 F + λ||V || 2 F ) (167)
where L = U V T is the low-rank matrix with known rank r. Cabral et al. [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF] showed that the existence of local minima in Equation 167 depends only on the dimension r imposed on matrices U and V . Equation 167 can be written as follows:

min L,U,V ||W 5 ⊙ (A -L)|| l1 + λ 2 (||U|| 2 F + λ||V || 2 F ) subj L = U V T (168)
where W 5 ∈ R m×n is a weight matrix that is used to denote missing data (i.e., w 5,(ij) = 0) and observed data (i.e., w 5,(ij) = 1). Furthermore, Cabral et al. [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF] developed an ALM (Augmented Lagrange Multiplier) to solve Equation 167.

Experimental results [START_REF] Cabral | Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition[END_REF] on the I2R dataset [177] show that UNNBF outperforms PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF], GRASTA [119] and PRMF [321] with less computation times. show that LOIRE outperforms the PCP solved via IALM [START_REF] Candes | Robust principal component analysis? International[END_REF] and SemiSoft GoDec [389] in terms of computation efficiency with a similar visual accuracy.

Experimental Evaluation

First, we remind the challenges met in video-surveillance, a brief description of the Background Models Challenges (BMC) dataset [315] and the measures used for the performance evaluation. Then, we provide the evaluation and the comparison on 32 algorithms of decomposition into low-rank plus additive matrices. We have chosen the BMC dataset as a recent large-scale dataset in order to compare these algorithms. Qualitative and quantitative results are provided and discussed among this dataset. Several codes are available in the LRSLibrary [293] and some software packages for partial SVD computation can be found in [187].

Challenges in Video Surveillance

Three main conditions assures a good functioning of the background subtraction in video surveillance: the camera is fixed, the illumination is constant and the background is static. In practice, several challenges appear and perturb this process. They are the following ones:

-Noise image: It is due to a poor quality image source such as images acquired by a web cam or images after compression. -Camera jitter: In some conditions, the wind may cause the camera to sway back and so it cause nominal motion in the sequence. Foreground mask show false detections due to the motion without a robust maintenance mechanism. -Camera automatic adjustments: Many modern cameras have auto focus, automatic gain control, automatic white balance and auto brightness control. These adjustments modify the dynamic in the color levels between different frames in the sequence. -Illumination changes: They can be gradual such as ones in a day in an outdoor scene or sudden such as a light switch in an indoor scene. -Bootstrapping: During the training period, the background is not available in some environments. Then, it is impossible to compute a representative background image. -Camouflage: A foreground objects pixel characteristics may be subsumed by the modeled background. Then, the foreground and the background can be distinguished.

-Foreground aperture: When a moved object has uniform colored regions, changes inside these regions may not be detected. Thus, the entire object may not appear as foreground. Foreground masks contain false negative detections. -Moved background objects: Background objects can be moved. These objects should not be considered part of the foreground. Generally, both the initial and the new position of the object are detected without a robust maintenance mechanism.

-Inserted background objects: A new background object can be inserted.

These objects should not be considered part of the foreground. Generally, the inserted background object is detected without a robust maintenance mechanism.

-Dynamic backgrounds: Backgrounds can vacillate and this requires models which can represent disjoint sets of pixel values. -Beginning moving object: When an object initially in the background moves, both it and the newly revealed parts of the background called "ghost" are detected.

-Sleeping foreground object: Foreground object that becomes motionless cannot be distinguished from a background object and then it will be incorporated in the background. How to manage this situation depends on the context. Indeed, in some applications, motionless foreground objects must be incorporated and in others it is not the case. -Shadows: Shadows can be detected as foreground and can come from background objects or moving objects [125].

Background Models Challenge Dataset

The BMC (Background Models Challenge 38 ) dataset consists of both synthetic and real videos to permit a rigorous comparison of background subtraction techniques for the corresponding workshop organized within Asian Conference in Computer Vision (ACCV). This dataset [315] consists of the following sequences:

-Synthetic Sequences: A set of 20 urban video sequences rendered with the SiVIC simulator. With this tool, the associate ground truth was rendered frame by frame for each video at 25 fps. Several complex scenarios are simulated such as fog, sun and acquisition noise for two environments (a rotary and a street).

A first part of 10 synthetic videos are devoted to the learning phase, while the 10 others are used for the evaluation.

-Real Sequences: The BMC dataset also contains 9 real videos acquired from static cameras in video-surveillance contexts for evaluation. This real dataset has been built in order test the algorithms reliability during time and in difficult situations such as outdoor scenes. So, real long videos about one hour and up to four hours are available, and they may present long time change in luminosity with small density of objects in time compared to previous synthetic ones. Moreover, this real dataset allows to test the influence of some difficulties encountered during the object extraction phase, as the presence of vegetation, cast shadows or sudden light changes in the scene.

Performance Evaluation Metrics

We used ground truth based metrics computed from the true positives (TP), true negatives (TN), false positives (FP) and false negatives (FN). FP and FN refer to pixels misclassified as foreground (FP) or background (FN) while TP and TN account for accurately classified pixels respectively as foreground and background. Then, we computed the metrics used in the BMC dataset [279] such as the detection rate, the specificity, the false positive rate, the false negative rate, the percentage of wrong classifications, the precision and the F-Measure. Detection rate gives the percentage of corrected pixels classified as background when compared with the total number of background pixels in the ground truth:

DR = T P T P + F N ( 173
)
The specificity is computed as follows:

Specif icity = T N T N + F P (174)
The false positive rate and the false negative rate are defined as follows:

F P R = F P F P + T N (175) F N R = F N T P + F N (176)
The percentage of wrong classifications is defined as follows:

P W C = 100(F N + F P ) T P + F N + F P + T N (177)
Precision gives the percentage of corrected pixels classified as background as compared at the total pixels classified as background by the method:

P recision = T P T P + F P (178)

A good performance is obtained when the detection rate is high without altering the precision. A precision score of 1.0 means that every pixel labeled as belonging to the class foreground in the mask does indeed belong to the corresponding class in the ground-truth but pixels classified as background in the mask can be labeled incorrectly whereas a recall of 1.0 means that every pixel from the class foreground was labeled as belonging to the class foreground but says nothing about how many other pixels were incorrectly also labeled as belonging to foreground. We also computed the F-Measure (or effectiveness measure) as follows:

F = 2 * DR * P recision DR + P recision (179)
The F-Measure characterizes the performance of classification in Precision-Detection Rate space. The aim is to maximize F closed to one.

Experimental Results

We made the experimental evaluation by using the quality metrics which are computable thanks to a free software named BMC Wizard. The results of the first workshop BMC 2012 are available at the related website. We evaluated the performance of the following 32 algorithms grouped by category:

-Basic (2): PCA [233] and RSL [310].

-Robust Principal Components Analysis (20): 20 show the evaluation results using the synthetic videos. Table 21 and Table 22 shows the evaluation results using the real videos for evaluation phase. First, we provide a short qualitative analysis in presence of illumination changes and dynamic backgrounds. Then, we give a full quantitative evaluation. For the experimental setup, we used the parameters set in each original paper of the corresponding algorithm. To reduce the computation time in the experiments, we initialized batch algorihms with 200 frames instead of 300 frames (taken in the survey restricted on the RPCA framework [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF]). Thus, the performance in terms of F-measure are lower for the shared algorithms between this paper and the survey on the PCA framework [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF]. But, the rank of the algorithms is the same showing the stability of the algorithms and the evaluation. Thus, the conclusions made in [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF] are preserved. Because visual results of each algorithm seem very similar (as in [START_REF] Bouwmans | Robust PCA via principal component pursuit: A review for a comparative evaluation in video surveillance[END_REF]), we have done a full quantitative evaluation.

Analysis of the Experimental Results

Synthetic Videos

From Table 19 and Table 20, we can see that the algorithms that gives the best recall (more than 0.9) are the following ones (in italic in the tables): FAM [259], MahNMF [START_REF] Guan | MahNMF: manhattan non-negative matrix factorization[END_REF], ROSL [287] and DECOLOR [391]. For the precision, the algorithms R2PCP [127], OR-PCA with MRF [147], GRASTA [119], pROST [115] and DRMF present the higher precision. For the F-Measure, OR-PCA with MRF [147] gives the best performance followed by Lag-SCP-QN [START_REF] Aravkin | A variational approach to stable principal component pursuit[END_REF], OR-PCA without MRF [144], DECOLOR [391] and R2PCP [127]. We have indicated in the tables for the fifth best score the rank of the algorithm between parenthesis. Figure 4 shows a visual overview of the F-Measure for the 32 algorithms and the adaptive MOG [284]. The algorithm OR-PCA with MRF gives the highest F-Measure. For the stable version of PCA, Lag-SCP-QN [START_REF] Aravkin | A variational approach to stable principal component pursuit[END_REF] shows the best performance. DECOLOR [391] shows the best robustness for the low-rank methods but it is very time consuming instead of DRMF [341] and PRMF [321].

Real Videos

From Table 21 and Table 22, we can see that for the F-Measure, OR-PCA with MRF [147] gives the best performance followed by OR-PCA without MRF [144], ROSL [287], PRMF [321] and DRMF [341]. Thus, the algorithm OR-PCA with and without MRF has stable performance both on synthetic and real videos by preserving their rank. Otherwise, ROSL, PRMF and DRMF are more robust on real video than on synthetic videos. We has not conducted experiments on real videes for DECOLOR and RegL1-ALM due to very expensive computation time, and for GOSUS due to an implementation problem in the original code when very long sequences appear. Figure 5 shows a visual overview of the F-Measure for the 32 algorithms and the adaptive MOG [284].

Another main conclusion is that most of the RPCA algorithms outperforms the adaptive MOG for the synthetic videos as for real videos too. Figure 4 and Figure 5 show the F-Measure of the evaluated algorithms for synthetic and real videos, respectively. and LADM. The SVDs and CPU time of each algorithm were computed for each sequence. Table 23, Table 24, Table 25 and Table 26 group the average times. The CPU times are reported in the form hh:mm:ss for images of size 144×176 and with 200 frames for the training to allow easy comparison with other RPCA algorithms as the previous publications in this field present these performances on the I2R dataset [177] in this data format. In this paper, the results for the NSA [12] and the PSPG [START_REF] Aybat | Efficient algorithms for robust and stable principal component pursuit[END_REF] come from their authors. We can see that EALM and ADM are very computational expensive due to the fact that theses algorithms compute full SVDs. On these problems of extremely low ranks, the partial SVD technique used in IALM, LADMAP, LSADM and LADM become quite effective and reduce significantly the computation time. For the SPCP, the PSPG solver is the most efficient follows by the NSA and the ASALM. The variational BRPCA is less computational expensive than the BRPCA. The GoDec algorithm is the one which requires less time computation time in the RPCA category, and then it makes large videos applications reachable in realtime. ROSL and ROSL+ reach time requirement under 15 seconds as OR-PCA algorithms. In the RST category, GOSUS achieved less computation time than GRASTA and pROST. DRMF is the best in the category PRMF. Practically, we observed than the computation time of several algorithms increased differently in function of the image size because of their difference in term of complexity and memory requirement. Thus, the rank in term of computation time is not strictly the same at other image sizes, and several algorithms are only suitable in the case of low resolution. Finally, a full visual ranking of both matrix and tensor algorithms is provided at the LRSLibrary website.

Conclusion

In this paper, we have firstly presented a full review of recent advances on problem formulations based on the decomposition into low-rank plus additive matrices which are robust principal component analysis, robust matrix factorization, robust matrix completion, robust subspace recovery, robust subspace tracking and robust low-rank minimization. Thus, we proposed a unified view of the decomposition into low-rank plus additive matrices that we called DLAM. We In conclusion, this review for a comparative evaluation of robust subspace learning via decomposition into low-rank plus additive algorithms highlights the following points:

-Decomposition into low-rank plus additive matrices offers a suitable framework for background/foreground separation with robustness against illumination changes and dynamic backgrounds if spatial and temporal constraints are taken into account in the optimization problem with structured norms [275] or MRF [147]. Furthermore, several DLAM models can outperform state-of-theart models in background/foreground separation such as the MOG.

-The main disadvantages of the DLAM models is that their original version used batch algorithms that oftenly need too expensive computation time to reach real-time requirements. Thus, many effort have been done to reach realtime performance and to develop incremental algorithms as in the works of Rodriguez et al. [259][260][258][262],and Vaswani et al. [242][249][101].

-As images are stored in vectors which are oftenly exploited as is, DLAM models in their original version loss the spatial and temporal constraints. Thus, it is more appropriate 1) to add the use of Markov Random Fields [147],

2) to use structured norms aiming to preserve the spatial structures of images while being insensitive to outliers and missing data [343][275][302][145], or 3) to formulate RPCA in the two-dimensional case rather than via image to vector conversion, which enables the preservation of the image spatial information with reduced computational time. In the basic formulation, the resulting two-dimensional subspace methods mainly include two-dimensional PCA (2dPCA) [353], two-dimensional SVD (2dSVD) [START_REF] Ding | 2-dimensional singular value decomposition for 2D maps and images[END_REF], two-dimensional LDA (2dLDA) [364], two-directional two-dimensional PCA ((2d) 2 PCA) [371], Generalized Low Rank Approximations of Matrices (GLRAM) [362][363]. Robust formulations can be found for robust two dimensions (R2DRPCA) [301][302] and for Robust GLRAM (RGLRAM) [283].

Future research may concern less computational SVD algorithms such as LMSVD [195] for batch algorithms, and DLAM models which would be both incremental and real-time to reach the performance of the state-of-the-art algorithms [269][285] in term of computation time and memory requirements. Finally, DLAM models show a suitable potential for background modeling and foreground detection in video surveillance [152][153]. Furthermore, DLAM can been extended to the measurement domain, rather than the pixel domain, for use in conjunction with compressive sensing. Moreover, other research may concern the extension of DLAM in tensor-wise way to exploit fully spatial and temporal constraints [304][305][313][178][182][379]. The interest of the tensor approch over the matrix approach is investigated by Anandkumar et al. [START_REF] Anandkumar | Tensor vs matrix methods: Robust tensor decomposition under block sparse perturbations[END_REF]. Furthermore, efficient incremental algorithms have been recently developed for background/foreground separation [289][250][142][294][224].

11 Acknowledgment

Oh

  || * + µ||ω -1 (S)||2 A = L + S Yes Wang and Feng [322] Modified PCP (Implicit Regularizers) A = L + S min L,S λ||L|| * + ϕ(S) A = L + S He et al. [126] Yes

Fig. 3

 3 Fig. 3 Unified View of the Different Decomposition into Low-rank plus Additive Matrices (DLAM).

W

  ≥0,H≥0 ||A -W H|| l1 [163], the robust Low-rank (local-search) which minimizes min W,H ||A -W H|| l1 [163], XRay-l 2 [164]

6. 2

 2 Transformed Grassmannian Subspace Tracking (t-GRASTA)He et al. [122][123] proposed t-GRASTA (transformed-GRASTA) which iteratively performs incremental gradient descent constrained to the Grassmannian manifold of subspaces in order to simultaneously estimate a decomposition of a collection of images into a low-rank subspace, a sparse part of occlusions and foreground objects, and a transformation such as rotation or translation of the image. Based on RASL (Robust Alignment by Sparse and Low-rank decomposition)[239] which poses the robust image alignment problem as a transformed version of RPCA,He et al. [122][123] adapted Equation 140 as follows: min U,w,S,τ

  164) Therefore, Zheng et al. [386] used a nuclear-norm regularizer ||B|| * = ||U V || * where ||U V || * = ||V || * = due to U T U = Ir. The regularized minimization problem is then formulated as follows: min

Fig. 4 F

 4 Fig. 4 F-Measure of the 32 Algorithms for the Evaluation Results using the Synthetic Videos for Evaluation Phase.

Fig. 5 F

 5 Fig. 5 F-Measure of the 32 Algorithms for the Evaluation Results using the Real Videos for Evaluation Phase.

Table 2

 2 Robust problem formulations based on decomposition into low-rank plus additive matrices: A complete overview (Part 2).

	Methods	Decomposition	Minimization	Constraints	Convexity
	PCP

Table 3

 3 Decompositions in low-rank plus additive matrices: An homogeneous overview (Part 1)

	Categories	Methods	Decomposition	Minimization	Constraints	Convexity
	RPCA-PCP	SPCP	A = L + S + E	min L,S	||L|| * + λ||S||l 1	||A -L -S||F < δ	Yes
		Zhou et al. [393]
		Modified SPCP (Bilateral Projection)

Table 4

 4 Decompositions in low-rank plus additive matrices: An homogeneous overview (Part 2)

	Categories	Methods	Decomposition	Minimization	Constraints	Convexity
	Bayesian RPCA	BRPCA Ding et al. [71]	A = D(ZG)W2 + BX + E	-log (p(Θ|A, H))	Distribution constraints	-
		VBRPCA	A = DB T + S + E	p(A, D, B, S)	Distribution constraints	-
		Babacan et al. [17]
		MoG-RPCA	A = L + S	min KL divergence	MOG distribution constraints for S	-
		Zhao et al.[385]
	Approximated RPCA	GoDec	A = L + S + E	min L,S	||A -L -S|| 2 F	rank(L) ≤ e, card(S) ≤ k	No
		Zhou and Tao [389]
		Semi-Soft GoDec	A = L + S + E	min L,S	||A -L -S|| 2 F	rank(L) ≤ e, card(S) ≤ τ	No
		Zhou and Tao [389]					τ is a soft threshold
	Sparse Additive Matrix Factorization	SAMF	A =	K k=0 S + E	r(Θ)		Distribution constraints
		Nakajima et al. [222]
	Variational Bayesian Sparse Estimator	VBSE

  mixed norm of the matrix M , and it corresponds to the l β -norm of the vector formed by taking the l α -norms of the columns of the underlying matrix. α and β are in the interval [0, 2]. For example, ||M || l2,0 corresponds to the number of non-zero columns of the matrix M[375]. ||M || l2,1 forces spatial homogeneous fitting in the matrix M[111], and it is suitable in presence of column outliers or noise[306][111][124]. ||M || l2,1 is equal to j ||M :j || l 2[375]. The influence of α and β on the matrices L and S was studied in[111]. This norm is also called structured norm.

  ||.|| d is the d dual norm of any norm ||.||norm previously defined, that is norm ∈ lα, l∞, l α,β , Lα, F, * . For example, the dual norm of the nuclear norm is ||.|| 2 called spectral norm which corresponds to the largest singular value of the matrix [322].

	-Matrix nuclear norm: ||M || * is the nuclear norm of the matrix M , and it corresponds to the sum of its singular values [375]. The nuclear norm is
	the l 1 -norm applied on the vector composed with the singular values of the
	matrix [66]. It is equivalen to the Ky Fan n-norm and the Schatten-1-norm
	[280].
	-Matrix Schatten-α norm with 0 < α ≤ 2: The Schatten-α norm ||M || Sα = ( min(m,n) k=1

-Matrix dual norm:

  -Matrix max-norm: The max-norm ||M ||max is defined as m

	Functions/Constraints	Original loss function	Surrogate loss functions
	flow ()	rank(.) [40]	PCP/SPCP/QPCP/BPCP/LPCP/OP
	Low-rank L		Nuclear norm [40]
			Truncated nuclear norm [128]
			Modified PCP
			Capped nuclear norm [299]
			Schatten-α norm [320][276], Rank-N [231]
			Heuristic Recovery
			Log-sum heuristic (LHR) [66]
			Stochastic Optimization
			Max-norm (MRMD) [280]
	fsparse()	l0-norm [40]	PCP/Modified PCP
	Sparsity S		l1-norm [40], capped l1-norm [299]
			lα-norm [320], Lα-seminorm [276]
			dual norm [322], M-estimator [126]
			SPCP
			Generalized fused Lasso [340]
			BPCP/LPCP
			l2,1-norm [306][333]
			OP
			l1,2-norm [205]
			LHR
			Log-sum heuristic [66]
			pROST
			Lp-quasi-norm [268]
	fnoise()	PCP	Modified PCP
	Error E
			ij	ax |M ij |. The
	max-norm is equivalent to the l∞-norm [280].
	4. Loss functions: Several loss functions can be used to enforce the low-rank,
	sparse and noise constraints on L, S, and E respectively. These functions are
	respectively noted f low (), fsparse() and f noise (). Most of the time, loss func-
	tions are defined on the previous defined norms such as: l 0 -loss function (||.|| l0 ),
	l 1 -loss function (||.|| l1 ), l 2 -loss function (||.|| l2 ), nuclear norm function, Frobe-
	nius loss function and Log-sum heuristic function [66]. Other loss functions
	can be used such as lσ-loss function [360], Least Squares (LS) loss function
	(||.|| 2 F ), Huber loss function [10], M -estimator based loss functions [126] and
	the generalized fused Lasso loss function [340]. Lipschitz loss function can also
	be used in a two-stage convex relaxation approach [116] by the majorization for
	a class of locally Lipschitz continuous surrogates of Equation 21, which solves
	the nuclear norm plus l 1 -norm minimization problem in the first stage and a
	nuclear semi-norm plus weighted l 1 -norm minimization problem in the second
	stage with theoretical guarantee.

Table 7

 7 Loss functions f low (), fsparse() and f noise () used for the low-rank, sparse and noise constraints in the different problem formulations.

Table 7

 7

	shows an

  λ 1 and λ 2 are the regularization parameters. C 2 is a constraint on A, L and S. f low (L) is a loss function which constrains the matrix L to be low-rank as the following ones: rank(.), or the surrogated norm ||.|| * . fsparse(S) is a loss function which constrains the matrix S to be sparse as the following ones: ||.|| l0 , or the surrogated norm ||.|| l1 . An overview of the different loss functions f low (.) and fsparse() which are used in the literature are shown in Table

  p 1 and p 2 are taken in the set {1, 2}. ||.||norm 1 and ||.||norm 2 could be any norm of the following set of norms: lα-norm, l∞, l α,β mixed norm, Lα-seminorm, Frobenius norm, nuclear norm, dual norm and Schatten norm. ||.||norm 1 and ||.||norm 2 are taken to enforce the low-rank and sparsity constraints of L and S, respectively. The constraint C 1 is generally based on

  3 + r 2 (m + n) + r 2 m + rmn)

	Liu et al. [198]	Opre=Unknown, Oconv =Unknown
	PRoximal Iterative SMoothing Algorithm (PRISMA)	Oiter (nm), Opre(log(ǫ)/ǫ), Oconv=Unknown
	Orabona et al. [235]
	Fast Alterning Minimization (FAM) 7	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Rodriguez and Wohlberg [259]
	Fast Alternating Direction Method of Multipliers (FADMM)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Yang and Wang [358]
	Fast Alternating Direction Method with Smoothing Technique (FADM-ST)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Yang [356]
	Online solvers
	Online Alternating Direction Method (OADM)
	Wang and Banerjee [318]	Oiter =Unknown, Opre=Unknown, Oconv (1/T )
	Non convex solvers
	Difference of Convex (DC)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Sun et al. [299]
	Fast Alternating Difference of Convex (FADC)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Sun et al. [299]
	Non-convex Alternating Projections(AltProj)	Oiter (r 2 mn), Opre(log(1/ǫ)), Oconv =Unknown
	Netrapalli et al. [227]
	2D solvers
	Iterative method for Bi-directional Decomposition (IMBD)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Sun et al. [301]

Table 10

 10 Solvers for RPCA-PCP: An overview of their complexity per iteration at running time O iter , their complexity Opre to reach an accuracy of ǫ precision and their convergence rate Oconv for T iterations. "Unknown" stands for not indicated by the authors.

	Solvers	Complexity
	RPCA via SPCP (RPCA-SPCP)	Zhou et al. [393]
	Alternating Splitting Augmented Lagrangian method (ASALM 8 )	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Tao and Yuan [307]
	Variational ASALM (VASALM 8 )	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Tao and Yuan [307]
	Parallel ASALM (PSALM 8 )	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	Tao and Yuan [307]
	Non Smooth Augmented Lagrangian Algorithm (NSA 9 )

Table 12

 12 Oiter (|Ω| d 3 + d |Ω| + nd 2 ), Opre=unknown, Oconv=unknown with the Grassmannian geodesic gradient descent He et al. [119][120] Boyd et al. [34], Edelman et al. [78] pROST 22

	Methods	Solvers	Complexity
	Robust Subspace Tracking (RST)
	GRASTA 21	Augmented Lagrangian algorithm
		Conjugate Gradient	Oiter =unknown, Opre=unknown, Oconv =unknown
	Hage and Kleinstauber [115]	Hage and Kleinstauber [115]
	GOSUS 23	ADMM	Oiter =unknown, Opre=unknown, Oconv =unknown
	Xu et al. [343]	Xu et al. [343]
	FARST Ahn [4]	ADMM Ahn [4]	Oiter (|Ω| d 3 + d |Ω| + nd 2 ), Opre=unknown, Oconv=unknown
	Robust Low Rank Minimization (RLRM)
	LRM with Contiguous Outliers Detection	Alternating Algorithm (SOFT-IMPUTE)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	(DECOLOR 24 ) Zhou et al. [391]	Mazumder et al. [215]
	LRM with DRMF	Block Coordinate Descent strategy	Oiter (mn(r + log(p)))
	(DRMF 25 ) Xiong et al. [341]	Xiong et al. [341]	K partial SVD at each iteration
	LRM with DRMF-R	Block Coordinate Descent strategy	Oiter (mn(r + log(p)))
	(DRMF-R 25 ) Xiong et al. [341]	Xiong et al. [341]	K partial SVD at each iteration
	PRMF	Conditional EM Algorithm (CEM)	Oiter =Unknown, Opre=Unknown, Oconv =Unknown
	(PRMF

Solvers for RNMF, RMC, and RSR: An overview of their complexity per iteration at running time O iter , their complexity Opre to reach an accuracy of ǫ precision and their convergence rate Oconv for T iterations. "Unknown" stands for not indicated by the authors. 26 ) Wang et al. [321] Jebara and Pentland [148] BRMF Conditional EM Algorithm (CEM) Oiter =Unknown, Opre=Unknown, Oconv =Unknown (BRMF 27 ) Wang and Yeung [366] Jebara and Pentland [148] PLRMF Inexact ALM with Gauss-Seidel iteration Oiter = r(max(m, n)) 2 , Opre=Unknown, Oconv=Unknown (RegL1-ALM 28 ) Zheng et al. [386] Zheng et al. [386] LRMF-MOG EM algorithm Oiter =Unknown, Opre=Unknown, Oconv =Unknown Meng et al. [216] Meng et al. [216] UNN-BF ALM Oiter (mnr + nr 2 ), Opre=Unknown, Oconv =Unknown Cabral et al.[37] Cabral et al.[37] RRF Alternative Direction Descent Algorithm (ADDA) Oiter =Unknown, Opre=Unknown, Oconv =Unknown Sheng et al. [282] Sheng et al. [282]

3 .

 3 Incremental algorithms: Incremental algorithms are needed to update the low-rank and additive matrices when a new data arrives. Several incremental algorithms can be found in the literature as follows: 1) in the RPCA framework(PCP [242][244]

[243][245][100][331][261][345]

, SPCP

[295]

,

[214], Approximated RPCA

[START_REF] Chen | Incremental low-rank and sparse decomposition for compressing videos captured by fixed cameras[END_REF]

), and 2) in the subspace tracking framework

(GRASTA [119]

[120], t-GRASTA

[122][123]

, pROST

[115][268]

, GOSUS

[343]

and FARST

  Waters et al. [329][330] developed an algorithm for solvingEquation 19, for sparse vector recovery and ADMiRA[173] for low-rank matrix recovery. To accelerate the convergence speed of SpaRCS, Kyrillidis and Cevher proposed an algorithm called Matrix ALPS and based on acceleration techniques from convex analysis and exploited well-known memory-based acceleration technique. As incorporating priori knowledge into the basic compressive sensing results in significant improvement of its performance,Zoonobi and Kassim [395] extended SpaRCS with partial known support.Jiang et al. [152] reformulated the problem 19) into an equivalent problem by introducing some splitting variables, and applied the ADM framework. Furthermore,

  Rachford splitting method. The convex optimization problem is solved by canceling the constraint of the variables, and the proximity operators of the objective function are computed alternately. The new algorithm can exactly recover the low-rank and sparse simultaneously, and it is proved to be convergent. Another approach developed byChai et al.

proposed to alleviate this ignorance by the Alternating Direction Method of Multipliers (ADMM) which minimizes the variables L and S serially. The convergence of ADMM for convex objective functions has been proved

[229][300][317][183][186][185]

. The iteration complexity was analyzed for Multi-Block ADMM

[184]

. ADMM achieves it with less computation cost than ALM. Recently, Chartrand

[START_REF] Chartrand | Non convex splitting for regularized low-rank and sparse decomposition[END_REF]

proposed a non convex splitting version of the ADMM

[368]

called NCSADM. This non convex generalization of [188] produces a sparser model that is better able to separate moving objects and stationary objects. Furthermore, this splitting algorithm maintains the background model while removing substantial noise, more so than the convex regularization does. The ALM neglects the separable structure in both the objective function and the constraint. Thus,

Zhang and Liu [376]

proposed a variant of the Douglas-Rachford splitting method (VDRSM) for accomplishing recovery in the case of illumination changes and dynamic backgrounds. In a similar way,

Zhu et al. [394]

proposed a Proximity Point Algorithm (PPA) based on the Douglas-

  . In another way,Orabona et al. [235] proposed an optimization algorithm called PRoximal Iterative SMoothing Algorithm (PRISMA) which is decomposed into three parts: a smooth part, a simple non-smooth Lipschitz part, and a simple non-smooth non-Lipschitz part. Furthermore, a time variant smoothing strategy is used to obtain a guarantee that does not depend on knowing in advance the total number of iterations nor a bound on the domain. Numerical experiments[235] show that PRISMA required less iterations than Fast-LSADM[START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex function[END_REF]. Another approach developed byRodriguez and Wohlberg [259] is able to compute a sparse approximation even after the first outer loop, (taking approximately 12 seconds for a 640 × 480 × 400 color test video) which is approximately an order of magnitude faster than IALM[188] with the same accuracy.Yang and Wang [358] proposed a Fast Alternating Direction Method of Multipliers (FADMM) algorithms which outperforms slightly IALM[188] and ADM[368] in term of computation times.Yang [356] improved FADMM by using a smoothing technique which is used to smooth the non-smooth terms in the objective function. The previous solvers are mainly bacth ones but online algorithms are better adapted for real-time application. So, proposed an efficient online learning algorithm named online ADM (OADM) which can solve online convex optimization under linear constraints where the objective could be nonsmooth.

-Online solvers:

-Non convex solvers:

Sun et al. [299]

developed for the non-convex RPCA formulation of RPCA with capped norms two algorithms called Difference of Convex (DC) and Fast Alternating Difference of Convex (FADC), respectively. DC programming treats a non-convex function as the difference of two convex functions, and then iteratively solves it on the basis of the combination of the first convex part and the linear approximation of the second convex part. Numerical measurements [299] demonstrate that DC approach performs better than both IALM [188] and NSA [12] in terms of the low-rank and sparsity. In another way, Netrapalli et al. [227] proposed a Non-convex Alternating Projections algorithm (AltProj) to solve a non-convex formulation of RPCA. The overall complexity of AltProj is O(r 2 mnlog(1/ǫ)). This is drastically lower than the best known bound of O(m 2 n/ǫ) on the number of iterations required by convex methods, and just a factor r away from the complexity of naive PCA. AltProj is around 19 times faster than IALM [188]. Moreover, visually, the background extraction seems to be of better accuracy.

1 .

 1 Recursive Projected Compressive Sensing (ReProCS): To address the two first issues,Qiu and Vaswani [242] proposed an online approach called Recursive Robust PCP (RR-PCP) in[242], and Recursive Projected Compressive Sensing (ReProCS) in[243][246][249][248]. The aim of ReProCS is to causally keep updating the sparse matrix S t at each time, and keep updating the principal directions. The t-th column of A, A t , is the data acquired at time t and can be decomposed as follows:

  the matrix U is an unknown m × m orthonormal matrix. The support of S t changes slowly over time. Let N t denote the support of x t which is assumed piecewise constant over time and given an initial estimate of P t = (U) Nt = Pt ,Qiu and Vaswani [242] solved for the sparse component S t by finding the orthogonal complement matrix Pt,⊥ , and then using the projection M t onto Pt,⊥ , denoted by y t :

y t = P T t,⊥ M t = P T t,⊥ L t + P T t,⊥ S t (25) to solve S t . The low-rank component is closed to zero if Pt ≈ P t , otherwise new directions are added. Furthermore, recent estimates of L t = A t -S t are stored and used to update P t . Confirming the first results obtained in [370][201], a correctness result for ReProCS is given by

5 .

 5 Practical ReProCS (Prac-ReProCS) :Guo et al. [100][102][103][101] designed a practically usable modification of the theoretical ReProCS algorithm. This practical ReProCS (Prac-ReProCS) requires much fewer parameters which can be set without any model knowledge and it exploits practically valid assumptions such as denseness for L t , slow subspace change for L t , and correlated support change of S t .6. Incremental Low-Rank (iLR) Algorithm: Wei et al.[331] proposed an incremental low-rank matrix decomposition algorithm that maintains a clean background matrix adaptive to dynamic changes with both effectiveness and efficiency guarantees. Instead of a batch RPCA which requires a large number of video frames (usually 200 frames) for each time period, 15 frames only are required with iLR. iLR algorithm is about 9 times faster than a batch RPCA.

7.

Incremental PCP (incPCP):

Rodriguez and Wohlberg [261][262]

proposed an incremental PCP which processes one frame at a time. Obtaining similar results to batch PCP algorithms, it has an extremely low memory footprint and a computational complexity that allows real-time processing. Furthermore, incPCP is also able to quickly adapt to changes in the background. A Matlabonly implementation of this algorithm [260] running on a standard laptop (Intel i7-2670QM quad-core, 6GB RAM, 2.2 GHz) can process color videos of size 640 and 1920 at a rate of 8 and 1.5 frames per second respectively. On the same hardware, an ANSI-C implementation [260] can deliver a rate of 49.6 and 7.2 frames per second for grayscale videos of size 640 and 1920 respectively. This algorithm has real-time performance on GPU [258]. Furthermore, Rodriguez and Wohlberg [263] developed a translational and rotational jitter invariant incPCP which reach real-time performanc on GPU [257]. 2.1.4 Methods for real time implementation of PCP

  norm relaxation,Sun et al. [301] developed an iterative method to solve Equation (37) efficiently via alternating optimization, by specific greedy algorithm for the l 0 -norm regularization. So, a robust 2D-PCA model by sparse regularization is then solved via an alternating optimization algorithms. Results on dynamic backgrounds from the I2R dataset[177] show the effectiveness of the Robust 2D-PCA (R2DPCA), compared with the conventional 2D-PCA [162] and PCP solved via IALM [188]. 9. Rank-N Soft Constraint:Oh [231] proposed a RPCA with Rank-N Soft Constraint (RNSC) based on the observation that the matrix A should be rank N without corruption and noise. So, the decomposition is formulated as estimating sparse error matrix and minimizing rank of low-rank matrix consisting of

N principal components associated to the N largest singular largest values. Thus, the miminization problem with rank-N soft constraint is the following one: min

Table 16

 16 Stable Principal Component Pursuit: A Complete Overview. The first column indicates the concerned category and the second column the name of each method. Their corresponding acronym is indicated in the first parenthesis. The third column gives the name of the authors and the date of the related publication.

  O(mr 2 + lr 2 ) time for r = max i k Ci , O(mr 2 + lr 2 ) time for the pseudoinversion of C i and O(mr 2 + lr 2 ) time for matrix multiplication with each C i in parallel. For the DFC-NYS, the per-iteration complexity O(mlr C ) where r C is the rank of the matrix C and O(mlr R ) where r R is the rank of the matrix R. The cost of combining requires O(lr 2 + dr 2 + min(m, n)r 2 ) time where r = max(r C , r R ). Mackey et al.

  2.4 RPCA via Block based Principal Component PursuitTang and Nehorai [306] proposed a block based PCP (BPCP) that enforces the lowrankness of one part and the block sparsity of the other part. This decomposition involves the same model than PCP in Equation(20), that is A = L + S, where L is the low-rank component but S is a block-sparse component. The low-rank matrix L and the block-sparsity matrix S can be recovered by the following optimization problem[389]: min

Table 17

 17 Bayesian Robust Principal Component Analysis: A Complete Overview. The first column indicates the concerned category and the second column the name of each method. Their corresponding acronym is indicated in the first parenthesis. The third column gives the name of the authors and the date of the related publication.

	real time applicability. This scheme was improved with dynamic feature selection
	[146][143], a depth extended version with spatiotemporal constraints [141], and for
	noisy videos [140].

  S. P (S; λ) = ij λ ij ||s ij || l 2 where ||.|| l 2 allows the minimization to address outliers in a row-wise manner. She et al. [277] used generalized M -estimator to solve this minimization. The computation is related to the orthogonality constraint, in addition to the non-smooth and possibly non-convex P , and She et al. [277] developed a fast alternating optimization algorithm on the basis of Stiefel manifold optimization and iterative nonlinear thresholdings.

	6 Robust Subspace Tracking

  7.9 Robust Rank Factorization (RRF)Sheng et al. [282] proposed a l 1 -regularized Outlier Isolation and REgression (LOIRE) model. The measurement process is written in the following decomposition form:A = L + S + E = BX + S + E (169)where A the observation through B, S denotes the outlier vector and E denotes a dense Gaussian noise. By adding a penalty term ||S|| 1 to the least mean squares on E, Sheng et al.[282] derived a minimization problem for X as follows:||E|| l2 + µ||S|| l1 subj µ > 0, A = BX + S + E ||S|| l1 subj µ > 0, A = BX + S + E(171)where the matrix B is generally unknown, then a simple way to find a most appropriate matrix B that fits the problem is to search one that minimizes the above optimization problem: ||S|| l1 subj λ > 0, A = BX + S +

				min S,E	(170)
	where µ is a regularization term. This formulation is then extended to realize
	robust rank factorization which can be applied to recover low-rank structures from
	massive contaminations.
	min S,X	||A -BX -S|| 2 l2 +
	min B	min S,X	||A -BX -S|| 2 l2 +

λ 2 λ 2 E (172) To ensure a unique solution for matrix B and X, each column of B should have a unit length. To solve LOIRE, Sheng et al. [282] developed an Alternative Matrix Descent Algorithm (AMDA). Experimental results [282] on the I2R dataset [177]

Table 19 and

 19 Table

Table 19

 19 BMC dataset: Evaluation Results using the Synthetic Videos for Evaluation Phase (Part 1). The number between parenthesis indicates the rank of the algorithm in terms of F-measure.9 Implementation and Computational CostAll the algorithms are implemented in Matlab and come from the LRSLibrary39 . The computational cost of the RPCA algorithms is mainly related to the singular value decomposition (SVD). It can be reduced significantly by using a partial SVD because only the first largest few singular values are needed. Practically, the implementation available in PROPACK 40 are used for the IALM, LADMAP, LSADM

	Algorithms	Measure			Street					Rotary			Average
			112	212	312	412	512	122	222	322	422	522
	Robust NMF
	MahNMF	Recall	0.982	0.982	0.934	0.909	0.943	0.979	0.977	0.944	0.894	0.969	-
	Guan et al. [99]	Precision	0.661	0.622	0.520	0.515	0.524	0.624	0.612	0.538	0.524	0.576	-
		F-measure	0.790	0.762	0.668	0.657	0.674	0.762	0.752	0.685	0.660	0.722	0.713
	Robust Subspace Recovery
	ROSL	Recall	0.984	0.985	0.948	0.924	0.948	0.984	0.983	0.950	0.912	0.976	-
	Shu et al. [287]	Precision	0.744	0.744	0.528	0.519	0.527	0.753	0.753	0.545	0.533	0.619	-
		F-measure	0.847	0.847	0.678	0.665	0.678	0.853	0.853	0.693	0.673	0.757	0.754
	Robust Subspace Tracking
	GRASTA	Recall	0.700	0.787	0.695	0.787	0.669	0.680	0.637	0.619	0.623	0.791	-
	He et al. [119]	Precision	0.980	0.847	0.965	0.843	0.960	0.902	0.548	0.530	0.778	0.714	-
		F-measure	0.817	0.816	0.807	0.814	0.789	0.776	0.589	0.571	0.692	0.751	0.618
	pROST		0.944	0.878	0.853	0.889	0.785	0.819	0.838	0.789	0.738	0.863	-
	Hage and Kleinsteuber [115]	Precision	0.844	0.937	0.968	0.931	0.961	0.903	0.847	0.953	0.730	0.815	-
		F-measure	0.891	0.906	0.907	0.909	0.864	0.859	0.842	0.863	0.734	0.838	0.718
	GOSUS	Recall	0.982	0.981	0.945	0.933	0.949	0.978	0.977	0.958	0.906	0.969	-
	Xu et al. [343]	Precision	0.617	0.609	0.526	0.526	0.528	0.616	0.612	0.555	0.534	0.577	-
		F-measure	0.758	0.751	0.676	0.672	0.678	0.756	0.752	0.702	0.672	0.723	0.714
	Low Rank Minimization
	DECOLOR	Recall	0.982	0.985	0.983	0.980	0.978	0.983	0.983	0.981	0.967	0.980	-
	Zhou et al. [391]	Precision	0.778	0.748	0.747	0.729	0.599	0.764	0.759	0.760	0.762	0.694	-
		F-measure	0.868	0.851	0.849	0.836	0.743	0.860	0.857	0.857	0.852	0.813	0.838 (4)
	DRMF	Recall	0.857	0.864	0.853	0.864	0.903	0.880	0.834	0.744	0.805	0.827	-
	Xiong et al. [341]	Precision	0.969	0.950	0.968	0.948	0.855	0.891	0.834	0.924	0.700	0.825	-
		F-measure	0.910	0.905	0.907	0.904	0.878	0.885	0.834	0.824	0.749	0.826	0.710
	PRMF	Recall	0.944	0.903	0.918	0.901	0.899	0.891	0.906	0.867	0.824	0.869	-
	Wang et al. [321]	Precision	0.819	0.919	0.879	0.922	0.862	0.887	0.849	0.888	0.707	0.821	-
		F-measure	0.877	0.911	0.898	0.911	0.880	0.889	0.877	0.878	0.761	0.845	0.727
	PLRMF (RegL1-ALM)	Recall	0.984	0.985	0.953	0.926	0.949	0.984	0.983	0.950	0.914	0.976	-
	Zheng et al. [386]	Precision	0.744	0.745	0.533	0.520	0.529	0.756	0.756	0.545	0.535	0.619	-
		F-measure	0.847	0.848	0.684	0.666	0.679	0.855	0.855	0.693	0.675	0.758	0.756
	Mixture of Gaussians
	Adaptive MOG	Recall	0.827	0.827	0.797	0.761	0.821	0.823	0.831	0.797	0.743	0.834	-
	Shimada et al. [284]	Precision	0.766	0.768	0.480	0.426	0.519	0.786	0.790	0.526	0.435	0.740	-
		F-measure	0.796	0.796	0.605	0.553	0.640	0.804	0.810	0.638	0.555	0.784	0.698

Table 20

 20 BMC dataset: Evaluation Results using the Synthetic Videos for Evaluation Phase (Part 2). The number between parenthesis indicates the rank of the algorithm in terms of F-measure.

Table 21

 21 BMC dataset: Evaluation Results using the Real Videos for Evaluation Phase (Part 1).evaluated its adequation to the application of background/foreground separation by investigated how these methods are solved and if incremental algorithms and real-time implementations can be achieved. Finally, experimental results on the Background Models Challenge (BMC) dataset show the comparative performance of these recent methods.

	Algorithms	Measure				Real Videos				Average
			001	002	003	004	005	006	007	008	009
	Robust NMF
	MahNMF	Recall	0.857	0.822	0.901	0.848	0.802	0.823	0.788	0.536	0.716	-
	Guan et al. [99]	Precision	0.646	0.671	0.739	0.529	0.520	0.620	0.672	0.505	0.599	-
		F-measure	0.737	0.739	0.812	0.651	0.631	0.708	0.726	0.520	0.652	0.686
	Robust Subspace Recovery
	ROSL	Recall	0.743	0.837	0.912	0.851	0.823	0.843	0.778	0.562	0.768	-
	Shu et al. [287]	Precision	0.865	0.731	0.779	0.531	0.512	0.680	0.684	0.508	0.852	-
		F-measure	0.799	0.781	0.840	0.654	0.631	0.753	0.728	0.534	0.808	0.725 (3)
	Robust Subspace Tracking
	GRASTA	Recall	0.719	0.767	0.852	0.823	0.533	0.802	0.751	0.673	0.730	-
	He et al. [119]	Precision	0.542	0.845	0.963	0.796	0.516	0.711	0.900	0.696	0.950	-
		F-measure	0.618	0.804	0.904	0.809	0.524	0.754	0.819	0.684	0.826	0.674
	pROST	Recall	0.824	0.672	0.923	0.835	0.760	0.797	0.596	0.741	0.850	-
	Hage and Kleinsteuber [115]	Precision	0.632	0.844	0.958	0.884	0.631	0.787	0.796	0.770	0.874	-
		F-measure	0.715	0.749	0.940	0.859	0.689	0.792	0.682	0.755	0.862	0.704
	Low Rank Minimization
	DRMF	Recall	0.828	0.719	0.934	0.874	0.772	0.823	0.617	0.762	0.842	-
	Xiong et al. [341]	Precision	0.600	0.856	0.949	0.877	0.618	0.801	0.811	0.766	0.929	-
		F-measure	0.696	0.782	0.941	0.875	0.686	0.812	0.701	0.764	0.883	0.714 (5)
	PRMF	Recall	0.830	0.720	0.937	0.874	0.790	0.805	0.596	0.746	0.839	-
	et al. [321]	Precision	0.665	0.853	0.949	0.880	0.627	0.790	0.824	0.779	0.962	-
		F-measure	0.738	0.781	0.943	0.877	0.699	0.797	0.691	0.762	0.897	0.718 (4)
	Mixture of Gaussians
	Adaptive MOG	Recall	0.849	0.580	0.859	0.829	0.754	0.780	0.691	0.723	0.828	-
	Shimada et al. [284]	Precision	0.682	0.546	0.780	0.580	0.435	0.636	0.603	0.495	0.790	-
		F-measure	0.757	0.562	0.818	0.785	0.558	0.702	0.644	0.591	0.809	0.680

Table 22

 22 BMC dataset: Evaluation Results using the Real Videos for Evaluation Phase (Part 2).

						RPCA-PCP
	EALM [188]	IALM [188]	ADM (LRSD) [368]	LADMAP [190]	LSADM [90]	LADM (LMaFit) [281]
	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU
	550	00:40:15	38	00:03:47	510	00:35:20	16	00:05:20	43	00:04:03	35	00:04:55

Table 23

 23 Time Performance of RPCA Algorithms: PCP.

			RPCA-SPCP					OR-PCA
	ASALM [307]	NSA [12]	PSPG [13]	OR-PCA [144]	OR-PCA with MRF [147]
	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU
	94	00:15:17	19	00:03:07	23	00:01:05	-	00:0:12	-	00:00:14

Table 24

 24 Time Performance of RPCA Algorithms: SPCP and OR-PCA.

		Bayesian RPCA		Approximated RPCA
	BRPCA [71]	VBRPCA [17]	GoDec [389]	Semi Soft GoDec [389]
	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU
	-	00:04:01	-	00:01:07	-	00:00:50	-	00:00:55

Table 25

 25 Time Performance of RPCA Algorithms: Bayesian RPCA and Approximated RPCA.

		RSR					RST				RLRM
	ROSL [287]	ROSL+ [287]	GRASTA [119]	pROST [115]	GOSUS [343]	DRMF [341]	PRMF [321]
	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU	SVDs	CPU
	-	00:00:15	-	00:00:06	-	00:01:06	-	00:01:05	-	00:00:23	-	00:00:17	-	00:00:36

Table 26

 26 Time Performance of RSR, RST and RLRM Algorithms.

http://perception.csl.uiuc.edu/matrix-rank/samplecode.html

http://www.caam.rice.edu/ yzhang/LMSVD/lmsvd.html

http://math.nju.edu.cn/ jfyang/LRSD/index.html

Available on request by email to the corresponding author

Available on request by email to the corresponding author

http://lmafit.blogs.rice.edu/

https://sites.google.com/a/istec.net/prodrig/Home/en/pubs

Available on request by email to the corresponding author

http://www2.ie.psu.edu/aybat/codes.html

http://www.uni-graz.at/imawww/ifb/r2pcp/index.html

https://sites.google.com/site/godecomposition/GreBsmo.zip

http://www2.ece.ohio-state.edu/ schniter/BiGAMP/BiGAMP.html

http://cvxr.com/tfocs/

http://guppy.mpe.nus.edu.sg/ mpexuh/publication.html

http://people.ee.duke.edu/ lcarin/BCS.html

http://www.dbabacan.info/software.html

http://sites.google.com/site/godecomposition/code

https://sites.google.com/site/nmfsolvers/

http://cs.tju.edu.cn/faculty/likun/projects/bf-separation/index.htm

https://sites.google.com/site/xianbiaoshu/

http://sites.google.com/site/hejunzz/grasta

http://www.gol.ei.tum.de/index.php?id=37 L=1

http://pages.cs.wisc.edu/ jiaxu/projects/gosus/

http://rsl-cv2015.univ-lr.fr/workshop/

https://sites.google.com/site/robustdlam/

http://www.vision.eecs.ucf.edu/projects/Turbulence/

http://bioinformatics.ust.hk/decolor/decolor.html

http://winsty.net/prmf.html

http://winsty.net/brmf.html

https://github.com/andrewssobral/lrslibrary

http://soi.stanford.edu/rmunk/PROPACK/

Solvers 1) Basic solvers

Singular Values Decomposition (SVT) Cai et al. (2008) [38] Iterative Thresholding (IT) Wright et al. (2009) [335] Accelerated Proximal Gradient (APG) Lin et al.(2009) [189] Dual Method (DM) Lin et al.(2009) [189] Exacted Augmented Lagrangian Method (EALM) Lin et al. (2009) [188] Inexact Augmented Lagrangian Method (IALM) Lin et al. (2009) [188] Alternating Direction Method (ADM) Yuan and Yang (2009) [368] Symmetric Alternating Direction Method (SADM) Goldfarb et al. (2010) [90] Non Convex Splitting ADM (NCSADM) Chartrand (2012) [50] Douglas-Rachford Splitting Method (DRSM) Gandy and Yamada (2010] [86] Variant of Douglas-Rachford Splitting Method (VDRSM) Zhang and Liu (2013) [376] Proximity Point Algorithm (PPA) Zhu et al. (2014) [394] Proximal Iteratively Reweighted Algorithm (PIRA) Wang et al. (2014) [320] Alternating Rectified Gradient Method (ARGM) Kim et al. (2014) [158] Parallel Direction Method of Multipliers (PDMM) Wang et al. (2014) [319] Generalized Singular Value Thresholding (GSVT) Lu et al. (2014) [204] Generalized Accelerated Proximal Gradient (GAPG) He et al. (2013) [126] Improved alternating direction method (IADM) Chai et al. (2013) [46] Optimal Singular Values Shrinkage (OptShrink) Moore et al. (2014) [217] Iterative Thresholding with Primal-Dual Method (IT-PDM) Fan et al. (2014) [81] 2) Linearized solvers Linearized ADM (LADM) Yang and Yuan (2011) [352] Linearized ADM with Adaptive Penalty (LADMAP) Lin et al. (2011) [190] Linearized Symmetric ADM (LSADM) Goldfarb et al. (2010) [90] Fast Linearized Symmetric ADM (Fast-LSADM) Goldfarb et al. (2010) [90] Linearized IAD Contraction Methods (LIADCM) Gu et al. (2013) [97] 3) Fast solvers Randomized Projection for ALM (RPALM) Mu et al. (2011) [221] l1-filtering (LF) Liu et al. (2011) [194] Block Lanczos with Warm Start Lin and Wei(2010) [191] Exact Fast Robust Principal Component Analysis (EFRPCA) Abdel-Hakim and El-Saban (2012) [1] Inexact Fast Robust Principal Component Analysis (IFRPCA) Abdel-Hakim and El-Saban (2012) [1] Matrix Tri-Factorization (MTF) Liu et al. (2013) [199] Fast Tri-Factorization(FTF) Liu et al. (2013) [198] PRoximal Iterative SMoothing Algorithm (PRISMA) Orabona et al. (2012) [235] Fast Alterning Minimization (FAM) Rodriguez and Wohlberg (2013) [259] Fast Alternating Direction Method of Multipliers (FADMM) Yang and Wang (2014) [358] Fast Alternating Direction Method with Smoothing Technique (FADM-ST) Yang (2014)[356] Fast Randomized Singular Value Thresholding (FRSVT) Oh et al. (2015)[232] 4) Online solvers Online Alternating Direction Method (OADM) Wang and Banerjee (2013) [318] 5) Non convex solvers Difference of Convex (DC) Sun et al. (2013) [299] Fast Alternating Difference of Convex (FADC) Sun et al. (2013) [299] Non-convex Alternating Projections(AltProj) Netrapalli et al. (2014) [227] Iterative Shrinkage-Thresholding/Reweighted Algorithm (ISTRA) Zhong et al. (2015) [387] 6) 2D solvers Iterative method for Bi-directional Decomposition (IMBD) Sun et al. (2013) [301] Thierry Bouwmans et al. Categories Methods Authors -Dates Incremental Algorithms Recursive Robust PCP (RR-PCP) Qiu and N. Vaswani (2010) [242] Recursive Project Compressive Sensing (ReProCS) Qiu and Vaswani (2011 [243] Support-Predicted Modified-CS RR-PCP) Qiu and Vaswani (2011) [244] Support-Predicted Modified-CS) Qiu and Vaswani (2011) [243] Automated ReProCS Qiu and Vaswani (2012) [245] Practical ReProCS (Prac-ReProCS) Guo et al. (2013) [100] Incremental Low-Rank (iLR) Algorithm Wei et al. (2014) [331] Incremental PCP (incPCP) Rodriguez and Wohlberg (2014) [261] Online RPCA (ORPCA) Xu (2014) [345] Online RPCA via Stochastic Gradient Descent (ORPCA-SGD) Song et al. (2015) [322] Projection based RPCA (ProjectionRPCA) Lee and Lee (2015) [172] Online RPCA with Truncated Nuclear Norm Regularization (OTNNR) Hong et al. (2015) [128] Real time Implementations CAQR Anderson et al. (2010) [8] Real-Time PCP Pope et al. (2011) [240] LR Submatrix Recovery/Reconstruction (LRSRR) Guo et al. (2014) [106] Multi-Features Algorithms Multi-Features Algorithm (MFA) Gan et al. (2013) [85] Multi-Task RPCA (MTRPCA) Wang and Wan (2014) [324] Spatial-Temporal Algorithms Dense Optical Flow Gao et al. (2012) [87] Consistent Optical Flow Huang et al. (2013) [138] Smoothness and Arbitrariness Constraints (RFDSA) Guo et al. (2014) [108] Total Variation (TV) Regularizer Gao et al. (2015) [43] Piece-wise Low-rank Model Newson et al. (2015) [228] Graphs Model Shahid et al. (2015) [270] Fast Graphs Model Shahid et al. (2015) [271] Shape and Confidence Map-based (SCM-RPCA) Sobral et al. (2015) [292] Connectivity and Saliency Map (MODSM) Pang et al. (2015) [237] Salient Motion Detection Chen et al. (2015) [53] Compressive Sensing Algorithms Sparsity Reconstruction for Compressive Sensing (SpaRCS) Waters et al. (2011) [329] SpaRCS with Partial Support Knowledge (SpaRCS-PSK) Zonoobi and Kassim (2013) [395] Adaptive Reconstruction Compressive Sensing (ARCS) Yang et al. (2013) [350] LRSD for Compressive Sensing (LRSDCS) Jiang et al. (2014) [169] Recursive Low-rank and Sparse decomposition (rLSDR) Li and Qi (2014) [181] Optimal PCP Solutions Minimum Description Length (MDL) Ramirez and Shapiro (2012) [254] Saliency Measure Gao et al. (2012) [87] SVD Algorithms Full SVD -Partial SVD -Linear Time SVD Yang and An (2013) [357] smaller-scale SVD Zhang and Tian (2013) [378] block-SVD Chai et al. (2013) [46] Limited Memory SVD (LMSVD) Liu et al. (2013) [195] Symmetric Low-Rank Product-Gauss-Newton (SLRPGN) Liu et al. (2014) [196] The authors would like to thank the following researchers: Zhouchen Lin (Visual Computing Group, Microsoft Research Asia) who has kindly provided the solver LADMAP [190] and the l 1 -filtering [194], Shiqian Ma (Institute for Mathematics and Its Applications, Univ. of Minnesota, USA) who has kindly provided the solver LSADM [205], Congguo Tang (Dept. of Electical System Engineering, Washington Univ., USA) who has kindly provided the solver of RPCA-BPCP [306].

Decompositions

Contiguous Outlier Detection (DECOLOR) Zhou et al. (2011) [391] Direct Robust Matrix Factorization (DRMF) Xiong et al. (2011) [341] Direct Robust Matrix Factorization-Rows (DRMF-R) Xiong et al. (2011) [341] Probabilistic Robust Matrix Factorization (PRMF) Wang et al. (2012) [321] Bayesian Robust Matrix Factorization (BRMF) Wang et al. (2013) [366] Markov Bayesian Robust Matrix Factorization (MBRMF) Wang et al. (2013) [366] Practical Low-Rank Matrix Factorization (PLRMF) Zheng et al. (2012) [386] Low Rank Matrix Factorization with MoG noise (LRMF-MOG) Meng et al. (2013) [216] Unifying Nuclear Norm and Bilinear Factorization (UNN-BF) Cabral et al. (2013) [37] Low Rank Matrix Factorization with General Mixture noise (LRMF-GM) Cao et al. (2013) [42] Robust Rank Factorization (RRF) Sheng et al. (2014) [282] Variational Bayesian Method (VBMF-l1) Zhao et al. (2015) [384] Robust Orthogonal Matrix Factorization (ROMF) Kim and Oh (2015) [160] Contiguous Outliers Representation via Online Low-Rank Approximation (COROLA) Shakeri and Zhang (2015) [272] Online Low Rank Matrix Completion (ORLRMR) Guo (2015) [104] Matrix Factorization -Elastic-net Regularization (FactEN) Kim et al. (2015) [159] Incremental Learning Low Rank Representation -Spatial Constraint (LSVD-LRR) Dou et al. (2015) [73] Online Robust Low Rank Matrix Recovery (ORLRMR) Guo (2015) [105] Solvers Alternating Algorithm (AA) Zhou et al. (2011) [391] Block Coordinate Descent Strategy (BCDS) Xiong et al. (2001) [341] Conditional EM Algorithm (CEM) Jebara and Pentland (1999) [148] Augmented Lagrangian Multiplier (ALM) Zheng et al. (2012) [386] Alternative Direction Descent Algorithm (ADDA) Sheng et al. (2014) [282]

Table 18 Robust Low Rank Minimization: A Complete Overview. The first column indicates the concerned category and the second column the name of each method. Their corresponding acronym is indicated in the first parenthesis. The third column gives the name of the authors and the date of the related publication.

where /alpha and /beta are regularization parameters. Then, Zhou et al. [391] defined the foreground support matrix of S, denoted F ∈ {0, 1} m×n as follows: 148). That is:

Let P F (X) represents the orthogonal projection of the matrix X onto the linear space of matrices supported by F :

and P F (X) is its complementary projection, i.e P F (X) + P F (X) = X. Thus, the Equation ( 152) is obtained: