Search

This paper considers the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on the algebraic variety of all $m$-by-$n$ real matrices of rank at most $r$. Several definitions of stationarity exist for this nonconvex problem. Among them, Bouligand stationarity is the strongest necessary condition for local optimality. Only a handful of algorithms generate a sequence in the variety whose accumulation points are provably Bouligand stationary. Among them, the most parsimonious with (truncated) singular value decompositions (SVDs) or eigenvalue decompositions can still require a truncated SVD of a matrix whose rank can be as large as $\min\{m, n\}-r+1$ if the gradient does not have low rank, which is computationally prohibitive in the typical case where $r \ll \min\{m, n\}$. This paper proposes a first-order algorithm that generates a sequence in the variety whose accumulation points are Bouligand stationary while requiring SVDs of matrices whose smaller dimension is always at most $r$. A standard measure of Bouligand stationarity converges to zero along the bounded subsequences at a rate at least $O(1/\sqrt{i+1})$, where $i$ is the iteration counter. Furthermore, a rank-increasing scheme based on the proposed algorithm is presented, which can be of interest if the parameter $r$ is potentially overestimated.

Optimization over the set of matrices $X$ that satisfy $X^\top B X = I_p$, referred to as the generalized Stiefel manifold, appears in many applications involving sampled covariance matrices such as the canonical correlation analysis (CCA), independent component analysis (ICA), and the generalized eigenvalue problem (GEVP). Solving these problems is typically done by iterative methods that require a fully formed $B$. We propose a cheap stochastic iterative method that solves the optimization problem while having access only to random estimates of $B$. Our method does not enforce the constraint in every iteration; instead, it produces iterations that converge to critical points on the generalized Stiefel manifold defined in expectation. The method has lower per-iteration cost, requires only matrix multiplications, and has the same convergence rates as its Riemannian optimization counterparts that require the full matrix $B$. Experiments demonstrate its effectiveness in various machine learning applications involving generalized orthogonality constraints, including CCA, ICA, and the GEVP., Comment: This v3 is a corrected version of the ICML 2024 paper (PMLR 235:49226-49248); see the errata at the end

We exhibit conjugate points on the Stiefel manifold endowed with any member of the family of Riemannian metrics introduced by H\"uper et al. (2021). This family contains the well-known canonical and Euclidean metrics. An upper bound on the injectivity radius of the Stiefel manifold in the considered metric is then obtained as the minimum between the length of the geodesic along which the points are conjugate and the length of certain geodesic loops. Numerical experiments support the conjecture that the obtained upper bound is in fact equal to the injectivity radius., Comment: v2: fixed MSC codes; fixed typo in the penultimate sentence of the proof of Theorem 6.1; added section 9 (Concluding remarks)

We propose a numerical method to obtain an adequate value for the upper bound on the rank for the tensor completion problem on the variety of third-order tensors of bounded tensor-train rank. The method is inspired by the parametrization of the tangent cone derived by Kutschan (2018). A proof of the adequacy of the upper bound for a related low-rank tensor approximation problem is given and an estimated rank is defined to extend the result to the low-rank tensor completion problem. Some experiments on synthetic data illustrate the approach and show that the method is very robust, e.g., to noise on the data.

In recent years, the proximal gradient method and its variants have been generalized to Riemannian manifolds for solving optimization problems with an additively separable structure, i.e., $f + h$, where $f$ is continuously differentiable, and $h$ may be nonsmooth but convex with computationally reasonable proximal mapping. In this paper, we generalize the proximal Newton method to embedded submanifolds for solving the type of problem with $h(x) = \mu \|x\|_1$. The generalization relies on the Weingarten and semismooth analysis. It is shown that the Riemannian proximal Newton method has a local quadratic convergence rate under certain reasonable assumptions. Moreover, a hybrid version is given by concatenating a Riemannian proximal gradient method and the Riemannian proximal Newton method. It is shown that if the switch parameter is chosen appropriately, then the hybrid method converges globally and also has a local quadratic convergence rate. Numerical experiments on random and synthetic data are used to demonstrate the performance of the proposed methods., Comment: Updates compared to the previous version: *) the updates for the published version. Additionally updates: *) local quadratic convergence rate *) update the proofs of Lemma 3.3. *) update the proofs of Proposition 3.13

