The problem of finding a subspace that captures the features of a given dataset and possesses certain properties is at the heart of many machine learning applications. One commonly encountered formulation of the problem, which is motivated largely by sparsity or robustness considerations, is given by
where Y ∈ R n×p is a given matrix and • r denotes the r norm of a vector. To better understand how problem (1) arises in applications, let us consider two representative examples.
• Orthogonal Dictionary Learning (ODL). The goal of ODL is to find an orthonormal basis that can compactly represent a given set of p (p n) data points y 1 , . . . , y p ∈ R n . Such a problem arises in many signal and image processing applications; see, e.g., [1], [2] and the references therein. By letting Y = [y 1 , . . . , y p ] ∈ R n×p , the problem can be understood as finding an orthogonal matrix X ∈ R n×n and a sparse matrix A ∈ R n×p such that Y ≈ XA. Noting that this means X Y ≈ A should be sparse, one approach is to find a collection of sparse vectors from the row space of Y and apply some post-processing procedure to the collection to form the orthogonal matrix X. This has been pursued in various works; see, e.g., [3]- [6]. In particular, the work [6] considers the formulation (1) and shows that under a standard generative model of the data, one can recover X from certain local minimizers of problem (1).
• Robust Subspace Recovery (RSR). RSR is a fundamental problem in machine learning and data mining [7]. It is concerned with fitting a linear subspace to a dataset corrupted by outliers. Specifically, given a dataset Y = [X, O]Γ ∈ R n×(p1+p2) , where the columns of X ∈ R n×p1 are the inlier points spanning a d-dimensional subspace S of R n (d < p 1 ), the columns of O ∈ R n×p2 are outlier points without a linear structure, and Γ ∈ R (p1+p2)×(p1+p2) is an unknown permutation, the goal is to recover the inlier subspace S, or equivalently, to cluster the points into inliers and outliers. One recently proposed approach for solving this problem is the so-called dual principal component pursuit (DPCP) [8], [9]. A key task in DPCP is to find a hyperplane that contains all the inliers. Such a task can be tackled by solving problem (1). In fact, it has been shown in [8], [9] that under certain conditions on the inliers and outliers, any global minimizer of problem (1) yields a vector that is orthogonal to the inlier subspace S.
Despite its attractive theoretical properties in various applications, problem (1) is numerically challenging to solve due to its nonsmooth objective and nonconvex constraint. Nevertheless, the manifold structure of the constraint set M := {x ∈ R n | x 2 = 1} suggests that problem (1) could be amenable to manifold optimization techniques [10]. One approach is to apply smoothing to the nonsmooth objective in (1) and use existing algorithms for Riemannian smooth optimization to solve the resulting problem. For instance, when tackling the ODL problem, Sun et al. [5], [11] and Gilboa et al. [12] proposed to replace the absolute value function t → |t| by the smooth surrogate t → h µ (t) = µ log(cosh(t/µ)) with µ > 0 being a smoothing parameter, while Qu et al. [13] proposed to replace the 1 norm with the 4 norm. They then solve the resulting smoothed problems by either the Riemannian trust-region method [5], [11] or the Riemannian gradient descent method [12], [13]. Although it can be shown that these methods will yield the desired orthonormal basis under a standard generative model of the data, the smoothing approach can introduce significant analytic and computational difficulties [6]. Another approach, which avoids smoothing the objective, is to solve (1) directly using Riemannian nonsmooth optimization techniques. For instance, in the recent work [6], Bai et al. proposed to solve (1) using the Riemannian subgradient method (RSGM), which generates the iterates via
Here, η k > 0 is the step size; ∂ R f (•) denotes the Riemannian subdifferential of f and is given by
where I n is the n × n identity matrix and ∂f (•) is the usual subdifferential of the convex function f [14, Section 5]. Bai et al. [6] showed that for the ODL problem, RSGM with a suitable initialization will converge at a sublinear rate to a basis vector with high probability under a standard generative model of the data. Moreover, by running RSGM O(n log n) times, each time with an independent random initialization, one can recover the entire orthonormal basis with high probability. Around the same time, Zhu et al. [9] proposed a projected subgradient method (PSGM) for solving (1). The method generates the iterates via
The updates (2) and (3) differ in the choice of the direction v k -the former uses a Riemannian subgradient of f at x k , while the latter uses a usual Euclidean subgradient. For the DPCP formulation of the RSR problem, Zhu et al. [9] showed that under certain assumptions on the data, PSGM with suitable initialization and piecewise geometrically diminishing step sizes will converge at a linear rate to a vector that is orthogonal to the inlier subspace S. The step sizes take the form η k = η (k-K0)/K +1 , where η ∈ (0, 1) and K 0 , K ≥ 1 satisfy certain conditions. In practice, however, the parameters η, K 0 , K are difficult to determine. Therefore, Zhu et al. [9] also proposed a PSGM with modified backtracking line search (PSGM-MBLS), which works well in practice but has no convergence guarantee.
Although the results in [6], [9] demonstrate, both theoretically and computationally, the efficacy of RSGM and PSGM for solving instances of (1) that arise from the ODL and RSR problems, respectively, two fundamental questions remain. First, while PSGM can be shown to achieve a linear convergence rate on the DPCP formulation of the RSR problem [9], only a sublinear convergence rate has been established for RSGM on the ODL problem [6]. Given the similarity of the updates (2) and (3), it is natural to ask whether the slower convergence rate of RSGM is an artifact of the analysis or due to the inherent structure of the ODL problem. Second, the convergence analyses in [6], [9] focus only on the ODL and RSR problems. In particular, they do not shed light on the performance of RSGM or PSGM when tackling general instances of problem (1). It would be of interest to fill this gap by identifying or developing practically fast methods that have more general convergence guarantees, especially since different applications may give rise to instances of problem (1) with different structures. In a recent attempt to address these questions, Li et al. [15] showed, among other things, that RSGM will converge at the sublinear rate of O(k -1/4 ) (here, k is the iteration counter) when applied to a general instance of problem (1) and at a linear rate when applied to a so-called sharp instance of problem (1). Informally, an optimization problem is said to possess the sharpness property if the objective function grows linearly with the distance to a set of local minima [16]. Such a property plays a crucial role in establishing fast convergence guarantees for a host of iterative methods; see, e.g., [15]- [17] and also [18]- [20] for related results. Since the ODL problem and the DPCP formulation of the RSR problem are known to possess the sharpness property under certain assumptions on the data [6], [9], the results in [15] imply that RSGM will converge linearly on these problems.
In this paper, we depart from the subgradient-type approaches (such as RSGM (2) and PSGM (3)) and present another method called the manifold proximal point algorithm (ManPPA) to tackle problem (1). At each iterate x k , ManPPA computes a search directon by minimizing the sum of f and a proximal term defined in terms of the Euclidean distance over the tangent space to M at x k . This should be contrasted with other existing PPAs on manifolds (see, e.g., [21], [22]), in which the proximal term is defined in terms of the Riemannian distance. Such a difference is important. Indeed, although the search direction defined in ManPPA does not admit a closed-form formula, it can be computed in a highly efficient manner by exploiting the structure of problem (1); see Section II-B. However, the search direction defined in the existing PPAs on manifolds can be as difficult to compute as a solution to the original problem. Consequently, the applicability of those methods is rather limited.
We now summarize our contributions as follows:
1) We show that ManPPA has a global sublinear convergence rate of O(k -1/2 ) when applied to a general instance of problem (1). Moreover, we show that if the instance has the sharpness property, then the local convergence rate of ManPPA is at least quadratic. Although the sublinear rate result follows from the results in [23], the quadratic rate result is new. Moreover, both rates are superior to those of RSGM established in [15]. Key to
Our interest in problem (4) stems from the observation that it provides alternative formulations of the ODL and RSR problems. Indeed, for the ODL problem, one can recover the entire orthonormal basis all at once by solving problem (4) with q = n. For the RSR problem, if one knows the dimension d of the inlier subspace S, then one can recover it by solving problem (4) with q = n-d. We show that a good feasible solution to problem (4) can be found in a column-by-column manner by suitably modifying ManPPA. Although the proposed method is only a heuristic, our extensive numerical experiments show that it yields solutions of comparable quality to but is significantly faster than existing methods on the ODL and RSR problems.
The rest of the paper is organized as follows. In Section II, we present ManPPA for solving problem (1) and describe a highly efficient method for solving the subproblem that arises in each iteration of ManPPA. We also analyze the convergence behavior of ManPPA. In Section III, we propose StManPPA, a stochastic version of ManPPA that is well suited for large-scale computation, and analyze its convergence behavior. In Section IV we discuss an extension of ManPPA for solving the matrix analog (4) of problem (1). In Section V, we apply ManPPA to solve the ODL problem and the DPCP formulation of the RSR problem and compare its performance with some existing methods. We draw our conclusions in Section VI.
Besides the notation introduced earlier, we use L to denote the Lipschitz constant of f ; i.e., |f (x)-f (y)| ≤ L x-y 2 for all x, y ∈ R n (note that L ≤ √ n max u∈R n Y u 1 / u 1 ). Given a closed set C ⊆ R n , we use Proj C (x) to denote the projection of x onto C and dist(x, C) := inf y∈C yx 2 to denote the distance between x and C. Given a proper lower semicontinuous function h : R n → R ∪ {+∞}, its proximal mapping is given by prox h (x) = argmin w∈R n h(w)
2 . Given two vectors x, y ∈ R n , we use x, y or x y to denote their usual inner product. Other notation is standard.
Since problem (1) is nonconvex, our goal is to compute a stationary point of (1), which is a point x ∈ M that satisfies the first-order optimality condition
(see [14]). In the recent work [23], Chen et al. considered the more general problem of minimizing the sum of a smooth function and a nonsmooth convex function over the Stiefel manifold and developed a manifold proximal gradient method (ManPG) for finding a stationary point of it. When specialized to solve problem (1), the method generates the iterates via
where the search direction d k is given by
and α k > 0, t > 0 are the step sizes. As the reader may readily recognize, without the constraint d x k = 0, the subproblem ( 6) is simply computing the proximal mapping of f at x k and coincides with the update of the classic proximal point algorithm (PPA) [24]. The constraint d x k = 0 in (6), which states that the search direction d should lie on the tangent space to M at x k , is introduced to account for the manifold constraint in problem (1) and ensures that the next iterate x k+1 achieves sufficient decrease in objective value. Motivated by the above discussion, we call the method obtained by specializing ManPG to the setting of problem (1) ManPPA and present its details below:
Algorithm 1 ManPPA for Solving Problem (1)
1: Input: x 0 ∈ M, β ∈ (0, 1), t > 0.
2: for k = 0, 1, 2, . . . do 3:
Solve the subproblem (6) to obtain d k .
Let j k be the smallest nonnegative integer such that
Set x k+1 according to (5) with α k = β j k . 6: end for Naturally, ManPPA inherits the properties of ManPG established in [23]. However, due to the structure of problem (1), many of the developments in [23] have to be refined when designing ManPPA. In particular, the SSN method used by ManPG for finding the search direction in each iteration requires the computation of the proximal mapping of the nonsmooth part of the objective function. However, due to the presence of the matrix Y , the objective function f of problem (1) does not have an easily computable proximal mapping. As such, the SSN method proposed in [23] cannot efficiently solve the subproblem (6). To circumvent this difficulty, we propose to use an inexact ALM, which can efficiently compute an accurate solution to (6); see Section II-B. Now, let us state the following result, which shows that the line search step in line 4 of Algorithm 1 is well defined. It simplifies [23,Lemma 5.2] and yields sharper constants. The proof can be found in Appendix B.
As a result, we have α k = β j k > β ᾱ for any k ≥ 0 in Algorithm 1, which implies that the line search step terminates after at most log β ᾱ + 1 iterations. In particular, if t ≤ 1/L, then we have ᾱ = 1, which implies that we can take j k = 0 in line 4 of Algorithm 1; i.e., no line search is needed.
In this subsection, we study the convergence behavior of ManPPA. Recall from [23,Lemma 5.3] that if d k = 0 in (6), then x k ∈ M is a stationary point of problem (1). This motivates us to call x k ∈ M an -stationary point of problem (1) with ≥ 0 if the solution d k to (6) satisfies d k 2 ≤ . By specializing the convergence results in [23,Theorem 5.5] for ManPG to ManPPA, we obtain the following theorem:
Theorem 1 shows that ManPPA has an iteration complexity of O( -2 ), which is superior to the O( -4 ) bound established for RSGM in [15]. Now, let us analyze the convergence rate of ManPPA in the setting where problem (1) possesses the sharpness property. Such a setting is highly relevant in applications, as both the ODL problem and DPCP formulation of the RSR problem give rise to sharp instances of problem (1) under certain assumptions on the data; see [6,Proposition C.8] and [15,Proposition 4]. To proceed, we first introduce the notion of sharpness. Definition 1 (Sharpness; see, e.g., [16]). We say that X ⊆ M is a set of weak sharp minima for the function f with parameters
From the definition, we see that if X is a set of weak sharp minima of f , then it is the set of minimizers of f over B(δ). Moreover, the function value grows linearly with the distance to X . In the presence of such a regularity property, ManPPA can be shown to converge at a much faster rate. The following result, which has not appeared in the literature before and is thus new, constitutes the first main contribution of this paper.
Theorem 2. Suppose that X ⊆ M is a set of weak sharp minima for the function f with parameters (α, δ). Let {x k } k be the sequence generated by Algorithm 1 with dist(x 0 , X ) < δ := min δ, α L and t ≤ min
Theorem 2 establishes the quadratic convergence rate of ManPPA when applied to a sharp instance of problem (1). Again, this is superior to the linear convergence rate of RSGM established in [15] for this setting. The proof of Theorem 2 can be found in Appendix C. Note that since M is nonconvex, one cannot directly apply standard convergence analysis techniques for PPA (see, e.g., [24]) to obtain Theorem 2. The key to overcoming this diffculty is the new Riemannian subgradient inequality we establish in Proposition 5 (see Appendix C), which provides a path for extending the convergence analysis of PPA to that of ManPPA.
It should be pointed out that the results in Theorems 1 and 2 do not assume any generative model of the data matrix Y . By contrast, the results developed in, e.g., [6] for the ODL problem and [9] for the DPCP formulation of the RSR problem do assume certain generative models of the data. Although the latter results may yield qualitatively sharper convergence guarantees for instances of (1) that arise from the ODL problem or the DPCP formulation of the RSR problem, the former apply to arbitrary instances of (1).
Observe that each iteration of ManPPA requires solving the subproblem (6) to obtain the search direction. Thus, the efficiency of ManPPA depends not only on its convergence rate (which has already been studied in Theorems 1 and 2) but also on how fast the subproblem (6) can be solved. To address the latter, we note that ( 6) is a linearly constrained strongly convex quadratic minimization problem. This motivates us to adopt the SSN-based ALM originally developed in [25] for LASSO-type problems to solve it. As we shall see, such an approach yields a highly efficient method for solving the subproblem (6).
To set the stage for our development, let us drop the index k from (6) for simplicity and set c = Y x. Then, the subproblem (6) can be equivalently written as
At this point, one may be tempted to use ADMM to solve problem (11). However, from a practical point of view, ADMM is often unable to return a high-accuracy solution in an efficient manner. Since ( 11) is a subproblem in ManPPA, a low-accuracy solution will adversely affect the convergence rate of ManPPA. In fact, this has been observed in our numerical experiments. Therefore, we propose to use an inexact ALM, which can solve problem (11) efficiently and accurately. This makes it possible for ManPPA to achieve fast convergence. To describe the algorithm, let us first write down the augmented Lagrangian function corresponding to (11):
Here, y ∈ R and z ∈ R p are Lagrange multipliers (dual variables) associated with the constraints in (11) and σ > 0 is a penalty parameter. Then, a generic iteration of the inexact ALM is given by
where the penalty parameters {σ j } j are chosen such that 0 < σ j σ ∞ ≤ +∞ [25], [26]. Since the subproblem (13a) can only be solved inexactly in general, we adopt the following stopping criteria, which are standard in the literature (see [25], [26]):
Here, Ψ * j is the optimal value of (13a). Equipped with conditions (14a)-(14c), we can show that starting from any initial point (d 0 , u 0 ; y 0 , z 0 ), the inexact ALM (13) will converge at an asymptotic superlinear rate to an optimal solution to problem (11). This result, which constitutes the second main contribution of this paper, is obtained from a new perturbation analysis of the solution set of the subproblem ( 6) and its dual. The details can be found in the companion technical report of this paper [27]. Now, it remains to discuss how to solve the subproblem (13a) in an efficient manner. Again, let us drop the index j in (13a) for simplicity. By simple manipulation, we have
.
Consider the function d → ψ(d) := inf u∈R p Ψ(d, u). Upon letting w = Y d + c + z σ ∈ R p and using the definition of the proximal mapping of u → h(u) := t u 1 , we have
It follows that ( d, ū) = argmin d∈R n ,u∈R p Ψ(d, u) if and only if
Using
Thus, we can find d by solving the nonsmooth equation
Towards that end, we apply an SSN method, which finds the solution by successive linearization of the map ∇ψ. To implement the method, we first need to compute the generalized Jacobian of ∇ψ [29, Definition 2.6.1], denoted by ∂(∇ψ). By the chain rule [29, Corollary of Theorem 2.6.6] and the Moreau decomposition, each element V ∈ ∂(∇ψ) takes the form
where Q ∈ ∂prox h/σ (w). Using the definition of h, it can be shown that the diagonal matrix Q = Diag(q) with
is an element of ∂prox h/σ (w) [25, Section 3.3] and hence can be used to define an element V ∈ ∂(∇ψ) via (16). Note that the matrix V so defined is positive definite. As such, the following generic iteration of the SSN method for solving (15) is well defined:
Here, ρ j > 0 is the step size. Moreover, since I p -Q is a diagonal matrix whose entries are either 0 or 1, the matrix V can be assembled in a very efficient manner; again, see [25].
In the companion technical report [27], we present a detailed implementation of the SSN method ( 17) and show that it will converge at a superlinear rate to the unique solution d to (15).
In this section, we propose, for the first time, a stochastic ManPPA (StManPPA) for solving problem (1), which is well suited for the setting where p (typically representing the number of data points) is much larger than n (typically representing the ambient dimension of the data points). To begin, observe that problem (1) has the finite-sum structure
where y j ∈ R n is the j-th column of Y . When p is extremely large, computing the matrix-vector product Y x can be expensive. To circumvent this difficulty, in each iteration of StManPPA, a column of Y , say y j , is randomly chosen and the search direction d k is given by
The key advantage of StManPPA is that the subproblem ( 18) admits a closed-form solution that is very easy to compute.
Proposition 2. Let µ = t(y j x k ). Then, the solution to (18) is given by
Proof. The first-order optimality conditions of (18) are
Suppose that d is a solution to (20). If
, and (20a) implies that 0 = d/t + y jλx k . This, together with (20b) and the fact that
This establishes the first case in (19). The other two cases in (19), which correspond to y j (x k + d) < 0 and y j (x k + d) = 0, can be derived using a similar argument.
We now present the details of StManPPA in Algorithm 2. It is worth noting that our proposed StManPPA is different from the ones developed in the recent work [30]. Indeed, in each iteration, the former only needs to solve a subproblem that involves a single component of the objective function, while the latter need to compute the proximal mapping of the entire objective function.
Algorithm 2 StManPPA for Solving Problem (1) 1: Input: x 0 ∈ M, t 0 , t 1 , . . . , t T > 0. 2: for k = 0, 1, . . . , T do 3:
Select j k ∈ {1, . . . , p} uniformly at random and solve the subproblem (18) with j = j k , t = t k to obtain d k .
Set x k+1 = Proj M (x k + d k ). 5: end for 6: Output: x = x k with probability t k / T k=0 t k .
In this section, we present our convergence results for StManPPA. Let us begin with some preparations. Define f j : R n → R to be the function f j (x) = y j x and let L j > 0 denote the Lipschitz constant of f j , where j = 1, . . . , p. Set L := max j∈{1,...,p} L j . Furthermore, define the Moreau envelope and proximal mapping on M by
respectively. The proximal mapping mprox is well defined since the constraint set M is compact. As it turns out, the proximal mapping mprox can be used to define an alternative notion of stationarity for problem (1). Indeed, for any λ > 0 and x ∈ M, if we denote x = mprox λf (x), then the optimality condition of (21) yields
Since I n -x x is a projection operator and hence nonexpansive, we obtain
In particular, if 1 λ xx 2 ≤ , then (i) x is -stationary in the sense that dist(0, ∂ R f ( x)) ≤ and (ii) x is close to the -stationary point x. This motivates us to use
as a stationarity measure for problem (1). We call x ∈ M an -nearly stationary point of problem (1) if Θ λ (x) ≤ . It is worth noting that such a notion has also been used in [15] to study the stochastic RSGM.
We are now ready to establish the following convergence rate result for StManPPA, which constitutes the third main contribution of this paper. Theorem 3. For any λ ∈ (0, 1/(p L)), the point x output by Algorithm 2 satisfies
where the expectation is taken over all random choices made by the algorithm. In particular, if the step sizes {t k } k satisfy
The proof of Theorem 3 can be found in Appendix D. Again, it makes crucial use of our newly established Riemannian subgradient inequality (see Appendix C, Proposition 5), which allows StManPPA to be analyzed in a similar way as various Euclidean stochastic methods. We remark that the iteration complexity bound O( -4 ) of StManPPA established in Theorem 3 is comparable to that of RSGM established in [15].
In this section, we consider the matrix analog (4) of problem (1), which also arises in many applications such as certain "one-shot" formulations of the ODL and RSR problems (see Section I-B). Currently, there are two existing approaches for solving problem (4), namely a sequential linear programming (SLP) approach and an iteratively reweighted least squares (IRLS) approach [8]. In the SLP approach, the columns of X are extracted one at a time. Suppose that we have already obtained the first columns of X ( = 0, 1, . . . , q -1) and arrange them in the matrix X ∈ R n× (with X 0 = 0). Then, the ( + 1)-st column of X is obtained by solving
This is achieved by the alternating linearization and projection (ALP) method, which generates the iterates via
Note that the z-subproblem is a linear program, which can be efficiently solved by off-the-shelf solvers.
In the IRLS approach, the following variant of problem (4), which favors row-wise sparsity of Y X and has been studied by Lerman et al. in [31], is considered:
Here, Y X 1,2 denotes the sum of the Euclidean norms of the rows of Y X. The IRLS method for solving (22) generates the iterates via
where w j,k = 1/ max{δ, X k-1 y j 2 } and δ > 0 is a perturbation parameter to prevent the denominator from being 0. The solution X k can be obtained via an SVD and is thus easy to implement. However, there has been no convergence guarantee for the IRLS method so far. Recently, Wang et al. [32] proposed a proximal alternating maximization method for solving a maximization version of (4), which arises in the so-called 1 -PCA problem (see [7]). However, the method cannot be easily adapted to solve problem (4).
The similarity between problems (1) and ( 4) suggests that the latter can also be solved by ManPPA, which is indeed the case. However, the SSN method for solving the resulting nonsmooth equation (i.e., the matrix analog of ( 15)) can be slow, as the dimension of the linear system (17) is high. Here, we propose an alternative method called sequential ManPPA, which solves (4) in a column-by-column manner and constitutes the fourth main contribution of this paper. The method is based on the observation that the objective function in (4) is separable in the columns of X = [x 1 , x 2 , . . . , x n ]. To find x 1 , we simply solve min
using ManPPA as it is an instance of problem (1). Suppose that we have found the first columns of X ( = 0, 1, . . . , q -1) and arrange them in the matrix Q ∈ R n× (with Q 0 = 0). Then, we find the ( + 1)-st column x +1 by solving
(23) As it turns out, problem ( 23) is equivalent to the deflation process discussed in [11] for sequentially recovering the columns of an orthogonal dictionary. Specifically, let V be an orthonormal basis of the orthogonal complement of Q . We can then find x +1 by solving
Note that ( 24) has the same form as (1) and thus can be solved by RSGM or PSGM. By contrast, problem (23) has an extra linear constraint and cannot be solved by RSGM or PSGM directly. Nevertheless, our ManPPA can solve both ( 23) and (24). Let us now briefly discuss how to use ManPPA to solve the former. To simplify notation, let us drop the index and denote x = x +1 , Q = Q . In the k-th iteration of ManPPA, we update the iterate by (5), where the search direction d k is computed by
This subproblem can be solved using the SSN-based inexact ALM framework in Section II. We omit the details for succinctness. The sequential ManPPA is guaranteed to find a feasible solution to problem (4). Moreover, as we shall see in Section V, it is computationally much more efficient than the ALP and IRLS methods on the ODL problem and the DPCP formulation of the RSR problem. However, it remains open whether the solution found by sequential ManPPA is a stationary point of problem (4). We leave this question for future research.
In this section, we compare our proposed ManPPA, StManPPA, and sequential ManPPA with the existing methods PSGM-MBLS, ALP, and IRLS. We do not include RSGM with diminishing step sizes in our comparison, as the numerical results in [9] show that they are slower than PSGM-MBLS and cannot achieve high accuracy. In all the tests, we used the step size t = 0.1 for ManPPA and set the maximum number of iterations of ManPPA, inexact ALM, and SSN to 100, 30, and 20, respectively. We stopped ManPPA if the relative change in the objective values satisfies |f (x k ) -f (x k-1 )|/f (x k-1 ) ≤ 10 -9 . In the k-th iteration of ManPPA, we stopped the inexact ALM if the primal and dual feasibility of problem (11)
For SSN, we used the same termination criteria as the ones given in [25] and solved the linear equation (17a) by Cholesky decomposition. We refer the reader to the companion technical report [27] for the setting of the parameters in the inexact ALM and SSN.
For StManPPA, we used the piecewise geometrically diminishing step sizes t k = β k/p t 0 for k = 0, 1, . . . , pT with t 0 = 0.6 and T = 500. Such step sizes are motivated by those used in PSGM [9]. We use StManPPA-β to specify the parameter β used in the algorithm. We stopped the algorithm if the relative change in the objective values satisfies |f (x k ) -f (x k-1 )|/f (x k-1 ) ≤ 10 -12 . For PSGM-MBLS, ALP, and IRLS, we used their default settings of the parameters. We stopped ALP if the change in the objective values satisfies |f (x k ) -f (x k-1 )| ≤ 10 -6 , while we stopped IRLS if the change in the objective values satisfies |f (x k ) -f (x k-1 )| ≤ 10 -11 . With these stopping criteria, the solutions returned by these algorithms usually achieve the same accuracy.
In this section, we consider the DPCP formulations of the RSR problem. For the recovery of a vector that is orthogonal to the inlier subspace (the vector case), we compared the performance of ManPPA, StManPPA, PSGM-MBLS, ALP, and IRLS on problem (1). For the recovery of the entire inlier subspace of known dimension d (the matrix case), we compared the performance of sequential ManPPA, PSGM-MBLS, ALP, and IRLS on problem (4) with q = n -d.
1) Synthetic Data: We first test the algorithms on synthetic data. The data matrix takes the form Y = [X, O] ∈ R n×(p1+p2) . We generated the inlier data X as X = QC, where Q ∈ R n×d is an orthonormal matrix and C ∈ R d×p1 is a coefficient matrix. The matrix Q was generated by orthonormalizing an n × d standard Gaussian random matrix via QR decomposition, while the matrix C was generated as a d × p 1 standard Gaussian random matrix. The outlier data O ∈ R n×p2 was generated as a standard Gaussian random matrix. Finally, the columns of the matrix Y were normalized. The d-dimensional subspace spanned by X is denoted by S and its orthogonal complement is denoted by S ⊥ .
Vector case. We set the initial point x 0 of all algorithms to be the eigenvector of Y Y corresponding to the eigenvalue with minimum magnitude. We compared the performance of the algorithms on problem (1) with different dimension n, number of inliers p 1 , and number of outliers p 2 in Figure 1. The first row of Figure 1 reports the principal angle 1 1 The principal angle is the distance between x k and S ⊥ . Any optimal solution x * to problem (1) is orthogonal to the inlier subspace S [9, Theorem 1]. Using the Lipschitz continuity of f , we know that θ also measures the function value gap f (x) -f (x * ).
θ between x k and S ⊥ versus the iteration number. The second row reports θ versus CPU time. Note that x k = sin(θ)n + cos(θ)s, where n = Proj S (x k )/ Proj S (x k ) 2 and s = Proj S ⊥ (x k )/ Proj S ⊥ (x k ) 2 . From Figure 1, we see that PSGM-MBLS is the fastest, while ManPPA is slightly slower. However, they are both much faster than other compared methods. Moreover, the principal angle θ of ManPPA decreases much faster than PSGM-MBLS in terms of iteration number. This can be attributed to the quadratic convergence rate of ManPPA (Theorem 2). In Figure 2 we report the quadratic fitting curves of the different algorithms. As shown in [9], the DPCP formulation can tolerate O((#inliers) 2 ) outliers; i.e., p 2 = O(p 2 1 ). For different p 2 ∈ {40, 80, 120, . . . , 600}, we find the smallest p 1 ∈ {60, 70, 80, . . . , 260} such that θ < 10 -1 . Here, the principal angle θ is the mean value of 10 trials; i.e., we find pairs (p 1 , p 2 ) such that for a fixed p 1 , p 2 is the largest number of outliers that can be tolerated. We then use a quadratic function to fit these pairs (p 1 , p 2 ). A higher curve indicates that more outliers can be tolerated and hence the algorithm is more robust. From Figure 2, we see that the curve corresponding to PSGM-MBLS is the lowest one and thus the least robust, while ManPPA and StManPPA are more robust. In Figure 3 we report the CPU time versus p 1 and p 2 . For each algorithm, the shadow area represents the standard deviation (std) of 10 random trials, while the line within the shadow is the mean of those trials. From the left two subfigures of Figure 3, we see that the stds of IRLS and ALP are quite significant. In particular, they are usually more than ten times larger than the stds of other compared algorithms. To better illustrate the stds of the other four algorithms, we plot their CPU times in the right two subfigures of Figure 3. From these two figures, we see that ManPPA has a larger std than those of the other three algorithms. Overall, we see that PSGM-MBLS is the fastest and ManPPA is second, and they are both much faster than the other compared algorithms. Figures 2 and3 suggest that ManPPA is slightly slower than PSGM-MBLS but is more robust. Moreover, the choice of the parameter β for StManPPA is crucial and challenging, as StManPPA-0.9 is faster but less robust than StManPPA-0.8. We leave the determination of the best parameter β for StManPPA as a future work. Matrix case. We solved problem (4) using sequential ManPPA and compared its performance with PSGM-MBLS (applied to (24)), ALP, and IRLS. We report the results for q = 2 and q = 4 in Figures 4 and5, respectively. Since ALP has a very high std, for better illustration of the CPU time comparison, we exclude it from the right two subfigures of Figures 4 and5. The results suggest that sequential ManPPA is not as robust as IRLS but is the second most efficient one among the four compared algorithms. Moreover, we find that although PSGM-MBLS is very efficient, its fitting curve is not very good. This is due to the fact that PSGM-MBLS is very sensitive to the choice of step size.
2) Real 3D Point Cloud Road Data: Next, we compared ManPPA with PSGM-MBLS on the road detection challenge of the KITTI dataset [33]. This dataset contains image data together with the corresponding 3D points collected by a rotating 3D laser scanner. Similar to [9], we only used the 360°3 D point clouds to determine which points lie on the road plane (inliers) and which do not (outliers). By using homogeneous coordinates, this can be cast as a robust hyperplane learning problem (1) in R 4 . As reported in [9], PSGM-MBLS is the fastest algorithm when compared with other state-of-the-art methods. Thus, we only compared the performance of ManPPA and PSGM-MBLS on problem (1). Table I reports the area under the Receiver Operator Curve (ROC) and the CPU time. We see that all ROC values of ManPPA are better than those of PSGM-MBLS, with some sacrifice on the CPU time.
We generated instances of the ODL problem by first randomly generating an orthogonal matrix X ∈ R n×n and a Bernoulli-Gaussian matrix  ∈ R n×p with parameter γ (see, e.g., [3]), then setting Y = X Â.
Vector case. We first compared the performance of ManPPA and StManPPA with ALP, IRLS, and PSGM-MBLS on problem (1) with n = 30, p = 10n 1.5 . Figures 6 and 7 report the iteration numbers and CPU times of the compared algorithms. The quantity θ is the angle between x k returned by the algorithm and its nearest column in X. From Figures 6 and7, we see that PSGM-MBLS is the fastest algorithm in terms of CPU time, while ManPPA is slightly slower. However, they are both much faster than the other compared algorithms. Moreover, we see that ManPPA is much faster than PSGM-MBLS in terms of iteration number. This again can be attributed to the quadratic convergence rate of ManPPA (Theorem 2).
In Figures 8 and9 we report the linear fitting curves for log(n) and log(p) and CPU times of ManPPA, StManPPA-0.8, StManPPA-0.9, IRLS, ALP, and PSGM-MBLS. Note that for the ODL problem, it has been found empirically in [6] that the sample size p and dimension n should satisfy Here the principal angle θ is the mean value of 10 trials. We then use a linear function to fit the points {(log(n), log(p))} n,p . From Figures 8 and9, we find that PSGM-MBLS is the fastest but its fitting curve is high, which suggests that it is not robust. This is because PSGM-MBLS is very sensitive to the choice of step size. ManPPA appears to be the second fastest but is very robust based on the fitting curve. Matrix case. To find the entire orthogonal basis, we use sequential ManPPA to solve problem (4) with q = n. We compared sequential ManPPA with PSGM-MBLS (applied to (24)) and SLP based on ALP, and report the results in Figures 10 and11. We see that sequential ManPPA is not as robust as ALP, but it is much faster than ALP. Note that there is nothing for IRLS to do, as it tackles the objective function Y X 1,2 , which is a constant when X X = I n . We also find that the fitting curve of PSGM-MBLS is very high, which indicates that it fails to recover the dictionary in many cases. Again, this is due to the high sensitivity of PSGM-MBLS to the choice of step size.
In this paper, we presented ManPPA and its stochastic variant StManPPA for solving problem (1). By exploiting the manifold structure of the constraint set M, these methods not only are practically efficient but also possess convergence guarantees that are provably superior to those of existing subgradient-type methods. Using ManPPA as a building block, we also proposed a new sequential approach to solving the matrix analog (4) of problem (1). We conducted extensive numerical experiments to compare the performance of our proposed algorithms with existing ones on the ODL problem and DPCP formulation of the RSR problem. The results demonstrated the efficiency and efficacy of our proposed methods.
In this section, we collect some useful properties of the projector Proj M . Proposition 3. For any x ∈ M and d ∈ R n satisfying d x = 0 (i.e., d is a tangent vector at x), we have Moreover, if d 2 ≤ D for some D ∈ (0, +∞), then
Proof. It is straightforward to verify that
We then have (25) by using the fact that
we get (26). Proposition 4. For any x, z ∈ M and d ∈ R n satisfying d x = 0, we have Proof. We compute
where (27) follows from the fact that x + d
Since f is Lipschitz with constant L and d k x k = 0, we have by Proposition 3. Hence, for any α ∈ (0, ᾱ], we have
where (28a) follows from the convexity of f , (28b) holds because the strong convexity of the objective function in subproblem ( 6), together with the optimality of d = d k and feasibility of d = 0 for (6), implies that f and(28c) is due to α ≤ 1/(tL). If t ≤ 1/L, then ᾱ = 1. This completes the proof.
We begin with two preparatory results. The first states that the restriction of the objective function f in (1) on the nonconvex constraint set M satisfies a Riemannian subgradient inequality, which means that with respect to its Riemannian subgradient, the function f behaves almost like a convex function on M. Proposition 5. Let x ∈ M and d ∈ R n be such that d x = 0. Define x + = x + d. Then, for any z ∈ M and s ∈ ∂f (x + ), we have Proof. Since f is convex on R n , we have
where the second equality is due to x x = 1 and d x = 0, and the inequality follows from the fact that |x z -1| = The second establishes a key recursion for the iterates generated by ManPPA. Proposition 6. Let {x k } k be the sequence generated by Algorithm 1 with t ≤ 1/L. For any x ∈ M, we have
Proof. Since t ≤ 1/L, we have x k+1 = Proj M (x k + d k ) by from Proposition 1. From the optimality condition of the subproblem (6), there exists an
Denoting x k + = x k + d k , we have We are now ready to prove Theorem 2. We first prove (9) by induction. Let x * ∈ X be such that dist(x k , X ) = x k -x * 2 . By invoking Proposition 6 with x = x * , we have
where the last inequality follows from (8). Consider the function [0, δ] s → φ(s) = (1 + tL)s 2 -2tαs + t 2 L 2 . Observe that φ attains its maximum at s = δ if δ ≥ 2tα 1+tL . Given that t ≤ min In particular, we have dist(x k+1 , X ) ≤ δ whenever dist(x k , X ) ≤ δ. This establishes (9). Next, we prove (10). Again, let x * ∈ X be such that dist(x k , X ) = x kx * 2 . Since α ≤ L by (8), we have δ ≤ 1. This implies that x k x * = 2-x k -x * 2 2 2 ≥ 1 2 . Hence, the vector dk = x *
x k x * -x k is well defined and satisfies dk x k = 0 (i.e., dk is a tangent vector at x k ), Proj M (x k + dk ) = x * , and dk 2 ≤ √ 3. By the strong convexity of the objective function in subproblem (6) and noting the optimality of d k and feasibility of dk for (6), we have
Furthermore, by the Lipschitz continuity of f and Proposition 3, we get
Combining ( 31), (32a), and (32b), we have
Since t ≤ δ 2α-Lδ ≤ 1/L, we have 1 2t -L 2 ≥ 0. Moreover, since dist(x k+1 , X ) ≤ δ ≤ δ, we have f (x k+1 ) -f (x * ) ≥ α • dist(x k+1 , X ) by (8). It then follows from (33) and Proposition 3 that
This completes the proof. (34)
From the optimality condition of (18), we get
Hence, we compute
where (35a) follows from Proposition 5 and (35b) is due to the Lipschitz continuity of f j k . Upon taking the expectation on both sides of (35) with respect to j k conditioned on x k , we obtain
where (36a) follows from the fact that Upon summing the above inequality over k = 0, 1, . . . T and noting that λ < 1/(p L) and e λ (z) ≥ 0 for any z ∈ R n , we obtain
Upon dividing both sides of the above inequality by T k=0 t k and noting that the left-hand side becomes E Θ λ ( x) 2 , the proof is complete.