(PCA) problem that has been extensively studied in the literature [5]- [7].

A critical assumption of RMC is that the underlying matrix is of low rank. Therefore, RMC can be naturally formulated as a rank constrained optimization problem. However, the rank function is usually difficult to handle because of its combinatorial nature. To resolve this issue, there is a line of research focusing on convex relaxations of RMC formulations [5]- [7]. For example, a widely used technique in these works is to replace the rank function by the nuclear norm of the matrix under consideration. However, algorithms for solving convex relaxations involving nuclear norm usually require to compute a singular value decomposition (SVD) in each iteration, which can be very time consuming for large-scale problems. Recently, nonconvex RMC formulations based on low-rank matrix factorization were proposed in the literature [8]- [11]. In these formulations, the target low-rank matrix is factorized as the product of two much smaller matrices so that the low-rank property is automatically satisfied. The idea of factorization greatly reduces the price of promoting the low-rankness. Another critical assumption of RMC is that the matrix corresponding to the noise is sparse. Usually the 1 norm of the noise matrix is penalized to promote its sparsity.

In this paper, motivated by the recent developments on Riemannian optimization, we propose a new nonconvex formulation for RMC. Specifically, our new formulation is a nonsmooth optimization problem over the Grassmann manifold. The idea of modeling matrix completion as Riemannian optimization has been explored in the literature [12]- [16]. Modeling RMC as a Riemannian optimization is more challenging, because of the existence of the function promoting sparsity of the noise. Recently, Zhang and Yang [17] proposed to model RMC as a Riemannian optimization over the fixed-rank manifold. This formulation requires a priori knowledge of the level of sparsity of the sparse noise matrix, which might be hard to estimate in practice. Cambier and Absil [18] considered minimizing a nonsmooth objective over the fixed-rank manifold. However, their algorithm needs to smooth the objective using certain smoothing technique, and thus does not solve the original nonsmooth problem. The GRASTA algorithm [19], [20] proposed by He et al. minimizes a nonsmooth objective over the Grassmann manifold. The authors suggested to use the alternating direction method of multipliers (ADMM) to solve it, but it is known that ADMM lacks convergence guarantees in this case.

We propose an alternating manifold proximal gradient (AManPG) algorithm for solving our new RMC model. The AManPG algorithm alternatingly updates the low-rank matrix using a Riemannian gradient step and the sparse matrix using a proximal gradient step. This method is motivated by the ManPG algorithm recently proposed by Chen et al. [21]. We show that the iteration complexity of AManPG is O( -2 ) for obtaining an -stationary point. Note that our AManPG solves the nonsmooth problem directly, and it does not need any smoothing technique as used in [18]. Moreover, we adopt a continuation technique [22] to further improve the computational efficiency of AManPG. Our numerical experiments show that the proposed AManPG with continuation compares favorably to some existing methods for solving the RMC problem.

The rest of this paper is organized as follows. In Section II, we briefly review some existing related work on the RMC problem and give our new RMC formulation. In Section III, we review some basics of manifold optimization and propose the ManPG algorithm for our RMC formulation. In Section IV and Section V, we propose our AManPG with continuation algorithm and provide its convergence analysis. In Section VI, numerical results are presented to demonstrate the advantages of the proposed new formulation and algorithm. Finally, we draw some conclusions in Section VII.

In this section, we discuss some related work in the literature, and then present our RMC formulation.

Robust PCA [5], [6] is an important tool in data analysis and has found many interesting applications in computer vision, signal processing, machine learning, and statistics etc. The goal is to decompose a given matrix M ∈ R m×n into the superposition of a low-rank matrix L and a sparse matrix S, i.e., M = L + S. The works by Candès et al. [5] and Chandrasekaran et al. [6] formulate the problem as the following convex optimization problem:

where γ > 0 is a weighting parameter, the nuclear norm L * sums the singular values of L, and the 1 norm S 1 sums the absolute values of all entries of S, i.e.,

