ROUTE: Robust Outlier Estimation for Low Rank Matrix Recovery

1 Introduction

Low rank matrix recovery (LRMR) is a process of discovering underlying structures from given measurements, the inspiration and motivation of which are both that, in real cases, even very high-dimensional observations should be from a low-dimensional subspace but unfortunately with intrusion of outliers and/or noises. As a theoretic foundation in computer vision, pattern recognition and machine learning, the effectiveness of LRMR has been confirmed by numerous fundamental tasks, such as principal component analysis [Pearson, 1901; Cand`es et al., 2011], collaborative filtering [Zhang and Wang, 2016; Chen et al., 2014] and subspace clustering [Nie and Huang, 2016; Liu et al., 2013], as well as a wide spectrum of applications, like image denoising [Gu et al., 2014], reflection separation [Guo et al., 2014] and super-resolution [Jing et al., 2015], to name just a few.

Formally, the LRMR problem can be directly or indirectly written in the following form:

where Y ∈ Rm×n , L ∈ Rm×n and E ∈ Rm×n designate the given data, the desired structure and the error residue, respectively. The function Ψ(·) is a penalty on the residual between the observed and recovered signals, rank(·) stands for the low-rank constraint, and α is a non-negative parameter that provides a trade-off between the recovery fidelity and the low-rank promoting regularizer. Furthermore, PΩ(·) is the orthogonal projection operator on the support Ω ∈ {0, 1}m×n . From Eq. (1), we can find that the quality of recovery depends on both the models of rank(L) and Ψ(E).

As one of the two pivotal factors in LRMR, a proper lowrank promoting constraint on L is required to advocate the expected structure. It is computationally intractable (NPhard) to directly minimize the rank function, say rank(L), due to its non-convexity and discontinuity. A widely used scheme is employing its tightest convex proxy, i.e. the nuclear norm ‖L‖∗ [Recht et al., 2010; Cand`es et al., 2011; Zhou et al., 2013]. Nuclear norm minimization (NNM) based approaches can perform stably without knowing the target rank of recovery in advance. But, their applicability is often limited by the necessity of executing expensive singular value decomposition (SVD) for multiple times. At (much) less expense, bilinear factorization (BF) [Eriksson and van den Hengel, 2010; Srebro and T.Jaakkola, 2003; Meng and De la Torre, 2013; Salakhutdinov and Mnih, 2008; Lakshminarayanan et al., 2011] is an alternative by replacing L with UV, where the product of two factor matrices U ∈ Rm×r and V ∈ Rr×n implicitly guarantees that the rank of UV is never over r, typically r ≪ min(m, n). This factorization strategy, through getting rid of SVDs, can greatly release the pressure of computation and provide accurate results when the target rank is given. Unfortunately, in some tasks, the target rank is unknown beforehand. In such a situation, the performance of BF would sharply degrade because of its sensitivity to the guess of target rank, especially when the data are severely contaminated. For bridging NNM and BF, and inheriting their respective merits, some bridges are recently built [Zheng et al., 2012; Wang et al., 2012]. One representative is adding ‖U‖2F + ‖V‖2F into the objective of the BF [Cabral et al., 2013]. Although the techniques above have made great progresses, the tolerance to dirty data is expected to be further improved.