Orthogonality constraints naturally appear in many machine learning problems, from principal component analysis to robust neural network training. They are usually solved using Riemannian optimization algorithms, which minimize the objective function while enforcing the constraint. However, enforcing the orthogonality constraint can be the most time-consuming operation in such algorithms. Recently, Ablin & Peyr\'e (2022) proposed the landing algorithm, a method with cheap iterations that does not enforce the orthogonality constraints but is attracted towards the manifold in a smooth manner. This article provides new practical and theoretical developments for the landing algorithm. First, the method is extended to the Stiefel manifold, the set of rectangular orthogonal matrices. We also consider stochastic and variance reduction algorithms when the cost function is an average of many functions. We demonstrate that all these methods have the same rate of convergence as their Riemannian counterparts that exactly enforce the constraint, and converge to the manifold. Finally, our experiments demonstrate the promise of our approach to an array of machine-learning problems that involve orthogonality constraints.

We consider the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on a stratified set and present a first-order algorithm designed to find a stationary point of that problem. Our assumptions on the stratified set are satisfied notably by the determinantal variety (i.e., matrices of bounded rank), its intersection with the cone of positive-semidefinite matrices, and the set of nonnegative sparse vectors. The iteration map of the proposed algorithm applies a step of projected-projected gradient descent with backtracking line search, as proposed by Schneider and Uschmajew (2015), to its input but also to a projection of the input onto each of the lower strata to which it is considered close, and outputs a point among those thereby produced that maximally reduces the cost function. Under our assumptions on the stratified set, we prove that this algorithm produces a sequence whose accumulation points are stationary, and therefore does not follow the so-called apocalypses described by Levin, Kileel, and Boumal (2022). We illustrate the apocalypse-free property of our method through a numerical experiment on the determinantal variety.

Quadratic hypersurfaces are a natural generalization of affine subspaces, and projections are elementary blocks of algorithms in optimization and machine learning. It is therefore intriguing that no proper studies and tools have been developed to tackle this nonconvex optimization problem. The quadproj package is a user-friendly and documented software that is dedicated to project a point onto a non-cylindrical central quadratic hypersurface.

We consider the problem of minimizing a differentiable function with locally Lipschitz continuous gradient over the real determinantal variety, and present a first-order algorithm designed to find stationary points of that problem. This algorithm applies steps of a retraction-free descent method proposed by Schneider and Uschmajew (2015), while taking the numerical rank into account to attempt rank reductions. We prove that this algorithm produces a sequence of iterates the accumulation points of which are stationary, and therefore does not follow the so-called apocalypses described by Levin, Kileel, and Boumal (2022). Moreover, the rank reduction mechanism of this algorithm requires at most one rank reduction attempt per iteration, in contrast with the one of the $\mathrm{P}^2\mathrm{GDR}$ algorithm introduced by Olikier, Gallivan, and Absil (2022) which can require a number of rank reduction attempts equal to the rank of the iterate in the worst-case scenario., Comment: arXiv admin note: text overlap with arXiv:2201.03962

We address the problem of projecting a point onto a quadratic hypersurface, more specifically a central quadric. We show how this problem reduces to finding a given root of a scalar-valued nonlinear function. We completely characterize one of the optimal solutions of the projection as either the unique root of this nonlinear function on a given interval, or as a point that belongs to a finite set of computable solutions. We then leverage this projection and the recent advancements in splitting methods to compute the projection onto the intersection of a box and a quadratic hypersurface with alternating projections and Douglas-Rachford splitting methods. We test these methods on a practical problem from the power systems literature, and show that they outperform IPOPT and Gurobi in terms of objective, execution time and feasibility of the solution., Comment: 46 pages, 18 figures

We compare two first-order low-rank optimization algorithms, namely $\text{P}^2\text{GD}$ (Schneider and Uschmajew, 2015), which has been proven to be apocalypse-prone (Levin et al., 2021), and its apocalypse-free version $\text{P}^2\text{GDR}$ obtained by equipping $\text{P}^2\text{GD}$ with a suitable rank reduction mechanism (Olikier et al., 2022). Here an apocalypse refers to the situation where the stationarity measure goes to zero along a convergent sequence whereas it is nonzero at the limit. The comparison is conducted on two simple examples of apocalypses, the original one (Levin et al., 2021) and a new one. We also present a potential side effect of the rank reduction mechanism of $\text{P}^2\text{GDR}$ and discuss the choice of the rank reduction parameter.