When only a subset of the entries of M is observed, robust PCA becomes the robust low-rank matrix completion problem [18]. Similar to (1), a convex formulation for RMC can be cast as follows:

min

where Ω denotes the set of the indices of the observed entries and P Ω : R m×n → R m×n denotes a projection defined as

The convex formulations (1) and ( 2) have been studied extensively in the literature, and we refer to the recent survey paper [7] for algorithms for solving them.

By assuming that the rank of L is known (denoted by r), the idea of many nonconvex formulations is based on the fact that L can be factorized to L = UV , where U ∈ R m×r , V ∈ R r×n . In this paper, we also utilize this idea and propose to solve the following formulation of RMC:

where 0 < λ < 1, γ > 0, and U denotes an orthonormal basis of the subspace U. We show that the manifold proximal gradient method (ManPG) proposed by Chen et al. [21] can be applied to solve (3) and the corresponding convergence analysis applies naturally. We then propose a variant of ManPG, named alternating ManPG (AManPG) that can significantly improve the efficiency of ManPG for solving (3). We further rigorously analyze the convergence rate of AManPG. Compared with GRASTA [19], [20], our proposed algorithms for solving (3) have rigorous convergence guarantees and convergence rate analysis. Compared with (9) and the algorithm proposed by Cambier and Absil [18], our algorithms solve the nonsmooth problem (3) directly without any smoothing technique. Finally, to further accelerate the convergence of AManPG, we incorporate the so-called continuation technique on the weighting parameter γ in (3). Our numerical results on both synthetic data and real data on background extraction from surveillance video demonstrate that our AManPG with Continuation (AManPGC) algorithm, compares favorably with existing methods for RMC.

Replacing L by UV in (1) and (2) leads to various nonconvex formulations for robust PCA and RMC. In particular, Li et al. [23] suggested using subgradient method to solve the following nonsmooth robust matrix recovery model

where vector y denotes the linear measurements to L via the known linear operator A : R m×n → R |Ω| . Shen et al. [11] proposed the LMaFit algorithm that implements an ADMM for solving the following nonconvex formulation of robust PCA:

Note that if ( Û, V ) solves (5), then ( ÛQ, Q -1 V ) also solves (5) for any invertible Q ∈ R r×r . Since all matrices ÛQ share the same column space, Dai et al. [12], [13] exploited this fact and formulated the matrix completion problem as the following optimization problem over a Grassmann manifold:

Authorized licensed use limited to: Univ of Calif Davis. Downloaded on May 20,2021 at 00:39:57 UTC from IEEE Xplore. Restrictions apply.

where Gr(m, r) denotes the Grassmann manifold, i.e., the set of r-dimensional vector subspaces of R m , and U ∈ R m×r is any matrix whose column space is U and is often chosen to be orthonormal. However, it is noticed that the outer problem might be discontinuous at points U for which the V problem does not have a unique solution. To address this issue, Keshavan et al. [14], [15] proposed to optimize both the column space and row space at the same time, which results in the following so-called OptSpace formulation for matrix completion: min U∈Gr(m,r),V∈Gr(n,r)

where U and V are any orthonormal bases of U and V, λ > 0 is a weighting parameter, and the regularizer U ΣV 2 F is used so that the outer problem is continuous. Boumal and Absil [16], [24] proposed to study the following variant of (6):

where Ω is the complement of Ω. They proposed to use the Riemannian trust-region method to solve this problem. Comparing with OptSpace (7), formulation (8) has a much smaller search space. Note that in (8), λ is usually chosen to be very close to zero, as it indicates that we have a small confidence that the entries (UV ) ij for (i, j) / ∈ Ω are equal to zero. For RMC, Cambier and Absil [18] proposed the following Riemannian optimization formulation:

where M r denotes the fixed-rank manifold, i.e., M r := {X | rank(X) = r}. The algorithm proposed in [18] needs to smooth the 1 norm first to change the problem to a smooth problem, and then apply the Riemannian conjugate gradient method to solve the smoothed problem. As a result, the algorithm in [18] does not solve (9) exactly but only solves an approximation of it. Related to (6) and ( 5), He et al. proposed the GRASTA algorithm [19], [20] that can be used to solve the following formulation of RMC:

where U is again an orthonormal basis of U. GRASTA uses alternating minimization and ADMM to solve (10), which is efficient in practice but lacks convergence guarantees. Recently, Zhang and Yang [17] studied the following formulation of robust PCA over the fixed-rank manifold:

where g : R m×n → R m×n is a hard thresholding procedure that depends on a parameter specifying the level of sparsity of matrix L -M . This parameter might be difficult to estimate in practice.

A similar model for RMC was also considered in [17].

In the following we list a few other works related to robust PCA and RMC. Zhao et al. [25] proposed two smooth functions that can promote the robustness against two types of outliers. He et al. [26] derived a correntropy-based cost function and applied the half-quadratic technique to solve the matrix completion problem. Zeng and So [27] proposed the iterative p -regression algorithm and ADMM for a RMC model that uses p norm to promote the sparisty. Zhang and Wang [28] considered the low-rank Hankel matrix recovery problem with sparse noise, and proposed an alternating projection type method to solve it. Dutta et al. [29] formulated robust PCA as a nonconvex feasibility problem and proposed an alternating projection method to solve it. In another work, Dutta et al. [30] proposed to formulate the robust PCA as the best pair problem that aims to find a pair of points minimizing the distance between two disjoint sets. An accelerated proximal gradient method was designed in [30] to solve this problem. There is a vast literature on robust PCA and RMC and the ones that we discussed above are by no means exhaustive. We refer the reader to [7], [31]- [33] for more comprehensive surveys on these topics.

In this section, we show that the ManPG algorithm recently proposed by Chen et al. [21] can be naturally adopted to solve (3). Note that the ManPG algorithm was originally proposed for solving problems over the Stiefel manifold, but the manifold in ( 3) is Grassmann manifold. Therefore, we need to further elaborate on the details how ManPG works on Grassmann manifold. To this end, we first introduce some backgrounds on the geometry of Grassmann manifold.

In this subsection we briefly introduce concepts and properties of the Grassmann manifold. Most of the materials here are from [16], and we include them here for the ease of discussion later. The Grassmann manifold Gr(m, r) is the set of r-dimensional linear subspaces of R m endowed with quotient manifold structure, whose dimension is dim(Gr(m, r)) = r(mr) [34]. Each point U ∈ Gr(m, r) is a linear subspace spanned by the column space of a full-column-rank matrix U :

where R m×r * denotes the set of all m × r matrices with full column rank, and col(U ) denotes the subspace spanned by the columns of U . For numerical reasons, we always represent U by its orthonormal basis, i.e., U ∈ U(m, r) := {U ∈ R m×r : U U = I r }, where U (m, r) denotes the Stiefel manifold. From now on, by slightly abusing the notation, we use the matrix U ∈ U(m, r) to represent U ∈ Gr(m, r). Since multiplying by an r × r orthonormal matrix does not change the column space of U , we can regard Gr(m, r) as a quotient of R m×r * by the equivalent relation U = UQ, where Q is any r × r orthonormal matrices. Throughout this paper, we use the Riemannian metric that is induced from the Euclidean inner product, i.e., U, V = Tr(U V ). Endowed with this Riemannian metric, the Grassmann manifold is also a Riemannian quotient manifold, and it admits a tangent space to Gr(m, r) at each point U:

Note that here we again slightly abuse the notation. More specifically, we use the matrix H ∈ R m×r to represent the tangent vector H ∈ T U Gr(m, r). For the Riemannian manifold M, the Riemannian gradient of a smooth function f : M → R is defined as follows.

Definition 1: (Riemannian Gradient) Given a smooth function f : M → R, the Riemannian gradient of f at X ∈ M, denoted by gradf (X), is the unique tangent vector in