In practical scenarios, acquiring perfect data is never the case. Furthermore, "a little gall spoils a great deal of honey" is quite a common issue. This is to say, without an effective strategy to reduce the negative effect from outliers and/or noises, the low rank matrix recovery is very likely prevented from reasonable solutions. Hence, besides properly modeling the low-rank structure, how to handle the pollution, especially gross outliers, is core to the performance of recovery. Arguably, the square loss (a.k.a. ℓ2 loss) is the most commonly used penalty, which is optimal to Gaussian noises, like PCA [Pearson, 1901]. But, the square loss lacks robustness to outliers that are not unusual to find in real data. To be robust against gross corruptions, the ℓ1 loss becomes popular, e.g. Robust PCA [Cand`es et al., 2011]. For better fitting the residual, combining the ℓ1 loss for the outlier component and the ℓ2 for the small noise one is considered in [Zhou et al., 2010]. Although the ℓ1 loss can perform better than the ℓ2 in dealing with the outlier, it still suffers from the scale issue. The ideal option to model outliers is the ℓ0 loss, due to its scale invariance. The non-convexity and discreteness of the ℓ0 penalty make it not so preferred by the community, although many works have proven that the ℓ0 loss can improve the performance on different tasks, like [Pan et al., 2017] for image deblurring, [Zhou et al., 2013] for foreground detection and [Guo et al., 2013] for video editing.

Back to the general formulation (1), if the support of both outliers and missing elements is given, the problem turns out to be a simpler version, i.e. the low rank matrix completion (LRMC). Compared with LRMC, the difficulty of LRMR, because of the unknown outlier support, significantly increases, which corroborates the intuition and theoretical fact that knowing the corruption location is beneficial. Therefore, it is natural to ask that if we can connect the LRMR to the LRMC via robustly estimating outliers, since by doing so the LRMR will be conquered more easily.

To answer the question, this paper proposes a Robust OUTlier Estimation method, called ROUTE, for recovering low rank matrices. More concretely, the contributions of this work are summarized as follows:

1. We design a method to jointly estimate outliers and recover the low rank matrix, namely ROUTE-LRMR, which unifies the LRMR and LRMC by treating both the missing and estimated outliers as weights;
2. Compared with the hard binary support, our weighting scheme assigns real-valued weights [0, 1], which can be viewed as classification with confidence/probability;
3. Our design employs a maximum entropy regularization term to minimize the prediction bias, which behaves like a sigmoid function, arguably the most suitable classification function;
4. To seek the optimal solution for ROUTE-LRMR, we customize an Alternating Direction Minimization based algorithm. Theoretical analysis together with experimental results on both synthetic and real data are provided to show the efficacy of our ROUTE and reveal its superiority over other state-of-the-arts.

2 Methodology

2.1 ProblemFormulation

In the simplest case, the support of observed elements is at hand, and the data are clean or just with slight noises. An option for recovering the low rank component (LRMC) is to optimize the following problem:

where ◦ is the Hadamard product operator. As mentioned, the nuclear norm minimization requires to execute expensive SVDs on the whole data. To mitigate the computational pressure, Theorem 1 provides a bridge between NNM and BF models.

Theorem 1. For any matrix L ∈ Rm×n , the following relationship holds [Mazumder et al., 2010]:

If rank(L) = r ≤ min(m, n), then the minimum solution above is attained at a factor decomposition L = UV, where U ∈ Rm×r and V ∈ Rr×n .

In the sequel, applying Theorem 1 on (2) reads:

Compared with directly minimizing ‖Ω ◦ (Y - UV)‖2F , the model (3) inherits the advantage of (2), which avoids overfitting when r is larger than the intrinsic rank.

In the real world, however, the data are frequently polluted by, besides small noises, gross corruptions, which may prevent the recovery from reasonable results. Hence, some steps should be taken for reducing the negative effect of such pollution. Recall that the ℓ0 loss is ideal to host outliers, with slight modification, we have:

where β is a weight to the corresponding term and 1 represents an all-one matrix with comparable size. We can see that from Eq. (4), the support Ω is replaced by a weight matrix W that can contain both the given support and the estimated outlier support. Please note that, under the binary weighting, ‖W ◦ (Y - UV)‖2F = ∑i,j wij[Y - UV]2ij and ‖W‖1 equals to ‖W‖0 that imposes the sparsity on the outliers.

The hard weighting, for one thing, frequently leads the optimization to be stuck into bad local minima. For another thing, the pollution in data is often non-homogeneously distributed. To address the discreteness issue and reflect the importance of elements more faithfully, we employ an entropy term. The definition of entropy is - ∑kc=1 pc log pc with ∑kc=1 pc = 1. The principle of maximum entropy tells that, the probability distribution which best represents the current state of knowledge is the one with largest entropy subject to accurately stated prior data. In other words, it is able to minimize the prediction bias. Return to our problem, the weighting variable wij can be equally viewed as the probability of the corresponding entry being classified as an outlier. It is instructive to note that maximizing the entropy (concave) is equivalent to minimizing its negative (convex). Consequently, we have the final formulation of ROUTE-LRMR as follows:

where γ is a non-negative coefficient controlling the importance of the corresponding term. Further, due to the relaxation, √Wwith entries √wij is used to hold the equivalence:

∑i,j wij[Y - UV]2ij = ‖√W ◦ (Y - UV)‖2F . For (4),

2.2 Optimization

As we have seen in (5), it has embraced all the aforementioned concerns for simultaneously pursuing outliers and recovering the low rank matrix. The Augmented Lagrange Multiplier (ALM) with Alternating Direction Minimization (ADM) scheme [Lin et al., 2011] has proven to be an efficient and effective solver for problems like (5). To apply ALM-ADM on our problem, the objective is required to be separable. To this end, an auxiliary variable L is introduced to replace UV in the third term. Accordingly, L = UV performs as an additional constraint. It is worth noting that the constraints on W and W are enforced as hard constraints. The augmented Lagrangian function of (5) is defined as:

where 〈·, ·〉 designates the inner product, µ is a positive penalty and Z is a Lagrangian multiplier. The solver updates the variables in an iterative manner. For ease of exposition, we split the variables, except for Z and µ, into two groups, including Group I: {U,V,L} and Group II: {W,W}, connected by L.

Group I - Dropping the unrelated terms yields:

For all of U, V and L, their respective solutions in closedform are calculated via equating the derivatives of (7) in U, V and L to zero:

where I means the identity matrix with proper size, C(p+1) stands for I+µ(t)UT(p+1)U(p+1), and the division in updating L is element-wise.

Group II - Picking out the terms relevant to W and W results in the following optimization problem:

From the objective of (9), we find that the problem can be decomposed into a set of independent sub-problems. Now, without any loss of generality, let us take the (i, j)-th element for example. Casting the problem into the Lagrange Multiplier framework gives the following Lagrange function:

where λi is a Lagrange multiplier. Taking the derivative of Q(wi,wi, λi) to wi, wi and λi respectively and setting them to zero gives the optimal solution to wi as:

which is in a standard sigmoid form. And its complementary wij(p+1) ← 1 - wij(p+1).

Remarks (a) When wij ∈ {0, 1} adopted and the entropy term disabled (hard weighting), the solution to Eq. (9) is: if α2 [Y - L]2ij < β, then wij ← 1; if α2 [Y - L]2ij = β, then wij could be either of {0, 1}; otherwise wij ← 0. (b) When wij ∈ [0, 1] adopted and the entropy term disabled (relaxed version), the solution to Eq. (9) is: if α2 [Y-L]2ij < β, then wij ← 1; if α2 [Y-L]2ij = β, then wij could be any value in [0, 1]; otherwise wij ← 0.

Multiplier and µ - The Lagrange multiplier Z and µ are updated via:

The parameter µ is monotonically increased by ρ during iterations, leading the solution to the feasible region.

Algorithm 1: ROUTE-LRMR: Solver to Eq.(5)

Input: Observation matrix Y ∈ Rm×n ; support Ω ∈ {0, 1}m×n ; a guess/target rank r; non-negative parameters α, β and γ. Init.: µ(0) ← 1 and ρ ← 1.1; W(0) ∈ Rm×n ← Ω ◦ 1; L(0) ∈ Rm×n , U(0) ∈ Rm×r and V(0) ∈ Rr×n are all initialized randomly; Z(0) ∈ Rm×n ← 0; t ← 0. while not converged do p ← 0; while not converged do Update U(p+1), V(p+1) and L(p+1) via (8); for ∀ (i, j) & Ωij do Update wij(p+1) via (11); end p ← p + 1; end {L(t+1) , U(t+1) , V(t+1) } ← {L(p), U(p), V(p)}; Update Z(t+1) and µ(t+1) via (12); t ← t + 1 end Output: Optimal W∗ and L∗

For clarity and completeness, the procedure of solving (5) is outlined in Algorithm 1. The algorithm should not be terminated until the equality constraint L = UV is satisfied up to a given tolerance, that is ‖L -UV‖F ≤ ς‖Y‖F , or the maximal number of iterations is reached. In all our experiments, the tolerance factor ς is chosen as 1e-7. Please refer to the complete Algorithm 1 for other details that we can not cover in the text.

3 Theoretical Analysis

We first provide some useful theoretical results, including Lemma 1 and Proposition 1 for Group II, as well as Theorem 2 for the inner loop of Algorithm 1.

Lemma 1. At stage p with L(p) fixed, the solution to Eq. (9) (Group II), i.e. wij given in Eqn.(11), is global optimal to the corresponding intermediary problem.

Proof. First, having L(t) fixed, the objective function in (9) is convex with respect to wij ∈ [0, 1]. The solution in Eq.(11) is computed by the Lagrange multiplier method, which guarantees that the obtained solution is feasible and satisfies the KKT conditions for the convex problem (9). Thus, we reach the conclusion.

Proposition 1. The function defined in Eq. (9), containing three parameters including ˜ β := β/α, ˜ γ := γ/α and εij := [Y-L]2ij , has the following properties:

1. wij(˜ β, ˜ γ, εij) is monotonically decreasing with respect to εij, which holds limεij→0 wij(˜ β, ˜ γ, εij) = 1 1+exp(- ˜ β/˜ γ) and limεij→+∞wij(˜ β, ˜ γ, εij ) = 0;
2. wij(˜ β, ˜ γ, εij) is monotonically increasing with respect to ˜ β, which holds that lim˜ β→0 wij ( ˜ β, ˜ γ, εij) = 1 1+exp(εij /˜ γ) and lim˜ β→+∞ wij ( ˜ β, ˜ γ, εij ) = 1;
3. wij(˜ β, ˜ γ, εij) is an inverse-'S' shaped function, which approaches a binary function when ˜ γ → 0 and the constant 1/2 when ˜ γ → +∞.

Each statement takes care of one target parameter with the others fixed to be constants.

Proof. It can be easily verified by the definition.

Theorem 2. At stage t with Z(t-1) fixed, the inner loop of Algorithm 1 is guaranteed to converge to a partial minimum (either a stationary point or a local minimum) of the corresponding intermediary problem.

Proof. The updating of U, V, L andW in the inner loop follows the manner of alternate convex search (ACS) [Bazaraa et al., 1993]. As can be seen from Eq. (8) together with Lemma 1, the KKT conditions to all the involved variables are satisfied, which is sufficient to draw the conclusion [Gorski et al., 2007].

Next, we shall consider the following lemmas required by analysis on convergence and optimality of the designed ROUTE-LRMR algorithm.

Lemma 2. Let {(U(t) ,V(t) ,L(t) ,W(t) )} be a sequence generated by Algorithm 1. Then the sequence approaches to a feasible solution.

Proof. First, we come to prove the boundedness of {Z(t) }. According to Theorem 1 and the optimality condition for (5) with respect to ˆL := UV, we have:

Through applying Lemma 3 on the above:

Lemma 3. [Fazel, 2002] Let H be a real Hilbert space endowed with an inner product 〈·, ·〉 and a corresponding norm ‖·‖, and any y ∈ ∂‖x‖, where ∂‖·‖ denotes the subgradient. Then ‖y‖∗ = 1 if x ̸= 0, and ‖y‖∗ ≤ 1 if x = 0, where ‖ · ‖∗ is the dual norm of the norm ‖ · ‖.

we obtain that the sequence {Z(t) } is bounded via observing the fact that the dual norm of ‖ · ‖∗ is the spectral norm. Together with the boundedness of {Z(t) } and limt→∞µ(t) = ∞, the relationship L(t) - U(t)V(t) = Z(t)-Z(t-1) µ(t-1) gives limt→∞L(t) - U(t)V(t) = 0. Further, the constraints of W+W = 1 and W,W ∈ [0, 1]m×n are immediately satisfied at each update, please see Lemma 1. Thus the statement holds.

Having the above theoretical results, we finally come to the convergence and optimality of ROUTE-LRMR.

Theorem3. The proposed Algorithm1 converges to a partial minimum to the optimization problem (5).

Proof. By Lemmas 1 and 2, Theorem 2, and the updating rule of multiplier Z, the KKT conditions for the constraints as well as the solutions to the variables U, V, L, W and Z are all satisfied. The first order optimality of (5) is sufficient to guarantee that the ROUTE-LRMR converges to a partial optimum to the problem (5).

Further, based on Theorem 1, the problem (5) is equivalent to the following one:

which is biconvex in W and L. The convergence to a partial optimum holds for the problem (13) too. Our ROUTE-LRMR is free to switch modes between NNM and BF. Concretely, instead of separately refreshing U and V, the updating of ˆL := UV in (7) can be done by minimizing the problem:

which can be solved in closed-form by the singular value thresholding [Cai et al., 2010]. Except for this step, no other changes happen in Algorithm 1.

4 Experimental Verification

In this section, we assess the performance of ROUTE-LRMR in comparison with several state-of-the-art methods including RegL1 [Zheng et al., 2012], PRMF [Wang et al., 2012], MoG [Meng and De la Torre, 2013], factEN [Kim et al., 2015], PSMSV [Oh et al., 2016] and Unifying [Cabral et al., 2013], the codes of which are either downloaded from the authors' websites or provided by the authors. Their settings follow the suggestions by the authors or the given parameters.

4.1 Synthetic Data

Data Preparation and Quantitative Metrics Similar to [Cand`es et al., 2011; Cabral et al., 2013], we generate a matrix Y0 as a product Y0 = U0V0. The U0 and V0 are of size m × r and r × n respectively, both of which are randomly produced by sampling each entry from the Gaussian distribution N(0, 1), leading to a ground truth rank-r matrix. Then we corrupt the entries via replacing a fraction s of Y0 with errors drawn from a uniform distribution over [-20, 20] at random, and the rest entries are polluted by Gaussian noise N (0, 0.12 ). To quantitatively measure the recovery performance, we employ 1) root mean square error (RMSE): 1 √mn ‖Y0 - ˆU ˆV‖F and 2) mean absolute error (MAE): 1 mn ‖Y0 - ˆUˆV‖1.