We study a continuous-time system that solves optimization problems over the set of orthonormal matrices, which is also known as the Stiefel manifold. The resulting optimization flow follows a path that is not always on the manifold but asymptotically lands on the manifold. We introduce a generalized Stiefel manifold to which we extend the canonical metric of the Stiefel manifold. We show that the vector field of the proposed flow can be interpreted as the sum of a Riemannian gradient on a generalized Stiefel manifold and a normal vector. Moreover, we prove that the proposed flow globally converges to the set of critical points, and any local minimum and isolated critical point is asymptotically stable.

Tangent and normal cones play an important role in constrained optimization to describe admissible search directions and, in particular, to formulate optimality conditions. They notably appear in various recent algorithms for both smooth and nonsmooth low-rank optimization where the feasible set is the set $\mathbb{R}_{\leq r}^{m \times n}$ of all $m \times n$ real matrices of rank at most $r$. In this paper, motivated by the convergence analysis of such algorithms, we study, by computing inner and outer limits, the continuity of the correspondence that maps each $X \in \mathbb{R}_{\leq r}^{m \times n}$ to the tangent cone to $\mathbb{R}_{\leq r}^{m \times n}$ at $X$. We also deduce results about the continuity of the corresponding normal cone correspondence. Finally, we show that our results include as a particular case the $a$-regularity of the Whitney stratification of $\mathbb{R}_{\leq r}^{m \times n}$ following from the fact that this set is a real algebraic variety, called the real determinantal variety.

This paper considers the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on the algebraic variety of real matrices of upper-bounded rank. This problem is known to enable the formulation of several machine learning and signal processing tasks such as collaborative filtering and signal recovery. Several definitions of stationarity exist for this nonconvex problem. Among them, Bouligand stationarity is the strongest first-order necessary condition for local optimality. This paper proposes a first-order algorithm that combines the well-known projected-projected gradient descent map with a rank reduction mechanism and generates a sequence in the variety whose accumulation points are Bouligand stationary. This algorithm compares favorably with the three other algorithms known in the literature to enjoy this stationarity property, regarding both the typical computational cost per iteration and empirically observed numerical performance. A framework to design hybrid algorithms enjoying the same property is proposed and illustrated through an example., Comment: The main changes with respect to the first version are listed in section 1.3

The low-rank matrix completion problem can be solved by Riemannian optimization on a fixed-rank manifold. However, a drawback of the known approaches is that the rank parameter has to be fixed a priori. In this paper, we consider the optimization problem on the set of bounded-rank matrices. We propose a Riemannian rank-adaptive method, which consists of fixed-rank optimization, rank increase step and rank reduction step. We explore its performance applied to the low-rank matrix completion problem. Numerical experiments on synthetic and real-world datasets illustrate that the proposed rank-adaptive method compares favorably with state-of-the-art algorithms. In addition, it shows that one can incorporate each aspect of this rank-adaptive framework separately into existing algorithms for the purpose of improving performance., Comment: 22 pages, 12 figures, 1 table

The symplectic Stiefel manifold, denoted by $\mathrm{Sp}(2p,2n)$, is the set of linear symplectic maps between the standard symplectic spaces $\mathbb{R}^{2p}$ and $\mathbb{R}^{2n}$. When $p=n$, it reduces to the well-known set of $2n\times 2n$ symplectic matrices. We study the Riemannian geometry of this manifold viewed as a Riemannian submanifold of the Euclidean space $\mathbb{R}^{2n\times 2p}$. The corresponding normal space and projections onto the tangent and normal spaces are investigated. Moreover, we consider optimization problems on the symplectic Stiefel manifold. We obtain the expression of the Riemannian gradient with respect to the Euclidean metric, which then used in optimization algorithms. Numerical experiments on the nearest symplectic matrix problem and the symplectic eigenvalue problem illustrate the effectiveness of Euclidean-based algorithms.

We address the problem of computing the smallest symplectic eigenvalues and the corresponding eigenvectors of symmetric positive-definite matrices in the sense of Williamson's theorem. It is formulated as minimizing a trace cost function over the symplectic Stiefel manifold. We first investigate various theoretical aspects of this optimization problem such as characterizing the sets of critical points, saddle points, and global minimizers as well as proving that non-global local minimizers do not exist. Based on our recent results on constructing Riemannian structures on the symplectic Stiefel manifold and the associated optimization algorithms, we then propose solving the symplectic eigenvalue problem in the framework of Riemannian optimization. Moreover, a connection of the sought solution with the eigenvalues of a special class of Hamiltonian matrices is discussed. Numerical examples are presented., Comment: 24 pages, 2 figures