where Df denotes the directional derivative of f . For the Grassmann manifold, the orthogonal projector from R m×r onto the tangent space T U Gr(m, r) is given by (here we again slightly abuse the notation and use the orthonormal basis U to represent the subspace U):

The definition of retraction operation for manifold M is given below. Definition 2: (Retraction) Let Retr X (ξ) : TM → M be a mapping from the tangent bundle TM to the manifold M. We call Retr X (•) a retraction at X if the following hold for any X ∈ M and ξ ∈ T X M:

In our numerical experiments, we choose the QR decomposition as the retraction for Grassmann manifold:

where qf(X) denotes the Q-factor of the QR decomposition of X. Here we again use U to represent the subspace U and use H to represent the element in the tangent space. Note that because of ( 12), we know that U + H is always full-rank, and therefore the output of the QR decomposition is unique.

Recently, Chen et al. [21] proposed a novel ManPG algorithm for solving nonsmooth optimization problem over the Stiefel manifold St(n, r) in the following form:

in which F 1 is smooth with Lipschitz continuous gradient, F 2 is nonsmooth and convex. Here the smoothness, convexity and Lipschitz continuity are interpreted when the functions are considered in the ambient Euclidean space. A typical iteration of ManPG algorithm for solving (17) is as follows:

where t > 0 and α k > 0 are step sizes, and ∇F 1 denotes the Euclidean gradient of F 1 . Chen et al. [21] proved that the iteration complexity of ManPG (18) is O(1/ 2 ) for obtaining an -stationary solution of (18). They also demonstrated that ManPG is very efficient for solving sparse PCA and compressed modes problems. Now we discuss how to apply ManPG [21] to solve (3). For ease of presentation, we denote the smooth part of F in (3) by f , i.e.,

) Note that for fixed U and S, the optimal V of F (U, V, S) (and f (U, V, S)) is uniquely determined. Therefore, by denoting

and

we know that our RMC formulation (3) reduces to

where U is an orthonormal basis of U. Note that although ManPG was originally designed for problems over the Stiefel manifold, we find that it can also be applied to problems over the Grassmann manifold. Also note that as long as the tangent space of the manifold can be represented by a linear operator, the subproblem of ManPG can be solved efficiently using a semi-smooth Newton method (see Section 4.2 of [21]), and this happens to be the case for the Grassmann manifold (e.g., (12)). Therefore, ManPG can be applied to solve (22) and the subproblems can be solved relatively easily. It is easy to see that ManPG for solving (22) reduces to the following updates in the k-th iteration:

where t S , t U , α and β are all step sizes, and h(S) := γ P Ω (S) 1 . We now make some necessary remarks on this ManPG algorithm (23). First, ( 23) is actually slightly different with a direct application of ManPG for solving (22). For a direct application of ManPG, we have t S = t U and α = β. Here in (23) we allow these step sizes to be different so that we have more freedom to choose the best step sizes in practice. Also we note that (23c) and (23d) correspond to a Riemannian gradient step with respect to the U variable. The updates (23a) and (23b) correspond to a proximal gradient step for the S variable in the Euclidean space. This can be interpreted in the following way. First, in RMC (22), the U variable does not appear in the nonsmooth part of the objective, so it is natural to perform a Riemannian gradient step for U . Second, for fixed U , the S Algorithm 1: ManPG for Solving RMC (22) problem is only an unconstrained problem in the Euclidean space, so we take a proximal gradient step. Moreover, the two subproblems (23a) and (23c) are both easy to solve. Specifically, (23c) can be reduced to

i.e., it is the negative Riemannian gradient of f multiplied by the step size t U . The ΔS subproblem (23a) can be solved by a simple 1 norm shrinkage operation (note that we are only interested in the P Ω (ΔS k ) and PΩ(ΔS k ) can be simply set to 0):

Now, to implement the ManPG (23), the only remaining component is to calculate the Riemannian gradient grad U f (U, S) used in (24). The procedure for computing it is outlined in [24]. By assuming that the subspace of the Grassmann manifold is represented by orthonormal bases, which means U is restricted to the Stiefel manifold, the Riemannian gradient of the smooth function f (U, S) with respect to U is given by:

where C ∈ R m×n is the mask operator whose components are given by: C ij = 1 if (i, j) ∈ Ω, and C ij = 0 otherwise. Here V U,S is defined in (20), and it can be computed as follows:

where vec denotes the vectorization operator, A is defined as

and ⊗ denotes the Kronecker product. For more details about these calculation, we refer the reader to [24].

With these preparations, we can finally summarize the ManPG algorithm (23) for solving (3) (or, (22)) as in Algorithm 1.

Remark

It should be noted that ManPG updates S and U in parallel. That is, ManPG ( 23) is a Jacobi type iterative algorithm and for fixed (U k , S k ), (23a) and (23c) can be computed in parallel. One way that can possibly improve the speed of ManPG is to use a Gauss-Seidel type algorithm. This idea has also been adopted in [35], where the authors showed that the Gauss-Seidel type ManPG performs much better than the original Jacobi type ManPG. Motivated by this, here we also propose an alternating ManPG (AManPG), which updates S and U sequentially, instead of in parallel. However, one crucial thing to note here is that the V variable will need to be re-calculated, when we have a new S variable before we update the U variable. Our AManPG algorithm is summarized in Algorithm 2.

Remark 4: When we compute ΔU k in AManPG, we used the latest S k+1 , which requires us to compute the latest V k U k ,S k+1 . While in ManPG, we used S k in the updates of ΔU k , and this does not require us to compute another V k . This is the main difference between AManPG (Algorithm 2) and ManPG (Algorithm 1). In both algorithms, we always set α = β = 1. Noticing that the matrix A only depends on variable U , and there is no need to recalculate A when computing V k U k ,S k+1 . The continuation technique: There are two parameters in the model (3): λ and γ. Since λ indicates our confidence level of the entries of (UV ) being zero, it needs to be very small. In our numerical experiments, we choose λ = 10 -8 . The parameter γ in (3) controls the sparsity level of P Ω (S). A larger γ yields sparser P Ω (S). In practice, we usually do not know how sparse the matrix S should be, and it is therefore not easy to determine γ. A usual practice in the literature to deal with this issue is to conduct a continuation technique on γ. Roughly speaking, Algorithm 3: AManPG With Continuation (AManPGC) for Solving (22).

1: Input: Step sizes t S , t U , parameters γ 0 γ min , shrinking factors μ 1 < 1, μ 2 < 1, initial accuracy tolerance 0 2: Initialize: U 0 , S 0 . Set = 0 3: while γ > γ min do 4:

Call AManPG to solve (22) with γ = γ , and set the output of AManPG as (U +1 , S +1 ).

γ +1 = μ 1 γ 6: +1 = μ 2 7:

= + 1 8: end while 9: Output: U , S the continuation starts with solving (3) using a relatively large γ. Then the parameter γ is decreased and ( 3) is solved again. This process is repeated until γ is very small. This idea has been widely adopted in the literature, e.g., [22], [36]- [38]. Combining this continuation idea with our AManPG algorithm, we obtain the AManPGC algorithm which works greatly in practice as confirmed by our numerical results in Section VI. AManPGC is summarized in Algorithm 3. Note that in Algorithm 3, we also shrink the accuracy tolerance in each iteration, as we want to solve the problem more and more accurately.

In this section, we analyze the convergence behavior and iteration complexity of AManPG (Algorithm 2). To simplify the notation, we denote M = Gr(m, r) and h(S) = γ P Ω (S) 1 . We analyze the convergence of AManPG for solving the following problem:

where f (U, S) is smooth and h(S) is nonsmooth and convex.

Here the smoothness and convexity are interpreted when the functions are considered in the ambient Euclidean space. For simplicity, we rewrite the AManPG for solving (28) here. One typical iteration of AManPG for solving (28) is the same as (23) except that (23c) is replaced by

Note that we use S k+1 in (29) to replace S k in (23c).