Parameter Effect We here focus on the parameter γ that controls the entropy term, the other two parameters α and β are empirically set to 50 and 1 throughout this paper. In this experiment, without loss of generality, square matrices of dimension m = n = 100 and rank r = 4 are considered. The left picture in Fig. 1 depicts RMSE and MAE curves (averaged over 10 trials) with respect to different outlier ratios. From the plots, we see that when γ approaches to 0, the errors rapidly go up. This is because, as analyzed in Sec. 3, the smaller γ is, the harder the weighting carries out, say the risk of being stuck into bad minima gets higher. It is also the evidence to prove the soft weighting is beneficial. In opposite, if γ gets too large, the performance also drops. The reason is that, in this situation, the weighting becomes almost constant (0.5 for each entry), which degenerates ROUTE-LRMR to PCA. Although the work range of γ shrinks as s grows, γ in [0.005, 0.8] can perform stably and sufficiently well. For the rest experiments unless stated otherwise, we set γ = 0.01. To better reveal the advantage of our method over the competitors especially on heavily ruined data, Table 1 reports the numerical comparison. As can be seen from Tab. 1, ROUTELRMR wins for all the cases, and the closest performance to ours is from Unifying. Note that the method HW is ROUTE with γ = 0.001 for mimicking the hard weighting strategy.