We consider a Canonical Polyadic (CP) decomposition approach to low-rank tensor completion (LRTC) by incorporating external pairwise similarity relations through graph Laplacian regularization on the CP factor matrices. The usage of graph regularization entails benefits in the learning accuracy of LRTC, but at the same time, induces coupling graph Laplacian terms that hinder the optimization of the tensor completion model. In order to solve graph-regularized LRTC, we propose efficient alternating minimization algorithms by leveraging the block structure of the underlying CP decomposition-based model. For the subproblems of alternating minimization, a linear conjugate gradient subroutine is specifically adapted to graph-regularized LRTC. Alternatively, we circumvent the complicating coupling effects of graph Laplacian terms by using an alternating directions method of multipliers. Based on the Kurdyka-{\L}ojasiewicz property, we show that the sequence generated by the proposed algorithms globally converges to a critical point of the objective function. Moreover, the complexity and convergence rate are also derived. In addition, numerical experiments including synthetic data and real data show that the graph regularized tensor completion model has improved recovery results compared to those without graph regularization, and that the proposed algorithms achieve gains in time efficiency over existing algorithms.

The symplectic Stiefel manifold, denoted by $\mathrm{Sp}(2p,2n)$, is the set of linear symplectic maps between the standard symplectic spaces $\mathbb{R}^{2p}$ and $\mathbb{R}^{2n}$. When $p=n$, it reduces to the well-known set of $2n\times 2n$ symplectic matrices. Optimization problems on $\mathrm{Sp}(2p,2n)$ find applications in various areas, such as optics, quantum physics, numerical linear algebra and model order reduction of dynamical systems. The purpose of this paper is to propose and analyze gradient-descent methods on $\mathrm{Sp}(2p,2n)$, where the notion of gradient stems from a Riemannian metric. We consider a novel Riemannian metric on $\mathrm{Sp}(2p,2n)$ akin to the canonical metric of the (standard) Stiefel manifold. In order to perform a feasible step along the antigradient, we develop two types of search strategies: one is based on quasi-geodesic curves, and the other one on the symplectic Cayley transform. The resulting optimization algorithms are proved to converge globally to critical points of the objective function. Numerical experiments illustrate the efficiency of the proposed methods., Comment: 28 pages, 7 figures, 5 tables

This paper concerns the minimax center of a collection of linear subspaces. When the subspaces are $k$-dimensional subspaces of $\mathbb{R}^n$, this can be cast as finding the center of a minimum enclosing ball on a Grassmann manifold, Gr$(k,n)$. For subspaces of different dimension, the setting becomes a disjoint union of Grassmannians rather than a single manifold, and the problem is no longer well-defined. However, natural geometric maps exist between these manifolds with a well-defined notion of distance for the images of the subspaces under the mappings. Solving the initial problem in this context leads to a candidate minimax center on each of the constituent manifolds, but does not inherently provide intuition about which candidate is the best representation of the data. Additionally, the solutions of different rank are generally not nested so a deflationary approach will not suffice, and the problem must be solved independently on each manifold. We propose and solve an optimization problem parametrized by the rank of the minimax center. The solution is computed using a subgradient algorithm on the dual. By scaling the objective and penalizing the information lost by the rank-$k$ minimax center, we jointly recover an optimal dimension, $k^*$, and a central subspace, $U^* \in$ Gr$(k^*,n)$ at the center of the minimum enclosing ball, that best represents the data., Comment: 26 pages

We show that H\"older continuity of the gradient is not only a sufficient condition, but also a necessary condition for the existence of a global upper bound on the error of the first-order Taylor approximation. We also relate this global upper bound to the H\"older constant of the gradient. This relation is expressed as an interval, depending on the H\"older constant, in which the error of the first-order Taylor approximation is guaranteed to be. We show that, for the Lipschitz continuous case, the interval cannot be reduced. An application to the norms of quadratic forms is proposed, which allows us to derive a novel characterization of Euclidean norms.

