Member, IEEE Arnaud Breloy

Sandeep Kumar

Ying Sun

Fellow, IEEE Daniel P Palomar

Majorization-Minimization on the Stiefel Manifold with Application to Robust Sparse PCA

Keywords: Majorization-Minimization, Stiefel manifold, Grassmann manifold, Subspace learning, PCA, sparse PCA

This paper proposes a framework for optimizing cost functions of orthonormal basis learning problems, such as principal component analysis (PCA), subspace recovery, orthogonal dictionary learning, etc. The optimization algorithm is derived using the majorization-minimization framework in conjunction with orthogonal projection reformulations to deal with the orthonormality constraint in a systematic manner. In this scope, we derive surrogate functions for various standard objectives that can then be used as building blocks, with examples for robust learning costs and sparsity enforcing penalties. To illustrate this point, we propose a new set of algorithms for sparse PCA driven by this methodology, whose objective function is composed of an M -estimation type subspace fitting term plus a regularizer that promotes sparsity. Simulations and experiments on real data illustrate the interest of the proposed approach, both in terms of performance and computational complexity.

I. INTRODUCTION

O PTIMIZATION problems on the Stiefel manifold, i.e., involving the orthonormality constraint on the matrix variable, are ubiquitous in signal processing and machine learning. To name a few, such problems arise in principal component analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis and Factor Analysis[END_REF][START_REF] Tipping | Probabilistic principal component analysis[END_REF], sparse PCA [START_REF] Chen | Reduced rank stochastic regression with a sparse singular value decomposition[END_REF][START_REF] Chen | Sparse reduced-rank regression for simultaneous dimension reduction and variable selection[END_REF][START_REF] Bunea | Joint variable and rank selection for parsimonious estimation of high-dimensional matrices[END_REF][START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF][START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF][START_REF] Uematsu | SOFAR: Largescale association network learning[END_REF], structured covariance matrix estimation [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF][START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF], robust subspace recovery [START_REF] Lerman | An overview of robust subspace recovery[END_REF][START_REF] La | A framework for robust subspace learning[END_REF][START_REF] Maronna | Principal components and orthogonal regression based on robust scales[END_REF][START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF][START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF][START_REF] Maunu | A well-tempered landscape for non-convex robust subspace recovery[END_REF], and subspace clustering [START_REF] Vidal | Subspace clustering[END_REF]. To tackle the challenges from big data and high-dimensional settings the solutions are often desired with certain additional properties, e.g., sparsity. These additional properties are often enforced by including regularization penalties to the objective function. Except for some special cases where the solution appears as eigenvectors of a given matrix, these problems are usually nontrivial to deal with, owing to the nonconvex orthonormality constraints. The resolution of these problems thus generally calls for the use of iterative constrained optimization methods to search for their local minima. A very general framework to account for the orthonormality constraint is brought by Riemannian optimization algorithms on the Stiefel manifold [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Boumal | Manopt, a matlab toolbox for optimization on manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF][START_REF] Abrudan | Steepest descent algorithms for optimization under unitary matrix constraint[END_REF][START_REF]Conjugate gradient algorithm for optimization under unitary matrix constraint[END_REF][START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF]. However, the algorithms derived within this perspective can turn to be computationally costly. Additionally, ensuring sparsity and orthonormality simultaneously is a challenging goal, which generally calls for a trade-off between these two properties in Arnaud Breloy is with LEME (EA4416), University Paris Nanterrre, France (e-mail: abreloy@parisnanterre.fr). Sandeep Kumar is with Department of Electrical Enginerring, Indian Institute of Technology Delhi, India (email: sandeep0kr@ee.iitd.ac.in). Ying Sun is with the Pennsylvania State University, USA (email: ysun@psu.edu). Daniel P. Palomar is with the Hong Kong University of Science and Technology (HKUST), Hong Kong (e-mail: palomar@ust.hk). This work was supported by the Hong Kong GRF 16207019 research grant.

the solutions obtained by existing state-of-the-art approaches (cf. [START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF] and reference therein). A recent line of work based on orthogonal Procrustes reformulations, however, has enabled researchers to ensure both of these properties. These works addressed specific applications, e.g., [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF] for the sparse PCA and covariance estimation and [START_REF] Uematsu | SOFAR: Largescale association network learning[END_REF] for low-rank matrix estimation problems.

Motivated by the applicability of these reformulations for tackling the aforementioned issues, this work aims to develop a general low-complexity algorithmic framework for largescale sparsity-regularized optimization problems on the Stiefel Manifold. This framework will be based on the majorizationminimization (MM) algorithm [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF], which proceeds with two steps: i) (majorization) finding a function that locally upperbounds the objective function up to a constant, referred to as surrogate function; ii) (minimization) minimizing this surrogate function. This procedure generates a sequence that monotonically decreases the objective value, and its main interest is that, with properly chosen surrogate functions, it can yield a sequence of subproblems that are easy to deal with. The main idea will thus be to deal with the orthonormality constraint through a systematic formulation of the surrogates minimization subproblems as a Euclidean projection onto the Stiefel manifold [START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF][START_REF] Schönemann | A generalized solution of the orthogonal Procrustes problem[END_REF][START_REF] Eldén | A Procrustes problem on the stiefel manifold[END_REF], as done for specific cases in, e.g., [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF][START_REF] Koschat | A weighted Procrustes criterion[END_REF][START_REF] Kiers | Setting up alternating least squares and iterative majorization algorithms for solving various matrix optimization problems[END_REF][START_REF] Breloy | Robust rank constrained kronecker covariance matrix estimation[END_REF]. Thereby, we propose a unified approach with convergence guarantee and a set of practical guidelines for applying this methodology that generalizes to a large class of cost functions. Notably, an emphasis is put on robust data fitting cost functions and sparsity enforcing penalties. The advantages of the proposed methodology are multiple: i) it can be applied to a large class of standard cost functions; ii) it guarantees that the iterates are orthonormal and has standard MM convergence guarantees, i.e., monotonic decrement of the objective value, and convergence to the set of critical points of the problem [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF]; iii) it yields scalable algorithms, as the computational bottleneck of each iteration lies only in the computation of a thin-SVD.

As a main example to illustrate this framework, we propose a new class of algorithms for sparse PCA, referred to as RSPCA (for robust sparse PCA). The familly of objective functions proposed for this task combines a M -estimation-type subspace fitting function [START_REF] La | A framework for robust subspace learning[END_REF][START_REF] Maronna | Principal components and orthogonal regression based on robust scales[END_REF][START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF][START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF], plus a sparsity promoting penalty. The considered penalties leverage the proxies of the 0 -norm proposed in [START_REF] Song | Sparse generalized eigenvalue problem via smooth optimization[END_REF], and allow for the optimization to be conducted with the proposed methodology. Interestingly, the resulting sparse PCA algorithms do not involve a relaxation of the orthonormality constraint, which is usually required in the existing state-of-the-art methods [START_REF] Chen | Reduced rank stochastic regression with a sparse singular value decomposition[END_REF][START_REF] Chen | Sparse reduced-rank regression for simultaneous dimension reduction and variable selection[END_REF][START_REF] Bunea | Joint variable and rank selection for parsimonious estimation of high-dimensional matrices[END_REF][START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF]. Finally, simulation results illustrate the usefulness of the proposed approaches, both in terms of performance and computational complexity.