Convergence As regards convergence speed, for different cases, the right graph in Fig. 1 shows that the stop criterion quickly declines within 20 iterations, while the algorithm converges within 60 ∼ 80 iterations. Moreover, experimental findings here and follow-up tell that our algorithm has very stable convergence behavior even with respect to random initializations.

Tolerance to Outliers To more thoroughly show the tolerance to outliers, we fix m = n = 400 and test the tendency by varying outlier ratio s ∈ [0, 0.6] and rank r ∈ {20, 40, 60}. According to the results in Tab. 1, Unifying is the method chosen to compare. From the left picture of Fig. 2, we see that at the beginning, Unifying and our method are close in terms of RMSE, but as s increases, the margin between them enlarges. The second graph in Fig. 2 further confirms the first one. In the case of r = 20, both the RMSE and MAE of ROUTE-LRMF stay very low even when s reaches 0.6. The tolerance to outliers becomes weaker when r gets larger, not just for our method and Unifying but also for all the methods. The reason is that a higher-dimensional space requires more data to accomplish the recovery.

4.2 Real Data

Photometric Stereo Images of a static Lambertian object sensed by a fixed camera under a varying but distant point lighting source lie in a rank-3 subspace [Hayakawa, 1994]. This experiment aims to evaluate the effectiveness of the LRMR techniques on modeling the face under different illuminations. The cropped Extended YaleB-10 sequence, containing 64 faces of one subject with size 192 × 168, is adopted as the dataset. The light imbalance including shadows and highlights on the face significantly breaks the low-rank structure (please see the 1st column in Fig. 3 for example). In this part, we set the guess rank r to 5 for all the competitors.