The third edition of the "international - Traveling Workshop on Interactions between Sparse models and Technology" (iTWIST) took place in Aalborg, the 4th largest city in Denmark situated beautifully in the northern part of the country, from the 24th to 26th of August 2016. The workshop venue was at the Aalborg University campus. One implicit objective of this biennial workshop is to foster collaboration between international scientific teams by disseminating ideas through both specific oral/poster presentations and free discussions. For this third edition, iTWIST'16 gathered about 50 international participants and features 8 invited talks, 12 oral presentations, and 12 posters on the following themes, all related to the theory, application and generalization of the "sparsity paradigm": Sparsity-driven data sensing and processing (e.g., optics, computer vision, genomics, biomedical, digital communication, channel estimation, astronomy); Application of sparse models in non-convex/non-linear inverse problems (e.g., phase retrieval, blind deconvolution, self calibration); Approximate probabilistic inference for sparse problems; Sparse machine learning and inference; "Blind" inverse problems and dictionary learning; Optimization for sparse modelling; Information theory, geometry and randomness; Sparsity? What's next? (Discrete-valued signals; Union of low-dimensional spaces, Cosparsity, mixed/group norm, model-based, low-complexity models, ...); Matrix/manifold sensing/processing (graph, low-rank approximation, ...); Complexity/accuracy tradeoffs in numerical methods/optimization; Electronic/optical compressive sensors (hardware)., Comment: 69 pages, 22 extended abstracts, iTWIST'16 website: http://www.itwist16.es.aau.dk

We consider the minimization of a cost function $f$ on a manifold $M$ using Riemannian gradient descent and Riemannian trust regions (RTR). We focus on satisfying necessary optimality conditions within a tolerance $\varepsilon$. Specifically, we show that, under Lipschitz-type assumptions on the pullbacks of $f$ to the tangent spaces of $M$, both of these algorithms produce points with Riemannian gradient smaller than $\varepsilon$ in $O(1/\varepsilon^2)$ iterations. Furthermore, RTR returns a point where also the Riemannian Hessian's least eigenvalue is larger than $-\varepsilon$ in $O(1/\varepsilon^3)$ iterations. There are no assumptions on initialization. The rates match their (sharp) unconstrained counterparts as a function of the accuracy $\varepsilon$ (up to constants) and hence are sharp in that sense. These are the first deterministic results for global rates of convergence to approximate first- and second-order Karush-Kuhn-Tucker points on manifolds. They apply in particular for optimization constrained to compact submanifolds of $\mathbb{R}^n$, under simpler assumptions., Comment: 33 pages, IMA Journal of Numerical Analysis, 2018

The implicit objective of the biennial "international - Traveling Workshop on Interactions between Sparse models and Technology" (iTWIST) is to foster collaboration between international scientific teams by disseminating ideas through both specific oral/poster presentations and free discussions. For its second edition, the iTWIST workshop took place in the medieval and picturesque town of Namur in Belgium, from Wednesday August 27th till Friday August 29th, 2014. The workshop was conveniently located in "The Arsenal" building within walking distance of both hotels and town center. iTWIST'14 has gathered about 70 international participants and has featured 9 invited talks, 10 oral presentations, and 14 posters on the following themes, all related to the theory, application and generalization of the "sparsity paradigm": Sparsity-driven data sensing and processing; Union of low dimensional subspaces; Beyond linear and convex inverse problem; Matrix/manifold/graph sensing/processing; Blind inverse problems and dictionary learning; Sparsity and computational neuroscience; Information theory, geometry and randomness; Complexity/accuracy tradeoffs in numerical methods; Sparsity? What's next?; Sparse machine learning and inference., Comment: 69 pages, 24 extended abstracts, iTWIST'14 website: http://sites.google.com/site/itwist14

Optimal distribution of power among generating units to meet a specific demand subject to system constraints is an ongoing research topic in the power system community. The problem, even in a static setting, turns out to be hard to solve with conventional optimization methods owing to the consideration of valve-point effects which make the cost function nonsmooth and nonconvex. This difficulty gave rise to the proliferation of population-based global heuristics in order to address the multi-extremal and nonsmooth problem. In this paper, we address the economic load dispatch problem (ELDP) with valve-point effect in its classic formulation where the cost function for each generator is expressed as the sum of a quadratic term and a rectified sine term. We propose two methods that resort to piecewise-quadratic surrogate cost functions, yielding surrogate problems that can be handled by mixed-integer quadratic programming (MIQP) solvers. The first method shows that the global solution of the ELDP can often be found by using a fixed and very limited number of quadratic pieces in the surrogate cost function. The second method adaptively builds piecewise-quadratic surrogate under-estimations of the ELDP cost function, yielding a sequence of surrogate MIQP problems. It is shown that any limit point of the sequence of MIQP solutions is a global solution of the ELDP. Moreover, numerical experiments indicate that the proposed methods outclass the state-of-the-art algorithms in terms of minimization value and computation time on practical instances., Comment: 8 pages