We make the following assumptions to (28) throughout this section.

Assumption 5: (F is lower bounded) There exists a finite constant F * , such that

Note that for (22), it is easy to see that F * = 0.

The following assumption is about grad U f (U, S), which regards the regularity of the pullback function f (ΔU, S) = f (Retr U (ΔU ), S), and differs from the standard Lipschitz continuity assumption because of the retraction operator. This assumption was originally suggested in [39].

Assumption 7: (Restricted Lipschitz-type gradient for pullbacks) There exists L U ≥ 0 such that, for sequence (U k , S k ) k≥0 generated by AManPG (Algorithm 2), the pullback function

Note that since the Grassmann manifold is a compact submanifold of R m , for fixed S, by Lemma 2.7 of [39], Assumption 7 holds if f (U, S) has locally Lipschitz continuous gradient in the Euclidean space. Assumption 7 allows simple extensions of existing proofs in the Euclidean case. One limitation of this assumption is that it might be difficult to estimate L U in practice.

From the Theorem 4.1 in [40], we can define the stationary point of problem (28) as follows.

Definition 8: (Stationary point). A pair of (U, S) ∈ M × R m×n is called a stationary point of problem (28) if it satisfies the first-order necessary conditions:

According to Theorem 4.1 in [40], the optimality conditions of the subproblems ( 29) and (23a) are

If ΔU k = 0 and ΔS k = 0, then we know that S k+1 = S k , and (U k , S k ) satisfies ( 30) and thus is a stationary point of (28). Therefore, we can use the norm of (U k , S k ) to measure the closeness to stationary point, and we define the -stationary point of (28) as follows. Definition 9: ( -stationary point). We say that

given by (23a) and ( 29) satisfies

where L := min{L S , L U }. Now we are ready to analyze the iteration complexity of AManPG for obtaining an -stationary point of (28). The convergence analysis presented here mainly follows from [21], [35]. The main difference lies in the Assumption 7, which was not used in the previous work [21], [35]. This assumption helps us in simplifying and extending the analysis from [21], [35] to the Grassmann manifold. First, we present two lemmas, which show that there is a sufficient reduction of the objective value after each update of S and U . These two lemmas essentially follow Lemma 4 of [41].

Lemma 10: Assume Assumption 6 holds, and the sequence (S k , U k , ΔS k , ΔU k ) is generated by AManPG. By choosing t S = 1/L S and α = 1, the following inequality holds

Proof: See Appendix A for details. Lemma 11: Assume Assumption 7 holds, and the sequence (S k , U k , ΔS k , ΔU k ) is generated by AManPG. By choosing t U = 1/L U and β = 1, the following inequality holds

) Proof: See Appendix B for details. Now we are ready to present our main convergence result of AManPG, which essentially follows Theorem 1 of [35].

Theorem 12: Assume Assumptions 5, 6 and 7 hold. By choosing t U = 1/L U , t S = 1/L S , α = β = 1 in AManPG (Algorithm 2), every limit point of the sequence {U k , S k } generated by AManPG is a stationary point of problem (28). Moreover, AManPG returns an -stationary point of problem (28) in at most 2L(F (U 0 , S 0 ) -F * )/ 2 iterations, where L := min(L S , L U ).

Proof: See Appendix C for details.

In this section, we provide numerical results for both synthetic and real datasets to illustrate the performance of the proposed algorithms. We focus on comparing our ManPG and AManPG algorithms with some baseline algorithms for solving nonconvex RMC formulations, in particular, the subgradient method (SubGM) [23], the LMaFit method [11] and the Riemannian conjugate gradient method for the smoothed 1 -norm objective function (RMC) [18]. We use the same continuation framework for ManPG and call it ManPGC. For the SubGM method, we assume that the linear operator A is a simple projection. For the LMaFit algorithm, we turn off the rank estimation since we assume that the rank is known for all cases. We use the original setting for the RMC algorithm. All algorithms use the same C-Mex code for accelerating the matrix multiplication between a sparse matrix and a full matrix and some other bottleneck computations. Moreover, to guarantee a fair comparison, each algorithm is carefully tuned to achieve its best performance. All experiments were run on Matlab R2018b with a 2.3 GHz Dual-Core Intel Core i5 CPU.