Comparison Figure 3 gives several comparison. We can observe that PRMF, factEN, MoG and Unifying perform reasonably well, which are superior to PSMSV and L1Reg but inferior to ours. As shown in the 2nd and 4th rows of Fig. 3, PSMSV and L1Reg fail to remove shadows. The results by PRMF, factEN, MoG and Unifying, although recalling some details previously hidden in the dark, look unreal in the 2nd and 3rd cases. Our ROUTE-LRMR 1 provides visually pleasant and real results for all the given cases, the benefit of which mainly comes from the effective outlier detection. The 2nd column in Fig. 3 displays the estimated weights W (brighter regions indicate closer values to 1, while darker ones stand for those to 0), from which we can find our strategy successfully detects and thus eliminates outliers. On the right side of Fig. 3, we further provide several results by our method (only, due to space limit). One may wonder if the weights can be formed by treating as outliers the pixels with intensity greater (highlights) or lower (shadows) than predefined thresholds like [Zheng et al., 2012]. This way can reduce the problem to LRMC, but is too heuristic, at high risk of sacrificing much useful information for recovery. Taking the bottom-right original for example, the thresholding may determine all the pixels as outliers, while our strategy can finish the job wisely and nicely. Moreover, in many real-world applications, manually seeking appropriate thresholds is, if not impossible, very difficult. Being able to adaptively assign weights to data is definitely desired, which is the goal and motivation of our design.

1 In image/video data, the outliers, such as shadows and foregrounds, often appear coherently. Considering this, in this experiment, we employ a 2 × 2 median filter on W.

5 Conclusion

This paper has shown a method for jointly detecting outliers and recovering the underlying low-rank matrix, called ROUTE-LRMR. Our weighting strategy employs an entropy regularization term to minimize the prediction bias, which behaves like a sigmoid function. To seek the optimal solution for ROUTE-LRMR, we have developed an Alternating Direction Minimization based algorithm. The theoretical analysis and the experimental results compared to the state-of-the-arts, have demonstrated the advantages of the proposed ROUTELRMR. Our strategy can be applied to numerous tasks such as regression, clustering, inpainting and foreground detection. It is also ready to embrace specific domain knowledge, like graph regularizer on the weight, for further boosting the performance on different applications.

Acknowledgments

X. Guo is supported by National Natural Science Foundation of China (grant no. 61402467). Z. Lin is supported by National Basic Research Program of China (973 Program) (grant no. 2015CB352502) and National Natural Science Foundation of China (grant nos. 61625301 and 61231002).