The presentation is organized as follows. Section II presents relevant background and the motivations for this study. Section III presents a generic MM framework for optimization over the Stiefel manifold, where a systematic linear majorization is applied in order to derive simple iterations. Such majorization can be obtained for various standard cost functions, for which a catalog is presented in Section IV. Section V presents several examples of algorithms derived following this framework, and exhibits some links between them and existing state-of-the-art algorithms. Section VI illustrates the application of this framework to sparse-PCA for a flexible class of sparsity promoting penalties. This adds to the list of functions covered by Section IV while discussing methodological insights for applying the framework to non-standard costs functions. Finally, Section VII presents some simulations and experiments on real data to illustrate the interest of the proposed approach and algorithms, both in terms of performance and computational complexity.

The following notation is adopted: italic, lower case boldface, and upper case boldface indicate respectively, scalar, vector, and matrix quantities. The upperscript H denotes the transpose conjugate operator. I is the identity matrix of the appropriate dimension. Tr{} is the Trace operator.

{w i } n i=1 denotes the set of elements w i , ∀i ∈ [[1, N ]]. diag({a i } n i=1 ) is the n × n diagonal matrix with diagonal entries a n . For tall matrices M ∈ C p×k (p > k) the thin-SVD (TSVD) is denoted M TSVD = U k D k V H ,

II. BACKGROUND AND POSITIONING

A. Majorization-Minimization (MM) algorithms

The MM framework is briefly reviewed below. For more complete information, we refer the reader to [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF]. Consider the following optimization problem:

minimize x∈X f (x), (1)
where f : X → R is a continuous function and X is a closed set. Given an initial point x 0 ∈ X , the MM procedure minimizes f over X by updating x iteratively as

x t+1 ∈ argmin x∈X g(x|x t ), (2)
where g(•|x t ) : X → R is a surrogate function of f satisfying the following property:

x t ∈ argmin x∈X g(x|x t ) -f (x). (3)
In other words, g(•|x t ) upperbounds f globally over set X up to a constant:

g(x|x t ) -f (x) ≥ c t {g(x t |x t ) -f (x t )}, ∀x ∈ X . ( 4
)
The sequence {f (x t )} t∈N generated by ( 2) is non-increasing since

f (x t+1 ) (4) ≤ g(x t+1 |x t ) -c t (2) ≤ g(x t |x t ) -c t = f (x t ). ( 5
)
The MM procedure suggests thus the possibility of minimizing f by iteratively seeking for a sequence of surrogate functions {g(•|U t )} t∈N that are easy to minimize over the feasible set.

B. Stiefel manifold and orthogonal projection problem

First, define the complex Stiefel manifold as:

St(p, k) = U ∈ C p×k | U H U = I for k ≤ p. (6)
A point U ∈ St(p, k) is a semi-unitary1 matrix, referred to as k-orthogonal frame. This point U is also an orthonormal basis that spans a k-dimensional subspace. However, notice that this representation of a subspace is not unique as UQ spans the same subspace for any

Q ∈ St(k, k).
Let X ∈ C p×k be a rank p matrix, the Euclidean projection of X onto the Stiefel manifold is the solution of the following problem minimize

U∈St(p,k) ||X -U|| 2 F . ( 7
)
When X is full rank the problem (7) has a unique global minimizer, given by [START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF]Prop. 7], as

U = P St {X} , (8)
where the operator P St is defined in Algorithm 1.

C. Optimization on St(p, k): some existing solutions Consider a generic optimization problem where the variable is constrained to the Stiefel manifold:

minimize U∈St(p,k) f (U) , (9)
where f : C p×k → R is a differentiable objective function suited to an application of interest. Notice that optimization problems over the Stiefel manifold St(p, k) are nonconvex due to the orthonormality constraint. Hence, they are usually hard to deal with, even for apparently simple objective functions f . A natural way to handle the orthonormality constraint is to turn to the framework of Riemannian optimization. A review of this general topic can be found in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], where St(p, k) is used as example throughout the book. More focused reviews of manifold-oriented counterparts of classical optimization algorithms on St(p, k) (steepest descent, conjugate gradient, Newton's method, etc.) can be found in [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF][START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF] ( [START_REF] Abrudan | Steepest descent algorithms for optimization under unitary matrix constraint[END_REF][START_REF]Conjugate gradient algorithm for optimization under unitary matrix constraint[END_REF] for St(p, p)). Notice that when the function is rotation invariant, i.e., f (U) = f (UQ), ∀Q ∈ St(k, k), it can be minimized on the Grassmann manifold Gr(p, k) (the set of p-dimensional subspaces of C p ), with corresponding Riemannian optimization algorithm [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF][START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF], which has the advantage of reducing the dimension of the problem by exploiting its invariance. In this case, one can still apply optimization algorithms developed for St(p, k), but the variable U ∈ St(p, k) simply acts as a representation of its equivalence class (a point in Gr(p, k)).

A second-more specific-option, can be to split f as f = f u + f v and rewrite [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF]

as minimize V,U∈St(p,k) f u (U) + f v (V) subject to U = V, (10)
then relax the equality constraint of this problem using the augmented Lagrangian method. The aim of this reformulation, e.g., used in [START_REF] Uematsu | SOFAR: Largescale association network learning[END_REF], is mainly to derive a block-coordinate algorithm where the updates with respect to U and V are separately easy to obtain.

D. Aim of this paper

In this paper, we explore the use of the MM framework to solve problems on St(p, k) formulated as in [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF]. Generally speaking, a majorization can be applied to cast subproblems that are more practical for the aforementioned algorithms to be applied. For example, the subproblems could have easy Newton steps, or a simpler update of the block U in problem [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF]. Here, our main focus will be to derive surrogates functions that formulate minimization subproblems with closed-form solutions. Interestingly this systematic approach yields lowcomplexity algorithms for solving problem as in [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF], which can also be used as subroutine for more complex problems where the formulation (10) is involved.

III. MM ALGORITHMS ON THE STIEFEL MANIFOLD WITH SYSTEMATIC PROJECTION REFORMULATIONS

A. Generic algorithm formulation

To tackle the optimization of the cost function in (9) under the orthonormality constraint, we consider applying the MM [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF] framework and minimizing f by solving a sequence of Euclidean projections on St(p, k). In short, we will construct surrogate functions that are linear when restricted to the feasible set St(p, k). The corresponding subproblem can then be recast as a Euclidean projection onto St(p, k) (cf. section II-B), which has a unique solution that can be obtained via thin-SVD. In the following, we assume that the objective in ( 9) is majorized at point U t by a surrogate g(U|U t ) that satisfies the following properties (examples are given later).

Assumption 1. The surrogate function g : C p×k × C p×k → R satisfies the following conditions:

i) Tightness: g(U t |U t ) = f (U t ), ii) Continuity: g(• | •) is continuous on C p×k × C p×k , iii) Upperbound: g(U|U t ) ≥ f (U), ∀U ∈ St(p, k),
iv) Linearity: restricting to St(p, k), g can be expressed as

g U|U t = -Tr R H U t U -Tr U H R U t + const., = -2Re Tr{U H R U t } + const., (11)
where R : C p×k → C p×k is a matrix function of U t . Following the MM procedure described in Section II-A, an update of the parameter U is given by

U t+1 ∈ argmin U∈St(p,k) g U|U t . ( 12
)
Since g is linear (cf. [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]) and U ∈ St(p, k), it is not hard to see that obtaining this update is equivalent to solving minimize

U∈St(p,k) ||R U t -U|| 2 F , (13)
1: initialize t = 0, U (0) ∈ C N ×R 2: repeat 3:
Compute R (U t ) from surrogate [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] 4:

Update U t+1 = P St {R (U t )} with Algorithm 1 5: 2 , the problem (13) admits a unique solution, leading to the MM update:

