Online Robust Low Rank Matrix Recovery

1 Introduction

Low Rank Matrix Recovery (LRMR) aims to seek a lowdimensional subspace from given data such that the residual between the given and the recovery is minimized with respect to some optimality criterion [Recht et al., 2010; Chandrasekaran et al., 2011]. Its benefit, besides being a standard tool for dimensionality reduction, has been witnessed by a wide spectrum of low level signal processing applications, such as image denosing [Gu et al., 2014], colorization [Wang and Zhang, 2012], reflection separation [Guo et al., 2014a] and rectification [Zhang et al., 2012], as well as higher level tasks, like collaborative filtering [Jamali and Ester, 2011; Park et al., 2013; Shi et al., 2013; Chen et al., 2014], visual salience detection [Shen and Wu, 2012], video editing [Guo et al., 2013] and dictionary learning [Mairal et al., 2010].

The quality of LRMR largely depends on the choice of optimality criteria for different tasks. Principle Component Analysis (PCA) [Pearson, 1901] adopts the ℓ2 loss (also known as the squared loss) by assuming that the residual existing in the observation follows a Gaussian distribution. Although the ℓ2 is arguably the most commonly used loss, its performance may sharply degenerate when dealing with data polluted by outliers, which is frequently encountered in realworld applications. To be robust against outliers, the ℓ0 loss is the ideal option. Unfortunately, the non-convexity, nonsmoothness and difficulty of being approximated (NP-hard) make the ℓ0 unsuitable for practical use. Alternatively, its tightest convex surrogate, the ℓ1 , is employed to achieve the robustness, which corresponds to the least absolute deviations technique and is actually optimal for Laplacian noises. In this way, Cand`es et al. proved that their proposed method, namely PCP [Cand`es et al., 2011], can exactly recover the underling low rank matrix from the observation even with the fraction of outlier up to about 1 3 under some conditions. In parallel, there are a few works carried out from probabilistic standpoints. Probabilistic Matrix Factorization (PMF) [Salakhutdinov and Mnih, 2008] and Probabilistic Robust Matrix Factorization (PRMF) [Wang et al., 2012] are two representatives, the residual penalties of which are equivalent to using the ℓ2 (PCA) and ℓ1 (PCP), respectively. Although the techniques above have made great progresses in LRMR, the tolerance to outliers is expected to be further improved.

Another issue should be concerned is the applicability to large-scale data. Conventional solvers of LRMR, like PCA [Pearson, 1901], PCP [Cand`es et al., 2011], PMF [Salakhutdinov and Mnih, 2008] and PRMF [Wang et al., 2012] process data in batch mode, which are both computationally expensive and memory exhausting. Moreover, they are inflexible to incrementally add samples. To overcome the shortcomings of batch mode, several online strategies have been developed recently. GROUSE [Balzano et al., 2011], an online extension of PCA, is derived by analyzing incremental gradient descent on the Grassmannian manifold of subspaces. GRASTA [He et al., 2012] follows the same technical line as GROUSE with further consideration of outliers, which can be seen as an online version of PCP. Wang et al. modified PRMF as Online PRMF for sequentially processing samples [Wang et al., 2012]. Online Robust PCA (OR-PCA) [Feng et al., 2013] has proven that it converges to the same optimal solution as PCP [Cand`es et al., 2011] does asymptotically. The methods mentioned above indeed enable the traditional schemes with the scalability, but the improvement space exists, especially in terms of the robustness to more corruptions.

This work proposes an online technique to robustly recover the low rank component from the observation with heavy outliers. Please notice that if the outlier support is given or precisely estimated, the problem reduces to the completion task that is much easier to conquer. Specifically, we first develop an online algorithm to solve the low rank matrix completion. Then, we step forward to a more challenging problem, say online low rank matrix recovery with robust outlier support estimation. As regards residual modeling, we employ a Gaussian component to fit the noise while a uniform distributed component to host the oultier. With this strategy, we are capable to handle more severely contaminated data. Extensive experiments on both synthetic and real data are conducted to verify the effectiveness of the proposed technique and demonstrate its advantages over the state-of-the-arts.