Optimization on manifolds is a rapidly developing branch of nonlinear optimization. Its focus is on problems where the smooth geometry of the search space can be leveraged to design efficient numerical algorithms. In particular, optimization on manifolds is well-suited to deal with rank and orthogonality constraints. Such structured constraints appear pervasively in machine learning applications, including low-rank matrix completion, sensor network localization, camera network registration, independent component analysis, metric learning, dimensionality reduction and so on. The Manopt toolbox, available at www.manopt.org, is a user-friendly, documented piece of software dedicated to simplify experimenting with state of the art Riemannian optimization algorithms. We aim particularly at reaching practitioners outside our field.

We consider two Riemannian geometries for the manifold $\mathcal{M}(p,m\times n)$ of all $m\times n$ matrices of rank $p$. The geometries are induced on $\mathcal{M}(p,m\times n)$ by viewing it as the base manifold of the submersion $\pi:(M,N)\mapsto MN^T$, selecting an adequate Riemannian metric on the total space, and turning $\pi$ into a Riemannian submersion. The theory of Riemannian submersions, an important tool in Riemannian geometry, makes it possible to obtain expressions for fundamental geometric objects on $\mathcal{M}(p,m\times n)$ and to formulate the Riemannian Newton methods on $\mathcal{M}(p,m\times n)$ induced by these two geometries. The Riemannian Newton methods admit a stronger and more streamlined convergence analysis than the Euclidean counterpart, and the computational overhead due to the Riemannian geometric machinery is shown to be mild. Potential applications include low-rank matrix completion and other low-rank matrix approximation problems.

Approximate matrix factorization techniques with both nonnegativity and orthogonality constraints, referred to as orthogonal nonnegative matrix factorization (ONMF), have been recently introduced and shown to work remarkably well for clustering tasks such as document classification. In this paper, we introduce two new methods to solve ONMF. First, we show athematical equivalence between ONMF and a weighted variant of spherical k-means, from which we derive our first method, a simple EM-like algorithm. This also allows us to determine when ONMF should be preferred to k-means and spherical k-means. Our second method is based on an augmented Lagrangian approach. Standard ONMF algorithms typically enforce nonnegativity for their iterates while trying to achieve orthogonality at the limit (e.g., using a proper penalization term or a suitably chosen search direction). Our method works the opposite way: orthogonality is strictly imposed at each step while nonnegativity is asymptotically obtained, using a quadratic penalty. Finally, we show that the two proposed approaches compare favorably with standard ONMF algorithms on synthetic, text and image data sets., Comment: 17 pages, 8 figures. New numerical experiments (document and synthetic data sets)

We consider the problem of approximating a multiple-input multiple-output (MIMO) $p\times m$ rational transfer function $H(s)$ of high degree by another $p\times m$ rational transfer function $\hat H(s)$ of much smaller degree, so that the ${\cal H}_2$ norm of the approximation error is minimized. We characterize the stationary points of the ${\cal H}_2$ norm of the approximation error by tangential interpolation conditions and also extend these results to the discrete-time case. We analyze whether it is reasonable to assume that lower-order models can always be approximated arbitrarily closely by imposing only first-order interpolation conditions. Finally, we analyze the ${\cal H}_2$ norm of the approximation error for a simple case in order to illustrate the complexity of the minimization problem., Comment: Paper ID sheet: http://www.inma.ucl.ac.be/~absil/Publi/H2-modred-mimo.htm

We propose an algorithm for solving nonlinear convex programs defined in terms of a symmetric positive semidefinite matrix variable $X$. This algorithm rests on the factorization $X=Y Y^T$, where the number of columns of Y fixes the rank of $X$. It is thus very effective for solving programs that have a low rank solution. The factorization $X=Y Y^T$ evokes a reformulation of the original problem as an optimization on a particular quotient manifold. The present paper discusses the geometry of that manifold and derives a second order optimization method. It furthermore provides some conditions on the rank of the factorization to ensure equivalence with the original problem. The efficiency of the proposed algorithm is illustrated on two applications: the maximal cut of a graph and the sparse principal component analysis problem., Comment: submitted