t = t + 1 6: until convergence criterion is met which is the projection problem discussed in Section II-B. When R (U t ) is full rank
U t+1 = P St R U t , (14)
where the operator P St is defined in Algorithm 1. Eventually, solving the sequence of projection problems results in an MM procedure to optimize f under the orthonormality constraint, which is summarized in Algorithm 2.

Remark 1. Importantly, the MM approach is also applicable to objective function consisting sum of functions of the form

f (U) = I i=1 f i (U). (15)
Then, if each f i can be majorized by a linear surrogate g i of the form

g i U|U t = -Tr R H i U t U -Tr U H R i U t + const., (16)
following the same steps as ( 12)-( 14), the MM updates can simply be obtained as

U t+1 = P St I i=1 R i U t . (17)
Obviously, this methodology cannot be applied to any arbitrary cost function. Still, Section IV presents a catalog of surrogate functions satisfying [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] for a large set of standard cost functions that can be used as building blocks (with several examples detailed in Section V). The method also suggests that simple algorithms can be obtained by designing meaningful proxies of the desired function that can be majorized by a linear surrogate on St(p, k). As a concrete example, we propose sparsity promoting penalties that satisfy this property in Section VI, leading to an original algorithm for sparse PCA.

Remark 2. Quadratic surrogate functions, if existing, can also be used within the considered framework. These surrogates should satisfy Assumption 1 with the following alternate condition: iv ) Restricting to St(p, k), g can be expressed as

g U|U t = -Tr U H M U t UD U t + const., (18)
where M :

C p×k → H + p and D : C p×k → D k are matrix functions of U t .
Indeed, such surrogates also yield tractable minimization subproblems under the orthonormality constraint: the corresponding updates are then obtained from the k eigenvectors associated to the largest eigenvalues of M (U t ) if D (U t ) 0 (smallest conversely). In this paper, we focused on linear surrogates since quadratic functions always admit linear surrogates on St(p, k) (cf. propositions 2 and 3). Thus, functions with quadratic majorizers can also be used in conjunction with function that only have linear ones within the proposed framework, while the converse is not true. Nevertheless, quadratic surrogates can in some cases provide tighter approximations that yield a faster convergence of the MM algorithm (at the cost of handling larger matrices at each iterations). A main example discussing the two options is detailed in section V-C and VII-A.

B. Convergence analysis

The convergence analysis of Algorithm 2 will be obtained by following the one of the successive upper-bound minimization (SUM) algorithm in [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF]. Note that the result of [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF] does not hold directly for Algorithm 2, as the SUM framework does not cover non-convex constraints. Nevertheless, this result can be adapted to St(p, k) as in [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF][START_REF] Fu | Scalable and flexible multiview MAX-VAR canonical correlation analysis[END_REF][START_REF] Kumar | A unified framework for structured graph learning via spectral constraints[END_REF], leading to the following proposition.

Proposition 1. Let {U t } ∞ t=0 be a sequence generated by Algorithm 2. Then the following hold:

1) The sequence {f (U t )} t∈N converges.

2) Every limit point U * of the sequence is a critical (also referred to as Karush-Kuhn-Tucker, or KKT) point of the problem (9). 3) The whole sequence converges to K, the set of KKT points of the problem (9).

Proof. First, remark that Algorithm 2 follows the MM procedure, by [START_REF] Bunea | Joint variable and rank selection for parsimonious estimation of high-dimensional matrices[END_REF] we have that the sequence {f (U t )} t∈N is nonincreasing, and thus converges. The compactness of the set St(p, k) then implies that the sequence {U t } t∈N is bounded.

Thus it admits at least one limit point U * . Following the same argument as [30, Theorem 1] we have

g (U * |U * ) ≤ g (U|U * ) , ∀U ∈ St(p, k), (19)
meaning that U * is a global minimizer of the problem minimize

U∈St(p,k) g (U|U * ) . ( 20
)
Since the linear independence constraint qualification (LICQ) holds at U * [7, Lemma 3], it satisfies the following conditions:

∇g (U * |U * ) = U * Λ (U * ) H U * = I, (21)
where Λ is the associated Lagrange multiplier (see Appendix of [START_REF] Breloy | Clutter subspace estimation in low rank heterogeneous noise context[END_REF] for more details). The continuity of f and g together with the upperbound condition of g implies that ∇g (U * |U * ) = ∇f (U * ), which, combined with [START_REF] Abrudan | Steepest descent algorithms for optimization under unitary matrix constraint[END_REF], implies that U * satisfies the KKT conditions of problem [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF]. Now, assume that the whole sequence {U t } ∞ t=0 does not converge to K. Then there exists a convergent subsequence, indexed by {t j }, such that lim j→∞ d (tj ) (K) ≥ c for c ∈ R * + , and where

d (tj ) (K) = min Y∈St(p,k) ||U tj -Y||. The subsequence {U tj } ∞
j=0 is bounded, and the compactness of the set St(p, k) again implies that it admits at least one limit point. As shown previously, this limit point is a KKT point of problem [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF], which is a contradiction. Thus, the whole sequence {U t } ∞ t=0 converges to K.

Remark 3. Note that the convergence to K does not imply the convergence of Algorithm 2 in terms of the variable U. Establishing this property requires a case-by-case analysis which goes beyond the scope of this paper. In some cases the monotonic decrement of the objective can directly imply the convergence in terms of variable [START_REF] Kiers | Maximization of sums of quotients of quadratic forms and some generalizations[END_REF]. For the case of rotation invariant costs, this convergence in variable requires to be expressed in terms of subspace, e.g., as in [START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF].

C. Computational cost

First, recall that p denotes the dimension of the data, k refers to the dimension of the subspace of interest, and t indexes the algorithm steps. The sample size (number of columns of the data matrix) will be denoted n. The computational cost will be studied per iteration, therefore, it is to be multiplied by the total number of iterations T to obtain the overall complexity. A single iteration in Algortithm 2 essentially involves two operations:

• The computation of the matrix R(U): this step usually involves functions of the p × n data matrix and/or multiplying this matrix with the current point U t . Thus, this step is generally O(npk) (cf. examples in Section V). Also notice that this computation can, most of the time, be parallelized. Hence it does not represent the major bottleneck of Algorithm 2, contrarily to the second step. • The computation of P St : this step requires to compute the thin-SVD of a tall matrix R ∈ C p×k which is O(pk 2 + k 3 ). Also notice that this projection can be obtained through

P St (R) = R(R H R) -1/2 , (22)
which does not improve the theoretical complexity but can be advantageous to implement in some practical cases. Comparing to the existing approaches, e.g., the steepest descent on the Stiefel manifold [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], an iteration requires computing the gradient (also generally O(npk)) and a retraction (local mapping between a point in St(p, k) and its tangent space). The choice of the retraction is not unique, which leads to several options, e.g., based on geodesic paths [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF], Euclidean projection [START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF], or QR decomposition [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Nevertheless, for all of the corresponding algorithms, the retraction step is O(pk 2 + k 3 ). Hence, the computational complexity of an iteration of Algortithm 2 is on par with standard first-order based methods. However, this MM procedure is step-size free, thus it does not require the knowledge of any global parameter (such as the Lipschitz constant), or its adaptive estimation using a line search-type method. Compared to the latter option, this property effectively reduces the computational burden of each iteration, as it does not involve multiple computations of the retraction step (the computational bottleneck).

Establishing the total complexity of Algorithm 2 implies to study its convergence rate, which would require additional assumptions on the objective function and goes beyond the scope of this paper (pointers for some special cases are discussed in Section V-D). The convergence speed of an MM algorithm is generally expected to be sub-linear. Nevertheless, this possibly slow convergence can be compensated by the low-complexity of each iteration, resulting in a faster algorithm when considering the total computation time. More details on these remarks are illustrated with experiments in Section VII.

D. Extensions to block-MM

In addition to U ∈ St(p, k) some problems may involve a set of side parameters θ ∈ C p (cf., e.g., [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF][START_REF] Uematsu | SOFAR: Largescale association network learning[END_REF][START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF]), and be formulated as

minimize U,A f (U, A) subject to (U, A) ∈ St(p, k) × Θ, ( 23
)
where Θ is a compact of C p . Such problems can be tackled by partitioning the variables and applying the block-MM algorithm [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF], i.e., performing cyclic updates of the blocks, while keeping the others fixed. In our context, U typically represents a block updated by using one or several iterations of Algorithm 2. In this case, the convergence of the block-MM algorithm can be established from the analysis of the BSUM algorithm [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF][START_REF] Hong | A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF] in conjunction with Proposition 1.