2 Our Method

Suppose we are given a finite set of data vectors Y = [y1, ..., yn] ∈ Rm×n . It will be helpful to consider that low rank matrices can be factorized by Y = UV T , where U ∈ Rm×r and V ∈ Rn×r . In such a way, the rank of recovered matrices is guaranteed to be never over r. This work formally concentrates on the following general problem:

To prevent U from being with arbitrarily large values, the common manner is to enforce the ℓ2 norm of each column to be less than or equal to one. For the rest of the paper, the discussion is always under this constraint.

Having W = [w1, ...,wn] ∈ {0, 1}m×n that indicates which elements of Y are observed, the associated problem is designated as Low Rank Matrix Completion (LRMC). The difficulty, however, greatly increases when the data are ruined by outliers, as the outlier support is unknown, we call this problem Robust Low Rank Matrix Recovery (RLRMR). In the next subsections, we will first develop an online algorithm based on stochastic optimization to handle LRMC, and then design a robust outlier support estimation scheme for solving RLRMF in an online fashion.

2.1 Online LowRankMatrix Completion

With the data matrix Y and its corresponding support W, solving the problem (1) is indeed equivalent to minimize the following empirical cost function:

where the loss function for each sample is defined as:

As pointed out by Bottou and Bousquet [Bottou and Bousquet, 2008], one may be interested in the minimization of the expected cost instead of the empirical cost, as the empirical is just an approximation of the expected. The corresponding expected cost function is:

1 Here, "a.s." means almost sure convergence.

where the expectation is taken with respect to the distribution of the samples and the known supports. However, classical batch optimization methods become impractical in terms of speed and/or memory requirements to do so. As a consequence, online techniques are desired.

In the following, we first establish an online algorithm to minimize the empirical cost function (2). At the end of this subsection, we show that the proposed algorithm converges to the solution of the expected cost function (4) asymptotically. When a new sample yt arrives, two things have to be dealt with, i.e. computing vt and updating U (t) . With the help of alternating minimizing strategy, we first focus on vt, which associates with the following optimization problem with fixed U (t-1) (obtained from previous iteration):

where Dt ≡ Diag(wt). As can be seen from (5), it is a classic least squares problem and has a closed-form solution:

where (DtU (t-1) )† denotes the Moore-Penrose pseudoinverse of DtU (t-1) .

For updating U in the tth time instance, say U (t) , we minimize the following objective function:

We can observe that St(U) is quadratic in U and aggregates the past data with a few sufficient vi. It is obvious that St(U) upper-bounds Ft(U) defined in (2). Similar to [Mairal et al., 2010], we can adopt U (t-1) as the warm restart for calculating U (t) . Because St-1(U) and St(U) become closer as t gets larger, and so do U (t-1) and U (t) . This inspires us to employ the block-coordinate decent with warm restarts [Bertsekas, 1999] to update U efficiently. Moreover, considering some entries in the samples may be missing leads us to handle different rows of U separately based on the visibility of the corresponding elements, for the sake of information balance. In other words, if an element in a sample is invisible, it should contribute nothing to the update of the corresponding row of U. Derived from (7), for each row Uj·, we have:

where [yij ]j denotes the j th entry of vector yij , tj represents the number of times that the j th row of samples is observed, and ij is the refreshed observation index. The details of solving the problem (7) are provided in Algorithm 1. Please note that, Aj (t) ≡ ∑tj ij=1 vij vTij , and bj (t) ≡ ∑tj ij=1 [yij ]jvij . For clarity, we summarize the whole procedure to solve the problem (2) in Algorithm 2.

Algorithm 1: The Basis Update

Algorithm 2: Online Low Rank Matrix Completion