We first test our ManPGC and AManPGC algorithms on different synthetic data. After picking the values of m, n, r, we generate U * ∈ R m×r , V * ∈ R r×n with i.i.d. normal entries of zero mean and unit variance. The target low-rank matrix is X * = U * V * . We test two different ways for generating the sparse matrix S * in which the positions of the nonzero entries are chosen uniformly random: (i) Gaussian (Cases 1-6 in Table I): the nonzero entries of S * follow the normal distribution N (0, 1); I): the nonzero entries of S * are uniform in [-500, 500]. Finally we set M = X * + S * and P Ω (M ) is the observed matrix, where Ω denotes the indices set of nonzero entries of S * . We set = 10 -8 L for the stopping criteria in Algorithms 1 and 2.

We test our algorithms for the cases in Table I, where sr denotes the sampling ratio and sp denotes the sparsity level, i.e., the ratio of nonzero entries.

Parameters: For Cases 1-6, we choose t S = 1, and for Cases 7-12, we choose t S = 0.8. We choose t U = 2/|Ω| for both AManPGC and ManPGC in all cases, except Case 4 in which we choose t U = 1.6/|Ω|, 0 = 20 for ManPGC. Moreover, we set γ 0 = 10, λ = 10 -8 , μ 1 = μ 2 = 1/10, 0 = 30 for Cases 1-6. For Cases 7-8, we set μ 1 = 1/10, μ 2 = 1/20, 0 = 200 for AManPG and μ 1 = 1/3, μ 2 = 1/8, 0 = 400 for ManPG. For Cases 9-12, we set μ 1 = 1/10, μ 2 = 1/20, 0 = 200 for AManPG and μ 1 = 1/6, μ 2 = 1/16, 0 = 200 for ManPG. We use the SVD of the observed matrix P Ω (M ) as the initial point for all algorithms. Note that the step sizes t S and t U can also be estimated using the line search technique, which should be used in practice. Here we chose to tune constant step sizes because we observed that they work better than the ones found by line search on the tested problems.

In the k-th iteration, we calculate the relative difference for each method as

and report the Relative Difference versus the CPU time (in seconds) in Fig. 1.

From Fig. 1 we can see that our AManPGC algorithm usually outperforms other baseline algorithms on CPU time. From Fig. 1, we see that AManPGC is the fastest algorithm in almost all 12 tested cases. Moreover, ManPGC also works well for Cases 1-6, while shows some disadvantage comparing with RMC for Cases 7-12. However, ManPGC is still better than SubGM and LMaFit in all cases. Here we discuss some possible reasons why ManPGC and AManPGC can be more efficient than RMC [18], SubGM [23] and LMaFit [11]. For RMC, it needs to solve a sequence of smooth Riemannian optimization problems using the Riemannian conjugate gradient method. As a contrast, our ManPGC and AManPGC only solve one (though nonsmooth) Riemannian optimization problem. Note that SubGM uses subgradient information, and requires a geometrically decreasing step size to guarantee the convergence. Therefore, it can be slow in practice. As a contrast, our ManPGC and AManPGC do not use the subgradient information, and allow constant step size, and therefore can be more efficient and achieve higher accuracy in practice. For LMaFit, note that it is essentially a nonconvex ADMM algorithm. Its convergence speed depends on the penalty parameter used in the augmented Lagrangian function, which is difficult to tune in practice. It is also worth pointing out that there still lacks a convincing convergence guarantee for LMaFit thus far, as the one given in [11] requires a strong and unverifiable condition. Therefore, it is not clear under what conditions LMaFit converges. Fig. 2. Recovery of full matrices from their entries. We declare X to be recovered if the relative difference satisfies UV -X * F / X * F < 10 -2 . Top row: phase transition for synthetic datasets by varying rank and outlier sparsity. We fix the sampling ratio = 0.1. Botton row: phase transition for synthetic datasets by varying rank and the sampling ratio. We fix the outlier sparsity = 0.1. For each rank-sampling ratio or rank-sparsity pair, we run three algorithms for 20 times and report the empirical recovery rate. White denotes perfect recovery in all experiments, and deep blue denotes failure for all experiments. We further plot the phase transition for the proposed AManPGC, ManPGC algorithms and the RMC formulation, which performs the best among all baseline methods. We set m = n = 1000, t S = 1, t U = 2/|Ω|, 0 = γ 0 = 1, μ 1 = μ 2 = 1/5, for both AManPGC and ManPGC. We set the maximum iteration number for all three algorithms to be 100. Fig. 2 shows the recovery rate for different ranks, outlier sparsities and sampling ratio. We see from Fig. 2 that the AManPGC and the RMC algorithms perform similarly and they are both better than the ManPGC algorithm.