IV. STANDARD COST FUNCTIONS AND THEIR SURROGATES

The key to apply Algorithm 2 is to obtain a linear surrogate of the objective on the set St(p, k). In this section, we derive such surrogates functions for several standard minimization problems. Furthermore, maximization problems can also be tackled by considering the negative of the objective function (or equivalently, Minorization-Maximization).

A. Quadratic forms

First, define the Brockett function [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]Sec. 4.8] for

U = [u 1 | • • • |u k ] ∈ St(p, k) as f B (U) = k r=1 d r u H r Mu r = Tr U H MUD , (24)
with M ∈ H + p , and D ∈ D k a diagonal matrix with

[D] r,r = d r satisfying 0 ≤ d 1 ≤ • • • ≤ d k .
In the following, functions of the form f B (resp. -f B ) are referred to as convex (resp. concave) QFs. Note that some other expressions of QFs exist, but they can usually be rewritten as special cases or combinations (e.g., sums) of Brockett functions. Proposition 2. (Majorization of concave QF) The function -f B as in (24) admits at point U t R a linear majorizing surrogate in the form of (18), with

R(U t ) = MU t D. ( 25
)
Equality holds at U t .

Proof. The function -f B as in [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF] is concave, so it can be majorized at point U t by its first order Taylor expansion (cf.

[24] section III.A), i.e.,

-f B (U) ≤ -Tr MU t D H U -Tr U H MU t D + const.. (26)
Remark 4. Majorizing a convex QF of U by a linear one seems counter-intuitive since it is not possible on the entire Euclidean space C N ×R . Nevertheless, the restriction to the set St(p, k) will make such an upperbound possible in Proposition 3. In order to give some insight, a visual example on R 2 is also presented in Figure 1.

Proposition 3. (Majorization of convex QF)

The function f B in (24) admits on St(p, k) and at point U t a linear majorizing surrogate in the form of (18), with

R(U t ) = -KU t D, (27)
where K = Mλ max M I and λ max M is the largest eigenvalue of M. Equality holds at U t .

Proof. The function f B in (24) can be expressed as

Tr U H MUD = Tr U H (M -λ max M I) UD + Tr U H (λ max M I) UD , (28)
where the second term is constant and equal to λ max M Tr {D} for the restriction U ∈ St(p, k). The first term of this expression is concave in U (U ∈ C M ×R ) so it can be upperbounded by its first order Taylor expansion, thus

f B (U) ≤ + Tr KU t D H U + Tr U H KU t D + const., (29)
with K defined as in Proposition 3.

B. Concave compositions of quadratic forms

Compositions involving inner QFs that yield concave functions are often used in order to build robust loss functions (examples are given in section V-C). The following proposition gives a linear majorizer of concave functions composed from the Brockett function. [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF], the function ρ(-f B ) admits at point U t a linear majorizing surrogate in the form of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], with

R(U t ) = ρ -f B U t MU t D. ( 30
)
Equality holds at U t .

Proof. The function -f B is concave in U and ρ is concave non-decreasing. It follows that the function ρ(-f B ) is concave [START_REF] Boyd | Convex Optimization[END_REF], so it can be upper-bounded at point U t by its first order Taylor expansion, i.e.,

ρ (-f B (U)) ≤ -ρ -f B U t Tr MU t D H U -ρ -f B U t Tr U H MU t D + const. (31)
Following this proof, other linear surrogates can be derived using compositions of concave (non-decreasing/nonincreasing) functions and the chain rule. It is also worth noting that one can apply the reformulation of Proposition 3 to express a quadratic QF as a concave term plus a constant in order to do so (the obtained majoration is then only valid on St(p, k)).

C. Quotients of quadratic forms

Various formulations of quotients of quadratic forms arise in generalized versions of PCA [START_REF] Kiers | Maximization of sums of quotients of quadratic forms and some generalizations[END_REF]. Most of them can be obtained as linear combinations of functions of the form

f q (U) = -Tr U H CU -1 U H AU , (32)
where C is a positive definite and A is positive semi-definite.

Proposition 5. (Majorization of quotient of QFs) The function f q as in (32) admits on St(p, k) and at point U t a linear majorizing surrogate in the form of (18), with

R(U t ) = T(U t ) -KU t T(U t ) , (33)
and

T(U t ) = AU t (U t ) H CU t -1 , T(U t ) = A -1/2 T(U t ) H A -1/2 T(U t ) , K = C -λ max C I, (34)
where λ max

C

is the largest eigenvalue of C. Equality holds at U t .

Proof. Starting from the inequality

U H CU -1/2 U H A 1/2 -U H CU 1/2 (U t ) H CU t -1 (U t ) H A 1/2 2 ≥ 0, (35)
we obtain

f q (U) ≤ -2Re T(U t ) H U + Tr U H CU T(U t ) , (36)
with T(U t ) and T(U t ) as in [START_REF] Breloy | Clutter subspace estimation in low rank heterogeneous noise context[END_REF], and where equality holds at U t . Following the proof of Proposition 3, we can majorize on St(p, k) the quadratic term in [START_REF] Hong | A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF] as

Tr U H CU T(U t ) ≤ + Tr KU t T(U t ) H U + Tr U H KU t T(U t ) + const., (37)
with K as in [START_REF] Breloy | Clutter subspace estimation in low rank heterogeneous noise context[END_REF], and where equality holds again at U t . Combining the inequalities (36) and (37) concludes the proof.

V. EXAMPLES AND LINKS WITH EXISTING ALGORITHMS

This section provides some examples of MM algorithms that leverage the proposed method with the surrogate functions from Section IV. When existing, some links with the existing state-of-the-art algorithms are highlighted.

A. Power iteration methods

Let u ∈ St(p, 1) and M ∈ H + p . Consider the classical problem maximize

u u H Mu subject to u H u = 1, (38)
whose local solutions are the eigenvectors of M. In order to compute those solutions, we can apply the framework proposed in Section III. First, we simply rewrite the problem as a minimization one by changing the sign of the objective:

minimize u∈St(p,1) u H (-M)u. (39)
Second, since the objective in ( 39) is a concave QF, we can apply Proposition 2 in Section IV-A to obtain a majorizing linear surrogate function at the point u t as:

g(u|u t ) = -u H Mu t -(u t ) H Mu + const. (40)
Finally, we recognize a surrogate function as in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], with R (u t ) = Mu t . Therefore, we can apply Algorithm 2, which leads to the following MM iterations:

u t+1 = P Proc Mu t . (41)
This algorithm simply corresponds to the normalization of the iterate Mu t . Thus, Algorithm 2 for this problem yields the well-known power iteration method [START_REF] Golub | Matrix Computations[END_REF]. This algorithm also appears as special case of a steepest-descent method on the Stiefel manifold (cf. Section 4.6.6 [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]).