Input: The observed data that will be revealed sequentially {y1 , ..., yT }, and the corresponding support {w1, ..., wT }, U0 ∈ Rm×r . ∀k ∈ {1, ..., m} Ak ← 0 ∈ Rr×r ; bk ← 0 ∈ Rr ; for t = 1 to T do vt ← argminvt ‖wt ⊙ (yt -Ut-1vt)‖22 ; for k = 1 to m do Ak (t) ←Ak (t-1) + [wt]k ×vtvTt ; bk (t) ← bk (t-1) + ([wt]k × [yt ]k)vTt ; end Update U (t) using U (t-1) as warm restart by Algorithm 1; end Output: U (T)

Prior to analyzing the convergence of the designed online solver for LRMC, we state some assumptions under which our analysis holds. The first assumption is that the data are bounded, which is natural for the realistic data, such as audio signals, images and videos. The second is the uniqueness of the solution to vt (5) is satisfied. The solution is unique if and only if DtU (t-1) is full column rank. In practice, this condition is easy to meet as r ≪ m and the fraction of missing elements in yi is not extremely large typically. The third is that the quadratic surrogate functions St(U (t) j· ) are strictly convex with lower-bounded Hessians. This hypothesis is experimentally verified after processing several data points with reasonable initialization (we will see the random initialization also works sufficiently well in Sec. 3), which has also been experimentally verified on the task of online dictionary learning [Mairal et al., 2010]. The convergence property of the proposed Algorithm 2 is given by Theorem 1.

Theorem 1. Assume the observations are always bounded. Given the rank of the optimal solution to (4) provided as r, and the solution U (t) ∈ Rm×r obtained by Algorithm 2 is full column rank, then U (t) converges to at least a stationary point of the problem (4) asymptotically.

Proof Sketch. From formula (7), we can obtain a set of row-wise sub-problems as given in (8), which are independent to each other. It indicates that the theorem can be recognized if each sub-problem holds. Please note that the separated sub-problems have the same form as the problem (6) shown in [Mairal et al., 2010], and thus can be proved in the way offered by [Mairal et al., 2010]. Due to space limit, instead of the complete proof, we here only provide the proof sketch including three steps: 1) the surrogate function St(U (t) j· ) converges almost surely; 2) the variations of U (t) j· are asymptotic, i.e. ‖U (t+1) j· - U (t) j· ‖F = O( 1t ); and 3) F(U (t) j· ) -St(U (t) j· ) → 0 almost surely.

2.2 Online Robust LowRankMatrix Recovery

In many real world applications, the observed data are very likely contaminated by outliers, the support of which is usually unknown in advance. For the purpose of accurate recovery of underling matrices, it is essential to eliminate or reduce the disturbance from these outliers. Now we move on from OLRMC to its upgraded problem, say Online Robust Low Rank Matrix Recovery (ORLRMR). It is sure that, if the outlier support is at hand, ORLRMR degenerates to OLRMC. Thus, compared with OLRMC, ORLRMR requires to detect outliers from observations automatically and robustly.

Recall that each entry [yi]j(i = 1, ..., n; j = 1, 2, ...,m) of the sample yi can be modeled as [yi]j = Uj·vi + [ei]j. Suppose observations are always bounded, an outlier can be with any value in the range with equal probability instead of either Gaussian (ℓ2-PCA and PMF) or Laplacian (ℓ1 - PCP and PRMF). Therefore, we assume the outliers follow the uniform distribution 1ϵ . As for the residuals caused by other factors, we simply assume they (approximately) follow a Gaussian distribution N(0, σ2 i ). Let πu i and πg i , such that πu i + πg i = 1, respectively denote the percentages of outlier and the other in the ith sample, which are also unknown in advance. As a result, each [ei]j can be viewed as a sample from a mixture model of distributions with probability p([ei]j) = πu i N(0, σ2 i ) + πg i 1ϵ . The likelihood of Y can be written in the following form:

where Θi = {σ2 i , πu i , πg i } is a parameter vector. The nonnegative log-likelihood function of (9) is:

So far, the outlier support has not been explicitly touched upon. But, as can be seen from (10), it is a typical mixture model that can be immediately minimized by ExpectationMaximization algorithms. This reminds us that there involve latent variables [zki ]j ∈ {0, 1} with ∑Kk=1 [zki ]j = 1, which indicate the assignment of [ei]j to a specific component of the mixture containing K models in total. In this task, zgi that corresponds the Gaussian component, can be employed as the desired support. To be explicit with respect to the latent variables, and make the model more compatible to additional constraints, by applying Proposition 1, we obtain the

Algorithm 3: Online Robust LRMR

cost function like:

Proposition 1. Given two functions πk(x) > 0 and pk(x) > 0, we have the following [Liu et al., 2013; Guo et al., 2014b]:

where Φ(x) = {Φ1(x), ...,ΦK(x)} are hidden variables, and ∆+ = {Φ(x) : 0 < Φk(x) < 1, and ∑Kk=1 Φk(x) = 1} is a convex relaxation of a characteristic function decomposition.

Proof. The proof can be easily done via Lagrangian Multiplier method. Due to the limited space, we omit it here.

With the modification on (10), the cost function (11) is explicit to the latent variables Φ, and more importantly, is ready to take into account extra constraints on Φ, such as sparsity of outlier and smoothness of support. In the following, for instance, we introduce the sparsity of outlier into (12), i.e. ‖Φu‖0, which is a commonly used prior to a variety of real cases. It is well known that the ℓ0 norm is non-convex and usually replaced with the ℓ1 . As the elements in Φ distribute between 0 and 1, by slight algebraic transformation, we have the following formula:

where λ is a coefficient controlling the weight of outlier sparsity. Please note that we use ˆ C to present ˆ C(U,V ,Θ,Φ) for brevity. Following the same spirit of solving the problem (2), we can seek the solution of (12) in an online manner.

When a new data vector yt comes, we first calculate Φt, Θt and vt over the previously updated U (t-1) by the alternating minimizing strategy as follows:

To further enforce the summation of [Φut ]j and [Φgt ]j to be 1, simply normalizing them gives:

It is worth noting that the above processing corresponds to the E step of EM. Next, we focus on updating the parameters of the mixed model including σ2 t , πu t and πg t , which can be achieved via [Dempster et al., 1977]:

In turn, vt can be obtained through optimizing the following:

where [Ωt]j ≡ √[Φgt ]j. This is in the same form with (5), and can be solved in the same way. Please notice that the updating of Θt and vt can be viewed as the M step of EM. Φt, Θt and vt are iteratively updated via the procedure above until converge. Moreover, the element of Φt is in (0, 1) instead of binary, although very close to either 0 or 1. To make the weight Φgt consistent with the binary support (LRMC), we pre-define a threshold (e.g. 0.5) to binarize Φgt as ˆ wt. With the well constructed support ˆ wt and vt, we are ready to update the basis U. For completeness, the solver of ORLRMF is summarized in Algorithm 3.

3 Experimental Verification

In this section, we conduct experiments on synthetic data to reveal the convergence speed and the robustness to outliers of the proposed ORLRMR, and on real data for showing the advantages of our method compared against alternative approaches. Because OLRMC is a special case of ORLRMR, we concentrate on ORLRMR. Please note that ORLRMR has one free parameter λ for balancing the outlier sparsity. From Eqs. (12) and (13), we can see that the effect of λ is actually to pull down 1ϵ to 1 ϵ exp(λ) . No doubt that when λ grows, the probability for a residual of being assigned to the outlier component decreases. But, the margin between the tail of Gaussian distribution and 1ϵ is usually large, therefore ORLRMR can perform stably with a relatively wide range of λ, especially for Gaussian noises with small variances. For the experiments shown in this paper, we uniformly set λ to 2. The experiments are carried out on a MacBook Pro running OS X 64bit operating system with Intel Core i.7 2.8GHz CPU and 16GB RAM. Our code is implemented in Matlab.