We now evaluate the performance of ManPGC and AManPGC for video background estimation using three data sets suggested in [42]. Each data set contains a sequence of grayscale frames each of which can be regarded a matrix. By stacking the columns of each frame of the video into a long vector, we obtain a matrix whose columns correspond to the frames. This matrix can be decomposed to the matrix corresponding to the static background plus the matrix corresponding to the moving foreground. The former one is expected to have low rank, because ideally the columns corresponding to the background are all the same; and the latter one is expected to be sparse, because the moving foreground usually only occupies a small portion of the frame. In our experiment, we observe that we are able to extract the background of the video using only partial observation. We apply our algorithms for this task and compare with existing methods. The three data sets are: We set r = 1 and test two cases for the ratio of observed entries: 50% and 10%. We terminate the algorithms when the following inequality holds:

which indicates that the algorithm does not make much progress. We plot the recovered background frames in Figs. 3, 4 and5. The CPU time (in and the number of iterations are reported in Tables II, III, and IV. From Tables II, III, and IV, we see that ManPGC and AManPGC are faster than the other three algorithms in all cases. Moreover, we see that even with

where X * denotes the ground truth, and

The d metric of the recovered backgrounds for different are shown in Fig. 6. From Fig. 6 we see that all the tested algorithms yield good results for recovering the backgrounds and the foregrounds. More specifically, from Fig. 6 we see when 50% entries of M are observed, SubGM is slightly worse than others for the "Escalator" data, and when 10% entries of M are observed, LMaFit is slightly worse than others.

In this paper, we have proposed a new formulation for RMC as a Riemannian optimization over the Grassmann Manifold. Inspired by the recent work ManPG, we have developed a new algorithm called AManPG for solving this nonconvex and nonsmooth manifold optimization problem. The iteration complexity of AManPG for obtaining an -stationary solution is rigorously analyzed, following a similar way as ManPG. In our numerical experiments, we have tested our proposed formulation and algorithms for both synthetic and real datasets, and the numerical results show that our proposed formulation and methods usually outperform some popular existing methods for RMC.

Proof: For fixed U ∈ M and S ∈ R m×n , define g U,S (T ) := ∇ S f (U, S), T + 1 2t S T 2 F + h(S + T ).

Note that g U,S is (1/t S )-strongly convex, therefore it holds that g U,S (T 1 ) ≥ g U,S (T 2 ) + ∂g U,S (T 2 ), T 1 -T 2

By letting T 1 = 0, T 2 = ΔS k in (39), we have

where the equality is from the optimality condition of (23a), i.e., 0 ∈ ∂g U k ,S k (ΔS k ). Using Assumption 6, we can get

, where the first inequality comes from (41), and the second inequality is from (40). This completes the proof.

Proof: From Assumption 7, we have

where the second inequality is due to (31). This completes the proof.

Proof: Combining (33) and (34) yields

Since F is decreasing and bounded below, we have

It follows that every limit point of {(U k , S k )} is a stationary point of (28). Moreover, if AManPG does not terminate after K