Let us now consider the general case U ∈ St(p, k), with corresponding problem:

maximize U U H MU subject to U H U = I. ( 42
)
The stationary points of this problem are k-dimensional invariant subspaces of M. Applying the same steps as in ( 39)-( 41), we obtain the MM algorithm:

U t+1 = P Proc MU t , (43)
which is a generalization of the power method similar to the orthogonal iterations [START_REF] Golub | Matrix Computations[END_REF] (that equivalently relies on the QR decomposition rather than P Proc to orthonormalize the updates).

B. MM for generic non-homogeneous QF (with a linear term)

Let U = [u 1 | • • • |u k ] ∈ St(p, k
) and a set of k positive semi-definite Hermitian matrices {M r }. Consider the following problem (appearing for example as a subproblems in [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF][START_REF] Ben Abdallah | Bayesian signal subspace estimation with compound gaussian sources[END_REF]):

maximize U k r=1 u H r M r u r + 2Re{u H r c r } subject to u H i u j = δ i,j ∀ i, j ∈ [[1, k]], (44)
which has no closed form solution related to the SVD of the matrices M r and the vectors c r . First, we rewrite the problem as

minimize U∈St(p,k) k r=1 u H r (-M r )u r -2Re{u H r c r } . ( 45
)
The objective of ( 45) can be split into two functions and allows us to follow ( 15)- [START_REF] Vidal | Subspace clustering[END_REF]. First, notice that the linear elements of the sum are already expressed as in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] since

-2Re{u H r c r } = -u H r c r -c H r u r , ∀r ∈ [[1, k]]. (46)
Second, for each concave quadratic element of the sum, we apply the same steps as in ( 39)- [START_REF] Jaggi | Revisiting frank-wolfe: Projection-free sparse convex optimization[END_REF]. This leads to a surrogate function for the total objective in [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF] in the form of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] with

R U t = M 1 u t 1 + c 1 | • • • |M k u t k + c k . (47)
Therefore, we can apply Algorithm 2, which leads to the following MM iterations:

U t+1 = P Proc M 1 u t 1 + c 1 | • • • |M k u t k + c k . ( 48
)
We also note that this algorithm corresponds to the one of [START_REF] Koschat | A weighted Procrustes criterion[END_REF], which is often called for minimizing linear-plus-quadratic terms under orthonormality constraints.

C. MM for nonconvex robust subspace recovery (RSR)

The aim of RSR is to estimate a low-dimensional subspace from a (demeaned) dataset {z i } n i=1 , while being robust to potential outliers within this set. There are many approaches to tackle this problem [START_REF] Lerman | An overview of robust subspace recovery[END_REF]. Among them, several RSR estimators have been proposed as minimizers of a nonconvex robust regression function [START_REF] La | A framework for robust subspace learning[END_REF][START_REF] Maronna | Principal components and orthogonal regression based on robust scales[END_REF][START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF][START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF][START_REF] Maunu | A well-tempered landscape for non-convex robust subspace recovery[END_REF], with orthonormality constraint.

Following [START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF], a problem for such RSR can be formulated as follows:

minimize U∈St(p,k) n i=1 ρ d 2 (U, z i ) , (49)
where

d 2 (U, z) = (UU H ) ⊥ z 2 F = z H z -Tr U H zz H U (50)
is the Euclidean distance between a vector z ∈ C p and the subspace spanned by U ∈ St(p, k), and ρ : R → R is a function that ensures the robustness to outliers. Here, ρ is assumed to be a concave nondecreasing function. However, this constraint is not restrictive, as it holds for a wide variety of usual robust formulations, as illustrated by the following examples:

Example 1. ( p -norm) For p > 0, p -norm nonconvex RSR estimators are defined as in ( 49) by using the function

ρ p (t) = t p/2 . ( 51
)
The least-square estimator is defined for p = 2. In this case, it is not hard to show that ( 49) is equivalent to [START_REF] Weber | Nonconvex stochastic optimization on manifolds via riemannian frank-wolfe methods[END_REF] with

M = n k=1 z k z H k .
The solution is the k the leading eigenvectors of the sample covariance matrix, which corresponds to the standard PCA estimator.

Example 2. (Huber-type) For a parameter T ≥ 0, Huber-type nonconvex RSR estimators are defined as in (49) by using the function

ρ H (t) = t/ √ T if t ≤ T, 2 √ t - √ T if t > T. ( 52
)
The median estimator, e.g., considered in [START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF][START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF], corresponds to the case limit case T → 0.

ρ GMC (t) = t/(T + t), (54)
which has been used in, e.g., [START_REF] La | A framework for robust subspace learning[END_REF].

From the problem in ( 49), we exhibit some links between algorithms from the state of the art and the proposed framework. First, notice that ρ is a concave nondecreasing and d in ( 50) is a concave QF of U. Therefore each element of the sum of the objective in ( 49) can be majorized by using Proposition 4 in Section IV. This leads to a surrogate function in the form of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] with R (U t ) = M (U t ) U t , where

M (U) = n i=1 ρ d 2 (U, z i ) z i z H i . (55)
Applying Algorithm 2 thus yields the MM iterations

U t+1 = P Proc M U t U t . (56)
The resulting algorithm corresponds to R1-PCA-type fixedpoint iterations, initially proposed as heuristic in [START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF].

Additionally, the MM perspective allows us to draw some links with other algorithms that aim at solving [START_REF] Mangasarian | Machine learning via polyhedral concave minimization[END_REF]. It is shown in [START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF] that the objective of ( 49) is also majorized with equality at the point U t by a concave quadratic surrogate function as

n i=1 ρ d 2 (U, z i ) ≤ -Tr U H M U t U + const. (57)
Minimizing the surrogate function in (57) is a problem equivalent to the one in [START_REF] Weber | Nonconvex stochastic optimization on manifolds via riemannian frank-wolfe methods[END_REF] (with M = M (U t )). An update can thus be obtained by computing the k leading eigenvectors of M (U t ). This algorithm corresponds to the fixed-point iteration proposed in [START_REF] Maronna | Principal components and orthogonal regression based on robust scales[END_REF]. An equivalent update can be obtained through the SVD of the data matrix

Z t = [z 1 | • • • |z n ] using normalized samples z t i = ρ (d 2 (U t , z i ))z i
, which corresponds to the FMS algorithm proposed in [START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF]. Following the example of Section V-A, a third option to compute this update is to perform an inner loop using an MM algorithm (i.e., the power method in [START_REF] Journée | Generalized power method for sparse principal component analysis[END_REF]). This procedure corresponds to the algorithm referred to as FMS with power method in [START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF]. Alternatively, the surrogate function in (57) is a concave QF so it can be majorized with Proposition 2: the resulting MM algorithm is again (56). Within this perspective, we can also note that Algorithm (2) in (56) (R1-PCA-type) also appears as the FMS with power method when only one power iteration is performed at each step. Some simulations comparing these approaches are presented in Section VII. Compute s = argmin s∈D s, ∇f (x t )

4:

Update x t+1 = (1γ)x t + γs for γ = 2/(t + 2)

5:

t = t + 1 6: until convergence criterion is met

D. Frank-Wolfe and first-order approximation methods