3.1 Synthetic Data

Data Preparation. For fair comparison, we generate data in a similar way with PCP [Cand`es et al., 2011] and ORPCA [Feng et al., 2013]. That is, a set of n clean data points is produced by Y ∈ Rm×n = UV T , where U ∈ Rm×r and V ∈ Rn×r . The entries of both U and V are i.i.d. sampled from the N(0, 1/n) distribution, and then the data are projected into [-1, 1]. The basis of the generated subspace is U, and r is the intrinsic rank (dimension) of the spanned subspace. As for the corruption, we contaminate the data by replacing a fraction po of Y with outliers drawn from a uniform distribution over the interval of [-1, 1].

Task and Quantitative Metric. We evaluate the ability of correctly recovering the subspace of corrupted observations under varying settings of the error density and intrinsic subspace rank. The U (0) that feeds to the competitors is randomly initialized in the same fashion as data generation. To measure the similarity between the recovered subspace ˆU and the ground-truth U, the Expressed Variance (E.V.) [Xu et al., 2010; Feng et al., 2013] is employed as our metric, the definition of which is:

where orth(·) and tr(·) stand for the orthogonalizing and trace operators, respectively. The averaged results from 10 independent trials are finally reported.

Competitors. We compare our ORLRMR with recently proposed online techniques, i.e. GRASTA and OR-PCA. The code of GRASTA is downloaded from the authors' website, the core of which is implemented in C++. Its most important parameter is the sub-sampling ratio, to obtain its best possible results, we turn off the sub-sampling to take into account whole information. The code of OR-PCA is not available when this paper is prepared. So we implement it in Matlab by strictly following the pseudo code given in [Feng et al., 2013], and adopt the parameters suggested by the authors.

Performance Analysis. We provide the performance curves of ORLRMR, OR-PCA and GRASTA against the number of samples in Fig. 1, under the setting of m = 400 and n = 104 , to see the robustness of the contestants to corruptions and how the performance is improved when more samples are revealed. As can be seen from the pictures, all the three methods achieve better results with more observations. For the four cases, the performance of GRASTA is much inferior to the others, possibly because the intrinsic dimensions versus the ambient are relatively large for GRASTA, as well as the basis updating by GRASTA only takes one gradient descent step while both ORLRMR and OR-PCA have closed-form solutions. We will see that GRASTA attains reasonably well results when the intrinsic rank is relatively small later. In the case of intrinsic rank r = 80 and outlier ratio po = 0.2, the plots of OR-PCA and ORLRMR are very close, both of which reach about E. V. = 0.9 after processing 104 samples. As the outlier ratio increases, the convergence speed drops. Even though, ORLRMR gives E.V. ≈ 0.7 and E.V. ≈ 0.6 for po = 0.3 and po = 0.4, respectively, which outperforms ORPCA by about 0.2. This indicates that ORLRMR can bear more corruptions than the others. In addition, we show a case of r = 120 and po = 0.3 to view the effect of r. The curves tell us that OR-PCA and ORLRMR perform similarly with those of r = 80 and po = 0.3 with very slight depression, which means the change of r influences not much to them. Please note that, the E.V. of the initialized basis increases as the intrinsic rank approaches to the ambient dimension.

3.2 Real Data

Dataset. The Star dataset consists of 9 real world surveillance videos, which has a variety of scenarios including Escalator (Esc: indoor, dynamic background, crowd), Hall (Hal: indoor, static background, crowd), BootStrap (BS: indoor, static background, crowd), Lobby (Lob: indoor, light switch, people), ShoppingMall (SM: indoor, static background, people), Campus (Cam: outdoor, dynamic background, cars), WaterSurface (WS: outdoor, dynamic background, lingering person), Fountain (Fou: outdoor, dynamic background, people), and Curtain (Cur: indoor, dynamic background, people). Task and Quantitative Metric. This task is to separate foregrounds from surveillance videos, in which the ground-truth foreground is given but the ground-truth background model is not available. Moreover, different to the simulation, the background may comprise dynamic factors, such as fountains and waving curtains. To assess the performance on foreground detection, we use the Receiver Operating Characteristic (ROC) curve as our metric. As the rank of the background of surveillance videos is typically small, we empirically set it to 5 for all the 9 sequences.

Competitors. Besides GRASTA and OR-PCA, we further employ a batch method, i.e. PRMF [Wang et al., 2012], as our reference. PRMF is a probabilistic approach to robust matrix factorization, which holds an equivalent relation with PCP. PCP [Cand`es et al., 2011] is not used, mainly because its costs in both memory and time are too expensive to process long sequences. The Matlab code of PRMF is downloaded from the authors' website, the parameters are set as the values provided in the implementation package released by the authors. In addition, to obtain best possible performance of GRASTA and OR-PCA, the first 5 frames of each sequence are used to initialize the basis for GRASTA and OR-PCA 2 . To show the insensitivity of ORLRMR to initialization, its basis is randomly generated as the previous experiment do.