The Frank-Wolfe method [START_REF] Frank | An algorithm for quadratic programming[END_REF] (cf. Algorithm 3) regained popularity in machine learning related fields [START_REF] Jaggi | Revisiting frank-wolfe: Projection-free sparse convex optimization[END_REF] for minimizing a convex function f (with gradient ∇f ) over a convex set D. This approach was also recently extended to the Riemannian optimiation framework in [START_REF] Weber | Nonconvex stochastic optimization on manifolds via riemannian frank-wolfe methods[END_REF] as an alternative to first order projection based method (e.g., Steepest descent [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]). In this method, a descent direction is obtained by minimizing a first-order linear approximation of the objective, which is reminescent of the proposed approach (though the linear approximation is not a majorizer of the concave function). However, the problems considered in this paper are not convex (nor g-convex) so Frank-Wolfe method is not trivially transposable on the Stiefel manifold. Yet, similar firstorder approximation methods have been proposed for this set. Notably, [START_REF] Journée | Generalized power method for sparse principal component analysis[END_REF]Sec. 3.4] considered minimizing a concave function with the following algorithm

U t+1 = argmin U∈St(p,k) U, ∇f (U t ) (58)
which can be interpreted as a Frank-Wolfe method with a constant (maximal) step γ = 1. This algorithm also corresponds to an MM approach since the first order Taylor approximation is a linear majorizing surrogate of a concave function. Interestingly, the proposed framework also generalizes to objectives that are not necessarily concave, as it only requires a linear surrogate of the objective on St(p, k). Conversely, some stronger convergence results (such as the convergence rate) for Algorithm 2 can be obtained from [START_REF] Journée | Generalized power method for sparse principal component analysis[END_REF] when the objective is concave.

VI. MAJORIZATION-MINIMIZATION FOR SPARSE PCA

In standard PCA, the estimated principal components are usually dense (i.e., a linear combinations of all entries of the variables). Since the principal components have an actual physical meaning in many applications, estimating sparse principal components can significantly help the interpretation, as well as the variable selection process. Many algorithms have been proposed to perform this task [START_REF] Chen | Reduced rank stochastic regression with a sparse singular value decomposition[END_REF][START_REF] Chen | Sparse reduced-rank regression for simultaneous dimension reduction and variable selection[END_REF][START_REF] Bunea | Joint variable and rank selection for parsimonious estimation of high-dimensional matrices[END_REF][START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF][START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF][START_REF] Uematsu | SOFAR: Largescale association network learning[END_REF]. Most of the proposed methods involve a variable U ∈ St(p, k), whose columns represent the principal components, and can be generically formulated through the problem minimize

U∈St(p,k) L(U, {z i } n i=1 ) + λξ (U) , (59)
where L is a data fitting term (orthogonal regression on the dataset {z i } n i=1 ), ξ is a sparsity promoting penalty, and λ is a regularization parameter. The introduction of the penalty usually makes the minimization in (59) hard to deal with under the orthonormality constraint. Thus, most algorithms in the literature relax this constraint and resort to a trade-off [START_REF] Chen | Reduced rank stochastic regression with a sparse singular value decomposition[END_REF][START_REF] Chen | Sparse reduced-rank regression for simultaneous dimension reduction and variable selection[END_REF][START_REF] Bunea | Joint variable and rank selection for parsimonious estimation of high-dimensional matrices[END_REF][START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF]. In the following, this issue is addressed by using the proxies of 0 -norm proposed in [START_REF] Song | Sparse generalized eigenvalue problem via smooth optimization[END_REF] to promote sparsity, while leveraging the MM framework of Section III to perform the optimization under the orthonormality constraint.

A. Sparse regularizers with linear surrogates on St(p, k)

One approach to force sparsity in the principal components

U = [u 1 | • • • |u k ]
is to use a penalty involving the 0 -norm3 :

U 0 = r,i sgn ( |[u r ] i | ) , (60)
where |•| stands for the modulus of a complex number. Another approach, that promotes a row-sparse structure, is the use of the 2,0 -norm, defined by:

(U) 2 2,0 = p i=1 k r=1 sgn ( |[u r ] i | ) 2 . (61)
However, both of these functions are too complex to deal with due to their discontinuity.

To alleviate this issue, we follow the approach proposed in [START_REF] Song | Sparse generalized eigenvalue problem via smooth optimization[END_REF], i.e., approximating the sign function by a smooth function denoted l γ , and defined as

l γ (x) = a|x| 2 , if |x| ≤ l γ (x) -b, if |x| > , (62)
with appropriate constants a and b so that the approximations l γ are continuous and differentiable (cf. [START_REF] Song | Sparse generalized eigenvalue problem via smooth optimization[END_REF]), and where l γ belongs to the family of functions: a) γ -norm [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF][START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF][START_REF] Lai | Improved iteratively reweighted least squares for unconstrained smoothed q minimization[END_REF]: [START_REF] Sriperumbudur | A majorization-minimization approach to the sparse generalized eigenvalue problem[END_REF][START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF]:

l γ (x) = |x| γ , γ ∈ (0, 1] b) 1 -norm approximation from
l γ (x) = ln(1 + |x| /γ) ln(1 + 1/γ), γ > 0 c
) lower bound of sign function from [START_REF] Mangasarian | Machine learning via polyhedral concave minimization[END_REF]: l γ (x) = 1e -|x|/γ , γ > 0, involving a tuning parameter γ for each case. Thus, this class covers most of standard 1-dimensional sparsity forcing penalties (i.e., a proxy of the sign function). Notice that this class is still valid for data with complex entries by reading |•| as the modulus function. We have the following proposition (cf. [31, Section III]): Proposition 6. The function l γ in (62) is majorized at point x t by the following quadratic surrogate:

l γ (x|x t ) ≤ φ(x t )|x| 2 + const. ( 63
)
where the function φ depends on l γ as: a) γ -norm [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF][START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF][START_REF] Lai | Improved iteratively reweighted least squares for unconstrained smoothed q minimization[END_REF]: [START_REF] Sriperumbudur | A majorization-minimization approach to the sparse generalized eigenvalue problem[END_REF][START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF]:

φ(x t ) = (γ/2) γ-2 , |x t | ≤ (γ/2) |x t | γ/2 , |x t | > b) 1 -norm approximation
φ(x t ) = ( 2 (γ + ) ln(1 + 1/γ) ) -1 , |x t | ≤ ( 2 ln(1 + 1/γ) |x t | (|x t | + γ) ) -1 , |x t | > c
) lower bound of sign function [START_REF] Mangasarian | Machine learning via polyhedral concave minimization[END_REF]:

φ(x t ) = e -/γ /2γ , |x t | ≤ e -|xt|/γ /2γ |x t | , |x t | > Equality is achieved at x t .
The functions l γ will serve as basic building blocks to build sparsity promoting penalties on U. In the following we propose two functions that extend the approach of [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF], with corresponding surrogates functions in the form of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Our aim is to give insight on how to build meaningful cost functions that can be easily handled within the framework of Section III: the key trick is to obtain a series of inequalities, notably using Proposition 6, that leads to a quadratic surrogate of U over St(p, k). At this step, Proposition 3 can be applied to obtain a linear majorizer on St(p, k).

1) Proxies of the 0 -norm on St(p, k): In order to mimic a weighted 0 -norm in (60), we consider the cost function:

r 0 (U) = k r=1 ω r p i=1 l γ ([u r ] i ) , (64)
where ω r are positive weights and with l γ in (62). Such type of penalty was initially proposed for Sparse PCA in [START_REF] Benidis | Orthogonal sparse PCA and covariance estimation via Procrustes reformulation[END_REF] using the 1 -norm approximation [START_REF] Sriperumbudur | A majorization-minimization approach to the sparse generalized eigenvalue problem[END_REF][START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF] (i.e., l γ in the class b) of the considered family)).

Proposition 7. The function r 0 in (64) is majorized on the set of unitary matrices St(p, k) at point