2 Without a good initialization, OR-PCA can not give reasonable results in this experiment. The reason possibly is that the subspace rank versus the ambient dimension is too small and the number of samples is inadequate for converging.

Performance Analysis. Figure 2 summarizes the ROC plots for the 9 videos. From the top row, we can see that our ORLRMR outperforms GRASTA, OR-PCA and even the batch approach PRMF. The reason is that these three videos are in presence of crowd, and ORLRMR can handle heavy outliers better than the others. For the rest videos, PRMF attains the best results (except for the Campus sequence, in which GRASTA slightly exceeds PRMF), while ORLRMR competes very favorably with PRMF in most cases. The only exception is the Curtain video, the background of which contains intensively moving curtains. Although the ground-truth foreground excludes these curtains, from the perspective of low rank recovery, they are more like outliers, and ORLRMR treats these as outliers. Please notice that GRASTA performs reasonably well on these real videos, since the subspace dimension is relatively small. As regards OR-PCA, its inferior performance for most cases may be because of the large fractions of foreground and/or the insufficient number of streaming samples for achieving convergence. Together with the simulation, we can conclude that ORLRMR is better at dealing with heavily polluted data and is of faster convergence speed than OR-PCA. In terms of time, PRMF costs, without surprise, the most (5.7 frames per second on average over the 9 sequences) among the contestants, due to its batch nature. ORLRMR achieves 48.3 fps, which outperforms ORPAC 40.6 fps and is competitive with GRASTA 52.9 fps, even though the core of GRASTA is implemented in C++. The fps numbers corresponding to each sequence can be found in the legends of Fig. 2. We would like to emphasize that the experiments in this part do not intend to verify our ORLRMR wins over those special designs on the task of foreground detection. The aim is to demonstrate the merits of our general model without involving domain knowledge. It is certain that equipping with proper priors for specific problems can further improve the performance of any of the general models like GRASTA, OR-PCA, ORLRMR and PRMF.

4 Conclusion

In this paper, we have proposed to address the low rank matrix completion and recovery tasks in an online fashion, which scales up gracefully to large data sets. The mixture model consisted of a Gaussian and a uniform distributed components has shown its significant improvement on the robustness to heavy outliers. The theoretical and experimental results have demonstrated the properties and advantages of our scheme compared alternative methods. It is worth noting that our technical framework is ready to embrace additional specific constraints for further boosting the performance on different tasks. For instance, the spatial smoothness can be introduced to accurately detect foregrounds in surveillance videos. It is positive that our proposed technique can be widely applied to tabulated data processing applications.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant 61402467 and in part by the Excellent Young Talent Program through the Institute of Information Engineering, Chinese Academy of Sciences.