U t = [u t 1 | • • • |u t k ] as r 0 (U) ≤ Tr H H U t U + Tr U H H U t + const., (65)
with matrices H(U t ) ∈ C p×k , W r ∈ C p×p , w r ∈ C p×1 defined as

         H U t = W 1 u t 1 | • • • | W k u t k W r = diag( w r -w max r 1 p ) [w r ] i = ω r φ u t r i w max r = max(w r ), (66)
with φ defined (depending on the chosen l γ ) in Proposition 6.

Proof. Proposition 6 allows to upperbound r 0 in (64) at point

U t as r 0 (U) ≤ k r=1 u H r (diag(w r )) u r + const., (67)
with w r defined in (66) and where equality is achieved at U t . Proposition 3 applied to each quadratic term of this majorizer leads to the surrogate in (65) and (66).

2) Proxies of the mixed-norm on St(p, k): To mimic the 2,0 -norm in (61), the main idea is to involve a nonlinear coupling among the elements in a row. We consider therefore the cost function:

r 2,0 (U) = p i=1 ln 1 + k r=1 l γ ([u r ] i ) . (68)
Notice that the function ln is essentially used in order to obtain a surrogate satisfying [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. However, the main goal that it exploits a tighter upperbound than the MM with linear majorization (cf. Section V-C). The speed of the steepest descents could be improved by the use of an exhaustive line search. Still, the Armijo rule offers a significant improvement compared to the constant step size. When it comes to the computation time, the MM with linear majorization appears faster. Indeed, each iteration of this algorithm requires only the computation of the thin-SVD of an p × k matrix, rather than the k strongest eigenvectors of a p × n matrix for the MM with quadratic majorization. The steepest descent the Armijo rule requires to run a line search for the step size. This involves evaluating the objective, therefore calls for multiple retractions, which is computationally more intense. The constant step size iterations are fast to compute, which can compensate its slow convergence in terms of number of iterations. Nevertheless, the step size had to be manually tuned for each scenario. Thus, this algorithm corresponds to a benchmark with little practical application Figure 3 displays total computation time, i.e., time to reach the stopping criterion ||U t -U t+1 || ∞ < 10 -10 , with respect to the problem dimensions (scaled proportionally to to p). As previously, we can notice that the MM with linear majorization offers an interesting solution in terms of scalability.

B. On RSPCA algorithm 1) Validation on synthetic data: In this section, we validate several properties of the RSPCA method on synthetic examples. In order to clarify the exposition, the term RSPCA will here denote the solution of problem (72) computed with Algorithm 4. We will focus on the use of the following parameters 4 : the RSR cost is ρ H (cf. Example 2) with T = 0.1, the sparse penalty is either r 0 or r 2,0 , where l γ is from the proxy of the 1 -norm (i.e., a)) in the familly of proposed proxies), with = 10 -2 and ω r = 1∀[[∈ 1, k]], and varying regularization parameter λ. the dimensions are p = 100, k = 15.

First, we study the performance of RSPCA in a standard PCA setting. Samples are drawn from z ∼ CN (0, SNR × UU H + I), where U are the first k eigenvectors of a the Topelitz matrix [Σ T ] i,j = (θ(1 + i)/ √ 2) |i-j| and SNR = 10. For a given estimator Û, the performance criterion is the average fraction of recovered energy:

AFE = E Tr{ ÛH UU H Û}/k , (78)
which is evaluated through 10 2 Monte-Carlo simulations. Figure 4 (resp. 5) displays the AFE of RSPCA with r 0 (resp. r 2,0 ) penalty in function of n and for various λ. Note that without regularization (λ = 0) the RSPCA corresponds to RSR with the Huber cost [START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF][START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF]. In a Gaussian setting, it is interesting to note that there is little loss due to the usage of the Huber cost function when comparing with the standard PCA, which corresponds to the maximum likelihood estimator in this case [START_REF] Jolliffe | Principal Component Analysis and Factor Analysis[END_REF]. If a sparse regularization penalty is involved, two regimes can be observed. When the underlying subspace basis is actually sparse (θ = 0 leads to canonical eigenvectors), this penalty significantly improves the recovery. When the underlying subspace basis is not sparse (θ = 0.5), the recovery can be severely impacted by a strong penalty promoting this property. However, we can notice that RSPCA appears robust to the mismatch when λ is set low enough. Second, we study the robustness of RSPCA to corruption by outliers in the sample set. The simulation setup is built around the so-called haystack model [START_REF] Lerman | Robust computation of linear models by convex relaxation[END_REF], which corresponds here to a mixture of orthogonal Gaussian distributions plus additive noise:

{z i } n i=1 = {{z in i } nin i=1 , {z out i } n i=nin+1 } z in ∼ CN (0, SNR × UU H + I) z out ∼ CN (0, ONR × U ⊥ U H ⊥ + I), (79)
where U and U ⊥ are built from the canonical basis such that U H U ⊥ = 0, SNR is signal to noise ratio, and ONR is outlier to noise ratio. We will compare four algorithms: i) RSPCA with least square fitting (cf. Example 1) and λ = 0 (i.e., the standard PCA); ii) RSPCA with least square fitting and λ = 100; iv) RSPCA with Huber cost (cf. Example 2) and λ = 0 (i.e., RSR); iii) RSPCA with Huber cost and λ = 1000.

The aim is to illustrate the contributions of both the robust fitting objective and the sparsity promoting penalty. Figure 6 displays the AFE of each algorithm with respect to ONR and the fraction of outliers in the sample set. We can notice that the use of a robust cost improves the performance compared to the standard PCA. Moreover, the introduction of the sparse penalty improves the results in terms of AFE (as seen in Figures 45), but also interestingly improves the robustness of the estimation process. Here, it is worth mentioning two critical points: a) RSPCA appears very robust when a valid signal subspace basis is actually sparse, which is probably because the sparse penalty contributes to naturally discard dense outliers. If the true subspace basis is dense, results can be degraded in practice (as also observed in Figures 45). b) The starting point plays an important role in the achieved robustness of all iterative algorithms (even without regularization). In these simulations, we used the spherical PCA (PCA applied on the normalized samples) as starting point. It has been observed that the achieved robustness can be lowered by using the standard PCA instead.

2) Experiment on real data: In this experiment, we compare the performance of RSPCA with state of the art sparse PCA algorithms on the Leukemia data set [START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF] 5 . RSPCA build with a GMC cost (cf. Example 4), r 0 penalty and l γ is from the proxy of the lower bound of the sign function (i.e., c) in the familly of proposed proxies). Again, similar conclusions can be drawn with other objectives, up to minor changes of the parameters. RSPCA is computed with Algorithm 4 using an outer loop, decreasing from 10 -1 to 10 -7 in order to avoid potential local minimums. This algorithm is compared to ALSPCA [START_REF] Lu | An augmented lagrangian approach for sparse principal component analysis[END_REF], SPCArt and rSVD-GP from [START_REF] Hu | Sparse principal component analysis via rotation and truncation[END_REF]. The studied which measure respectively the sparsity (expected to be close to 1, i.e., 100%), the explained variance (expected to be close to 1, i.e., 100%), and the non-orthonormality (expected to be low). Notice that the AFE criterion slightly favors the algorithms that relax the orthonormality constraint in this case. Figure 7 displays the AFE and NOR versus SP for the studied sparse PCA algorithms on the Leukemia data set. Interestingly, we can notice that RSPCA achieves state of the art performance when it comes to the explained variancesparsity trade-off, but without relaxing the orthonormality constraint, as done by the other algorithms.

VIII. CONCLUSION

This paper presented an optimization framework to deal with orthonormality constraints based on Majorization-Minimization. The core idea is to systematically obtain linear surrogates in order to formulate iterations as Euclidean projection problems. Such majorizers can be obtained for a large family of functions for which we presented a catalog. We also exhibited some links between this approach and wellknown algorithms. Finally, the presented framework drove the formulation of a novel sparse PCA algorithm (namely, RSPCA), combining a robust subspace recovery cost and sparsity promoting penalties. Several experiments illustrated the interest of the approach, both in terms of computational and estimation performance.

  where only the k first column vectors of U from the SVD are used. The set of p × p Hermitian positive semi-definite matrices is denoted H + p . The set of p×p diagonal matrices is denoted D p .

Algorithm 1 1 : 3 : 4 :Algorithm 2

 11342 Computation of P St (projection on St(p, k)) Entry: R ∈ C p×k 2: Compute the thin-SVD: R TSVD = V left DV H right Set U=V left V H right Output: U = P St (R) ∈ St(p, k) Generic MM Algorithm for St(p, k)

Fig. 1 .

 1 Fig. 1. Quadratic form on R 2 and its restriction to St(2, 1) (up). Linear majorization of this quadratic form at a given point over the set St(2, 1) (down).

Proposition 4 .

 4 (Majorization of concave function composed from concave QF) Let ρ : R → R be a concave nondecreasing function. For f B as in

Example 3 .

 3 (Cauchy-Lorentz-type) For a parameter T ≥ 1, Cauchy-Lorentz-type nonconvex RSR estimators are defined as in (49) by using the function ρ CL (t) = T ln(T + t). (53) Example 4. (Geman-McClure-type) For a parameter T ≥ 0, Geman-McClure-type nonconvex RSR estimators are defined as in (49) by using the function

Algorithm 3

 3 Frank-Wolfe method 1: initialize: x 0 ∈ D 2: repeat 3:

Fig. 2 .pFig. 3 . 2 F

 232 Fig.2. Gap to minimal value (attained by any converging algorithm) versus time and number of iterations for various algorithms called on problem (49): Steepest descent using Euclidean projection retraction with constant step size ) or Armijo rule ( * ). Steepest descent using QR retraction with Armijo rule (+). MM with quadratic majorizer (•). MM with linear majorizer ( ). Simulation carried out on MATLAB R2019a with Intel(R) Core(TM) i7-8850H CPU 2.60GHz processor.

Fig. 4 .

 4 Fig. 4. AFE n for various algorithms in a standard PCA context. p = 100, k = 10, SNR = 10. RPSCA built with r 0 penalty.

Fig. 5 .Fig. 6 .

 56 Fig. 5. AFE versus n for various algorithms in a standard PCA context. p = 100, k = 10, SNR = 10. RPSCA built with r 2,0 penalty.

Fig. 7 .

 7 Fig. 7. AFE and NOR versus SP for various sparse PCA algorithms on the Leukemia data set. p = 7129, k = 10, n = 72.

Or semi-orthogonal in the real valued case.

Rank deficiency of this matrix is a case we do not focus on: some pathological counter-examples can be build but they rely on either i) a cost function f that does not satisfy the initial regularity assumptions; ii) a subspace within U t that has reached a local stationary point. For the second point, the proposed method can still be applied by setting the stable subspace fixed and updating only the remaining portion of U t (i.e., recasting the problem with k < k). In practice, the issue has not been experienced with the cost functions (and starting points) considered in this paper.

Note that this function is not a proper norm, but we still adopt the common abuse of terminology.

We simply note that the general conclusions drawn from these examples can also be obtained with the other robust-fitting/sparse-promoting function presented in this paper, up to a good selection of the tuning parameters

From https://github.com/ramhiser/datamicroarray/wiki/[START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF]

Algorithm 4 MM for RSPCA 1: initialize t = 0, U 0 2: repeat 3:

Compute M (U t ) with equation (55).

4:

Compute the regularizer K (U t ) with (74)

with φ defined (depending on the chosen l γ ) in Proposition 6.

Proof. The function ln is concave, so it can be upperbounded by its first order Taylor expansion in (68), giving:

with ω t n in (70). The rest of the proof follows similarly to the one of Proposition 7.

B. Robust Sparse PCA (RSPCA): Proposed Algorithm 1) Formulation: From the generic formulation of (59) sparse PCA, we consider the following problem:

where the data fitting term is a nonconvex RSR problem (cf. Section V-C), and where ξ is either r 0 in (64) or r 2,0 in (68). Hence, the proposed objective offers a quite modular formulation that encompasses various RSR methods, and all the family of regularization penalty from Section VI-A.

2) Algorithm derivation: In this section, an MM algorithm is derived to aim at solving problem (72) by following the framework proposed in Section III. To majorize the nonconvex RSR cost, we apply the same steps as in Section V-C. The penalty ξ can be majorized according to either Proposition 7 or 8. Thus, the objective in (72) can be majorized at point U t by a surrogate satisfying [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] with

where M is defined in (55) and

Therefore, we can apply Algorithm 2, which leads to the following MM iterations:

This algorithm is summed up in Algorithm 4 and will be referred to as Robust-Sparse PCA (RSPCA).

Remark 5. Notice that we focused on the use of the RSR functions from Section V-C in order leverage their robustness properties. However, the proposed method can be generalized to any data fitting term (e.g., more complex likelihood functions) for which Algorithm 2 applies.

VII. SIMULATIONS

A. On MM algorithms with linear majorizers for St(p, k) This section illustrates the computational complexity advantage of the proposed Algorithm 2 in comparison to existing state-of-the-art techniques. As an example, we consider minimizing the RSR cost in [START_REF] Mangasarian | Machine learning via polyhedral concave minimization[END_REF] with the Huber-type function [START_REF] Lu | An augmented lagrangian approach for sparse principal component analysis[END_REF] with parameter T = 0.1. The data set is generated using a spiked Gaussian model, i.e., z ∼ CN (0, SNR×U 0 U H 0 +I), with SNR = 10 and where U 0 ∈ St(p, k) contains the first k vectors of the canonical basis (this setting is chosen as any rotation of the samples did not affect results in terms of computational speed). The following methods are compared: i) Steepest descent on the Stiefel manifold: for this algorithm, we compare both the retraction using the Euclidean projection [START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF]:

and the retraction through the Q factor of the the QR decomposition [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] (denoted qf{•}):

where γ t is a step size, and gradf | U is the Riemmanian gradient of f at U. The step size is computed at each iteration with the standard Armijo rule. We also compare an implementation of [START_REF] Manton | Optimization algorithms exploiting unitary constraints[END_REF] with a constant step size (selected as large as possible, while ensuring convergence); ii) MM with quadratic majorization: as discussed in section V-C, an iteration of this algorithm is solved either with the SVD of M(U t ) in (55) (fixed point [START_REF] Maronna | Principal components and orthogonal regression based on robust scales[END_REF]), or the SVD of appropriately normalized samples (FMS [START_REF] Lerman | Fast, robust and non-convex subspace recovery[END_REF]). In the considered scenarios n < p so the second option is chosen, as it is more computationally efficient; iii) MM with linear majorization (Algorithm 2): this algorithm is detailed in (56) and also corresponds to the fixed point iterations referred to as R1-PCA in [START_REF] Ding | R1-PCA: rotational invariant L1norm principal component analysis for robust subspace factorization[END_REF]. Figure 2 displays the convergence of the algorithms in terms of objective value for different scenarios with varying k (p = 1000 and n = p/2). The overall performance is compared both in terms of the number of iterations and the computation time. In terms of required number of iterations, the MM with quadratic majorization appears faster. This is due to the fact