A Brief Introduction to Manifold Optimization

Manifold optimization is ubiquitous in computational and applied mathematics, statistics, engineering, machine learning, physics, chemistry and etc. One of the main challenges usually is the non-convexity of the manifold constraints. By utilizing the geometry of manifold, a large class of constrained optimization problems can be viewed as unconstrained optimization problems on manifold. From this perspective, intrinsic structures, optimality conditions and numerical algorithms for manifold optimization are investigated. Some recent progress on the theoretical results of manifold optimization are also presented.

conditions as well as state-of-the-art algorithms for manifold optimization in section 3. For some selected practical applications in section 2, a few theoretical results based on manifold optimization are introduced in section 4.
2. Applications of manifold optimization. In this section, we introduce applications of manifold optimization in p-harmonic flow, max-cut problems, phase retrieval, eigenvalue problem, electronic structure calculations, Bose-Einstein condensates, cryo-electron microscopy (Cryo-EM), combinatorial optimization, deep learning and etc.
2.1. P -harmonic flow. P -harmonic flow is used in the color image recovery and medical image analysis. For instance, in medical image analysis, the human brain is often mapped to a unit sphere via a conformal mapping, see Figure 1. By Fig. 1. conformal mapping between the human brain and the unit sphere [57]. establishing a conformal mapping between an irregular surface and the unit sphere, we can handle the complicated surface with the simple parameterizations of the unit sphere. Here, we focus on the conformal mapping between genus-0 surfaces. From [77], a diffeomorphic map between two genus-0 surfaces N 1 and N 2 is conformal if and only if it is a local minimizer of the corresponding harmonic energy. Hence, one effective way to compute the conformal mapping between two genus-0 surfaces is to minimize the harmonic energy of the map. Before introducing the harmonic energy minimization model and the diffeomorphic mapping, we review some related concepts on manifold. Let φ N1 (x 1 , x 2 ) : R 2 → N 1 ⊂ R 3 , φ N2 (x 1 , x 2 ) : R 2 → N 2 ⊂ R 3 be the local coordinates on N 1 and N 2 , respectively. The first fundamental form on N 1 ∂x j , i, j = 1, 2. Given a smooth map f : N 1 → N 2 , whose local coordinate representation is f (x 1 , x 2 ) = (f 1 (x 1 , x 2 ), f 2 (x 1 , x 2 )), the density of the harmonic energy of f is where (g ij ) is the inverse of (g ij ) and the inner product between f * ∂ x i and f * ∂ x j is defined as: where E(f ) is called the harmonic energy of f . Stationary points of E are the harmonic maps from N 1 to N 2 . In particular, If N 2 = R 2 , the conformal map f = (f 1 , f 2 ) is two harmonic functions defined on N 1 . If we consider a p-harmonic map from n dimensional manifold M to n dimensional sphere S n ⊂ R n+1 , the p-harmonic energy minimization problem can be written as F (x) ∈ S n , ∀x ∈ M.
2.2. Max cut. Given a graph G = (V, E) with a set of n vertexes V (|V | = n) and a set of edges E. Denote by the weight matrix W = (w ij ). The max-cut problem is to split V into two nonempty sets (S, V \S) such that the total weights of edges in the cut is maximized. For each vertex i = 1, . . . , n, we define x i = 1 if i ∈ S and −1 otherwise. The maxcut problem can be written as It is NP-hard. By relaxing the rank-1 constraint xx to a positive semidefinite matrix X and further neglecting the rank-1 constraint on X, we obtain the following semidefinite program (SDP) (2.2) max X 0 tr(CX), s.t. X ii = 1, i = 1, · · · , n, where C is the graph Laplacian matrix divided by 4, i.e., C = − 1 4 (diag(W e) − W ). If we decompose X = V V with V := [V 1 , . . . , V n ] ∈ R p×n , a nonconvex relaxation of (2.1) is Vn] tr(CV V ), s.t. V i 2 = 1, i = 1, . . . , n.
It is an optimization problem over multiple spheres.
2.3. Low-rank nearest correlation estimation. Given a symmetric matrix C ∈ S n and a non-negative symmetric weight matrix H ∈ S n , this problem is to find a correlation matrix X of low rank such that the distance weighted by H between X and C is minimized: (2.4) min X 0 1 2 H (X − C) 2 F , s.t. X ii = 1, i = 1, . . . , n, rank (X) ≤ p.
Algorithms for solving (2.4) can be found in [80,34]. Similar to the maxcut problem, we decompose the low-rank matrix X with X = V V , in which V = [V 1 , . . . , V n ] ∈ R p×n . Therefore, problem (2.4) is converted to a quartic polynomial optimization problem over multiple spheres: . . , n.
2.4. Phase retrieval. Given some modules of a complex signal x ∈ C n under linear measurements, a classic model for phase retrieval is to solve (2.5) find where A ∈ C m×n and b ∈ R m . This problem plays an important role in X-ray, crystallography imaging, diffraction imaging and microscopy. Problem (2.5) is equivalent to the following problem, which minimizes the phase variable y and signal variable x simultaneously: In [89], the problem above is rewritten as where M = diag{b}(I − AA † )diag{b} is positive definite. It can be regarded as a generalization of the maxcut problem to complex spheres. If we denote X = uu * , (2.7) can also be modelled as the following SDP problem [21] min tr(M X) s.t. X 0, rank (X) = 1, which can be further relaxed as min tr(M X) s.t. rank (X) = 1, whose constraint is a manifold.
2.5. Bose-Einstein condensates. In Bose-Einstein condensates (BEC), the total energy functional is defined as where w ∈ R d is the spatial coordinate vector,ψ is the complex conjugate of ψ, L z = −i(x∂ − y∂x), V (w) is an external trapping potential, and β, Ω are given constants. The ground state of BEC is defined as the minimizer of the following optimization problem min where the spherical constraint S is The Euler-Lagrange equation of this problem is to find (µ ∈ R, φ(w)) such that Utilizing some proper discretization, such as finite difference, sine pseudospectral and Fourier pseudospectral methods, we obtain a discretized BEC problem where M ∈ N, β are given constants and A ∈ C M ×M is Hermitian. Consider the case that x and A are real. Since x x = 1, multiplying the quadratic term of the objective function by x x, we obtain the following equivalent problem x 2 = 1.
The problem above can be also regarded as the best rank-1 tensor approximation of a fourth-order tensor F [39], with For the complex case, we can obtain a best rank-1 complex tensor approximation problem by a similar fashion. Therefore, BEC is an polynomial optimization problem over single sphere.
2.6. Cryo-EM. The Cryo-EM problem is to reconstruct a three-dimensional object from a series of two-dimensional projected images {P i } of the object. A classic model formulates it into an optimization problem over multiple orthogonality constraints [81] to compute the N corresponding directions {R i } of {P i }, see Figure 2. EachR i ∈ R 3×3 is a three-dimensional rotation, i.e.,R iR i = I 3 and det(R i ) = 1. Letc ij = (x ij , y ij , 0) be the common line of P i and P j (viewed in P i ). If the data are exact, it follows from the Fourier projection-slice theorem [81], the common lines coincide, i.e.,R icij =R jcji .
Since the third column ofR 3 i can be represented by the first two columnsR 1 i and R 2 i asR 3 i = ±R 1 i ×R 2 i , the rotations {R i } can be compressed as a 3-by-2 matrix. Therefore, the corresponding optimization problem is where ρ is a function to measure the distance between two vectors, R i are the first two columns ofR i and c ij are the first two entries ofc ij . In [81], the distance function is set as ρ(u, v) = u − v 2 2 . An eigenvector relaxation and SDP relaxation are also presented in [81].
2.7. Linear eigenvalue problem. Linear eigenvalue decomposition and singular value decomposition are the special cases of optimization with orthogonality constraints. Linear eigenvalue problem can be written as where A ∈ S n is given. Applications from low rank matrix optimization, data mining, principal component analysis and high dimensionality reduction techniques often need to deal with large-scale dense matrices or matrices with some special structures. Although modern computers are developing rapidly, most of the current eigenvalue and singular value decomposition softwares are limited by the traditional design and implementation. In particular, the efficiency may not be significantly improved when working with thousands of CPU cores. From the perspective of optimization, a series of fast algorithms for solving (2.9) was proposed in [65,64,93,96], whose essential parts can be divided into two steps, updating a subspace to approximate the eigenvector space better and extracting eigenvectors by the Rayleigh-Ritz (RR) process.
The main numerical algebraic technique for updating subspaces is usually based on the Krylov subspace, which constructs a series of orthogonal bases sequentially. In [93], the authors propose an equivalent unconstrained penalty function model where µ is a parameter. By choosing an appropriate finite large µ, the authors established its equivalence with (2.9). When µ is chosen properly, the number of saddle points of this model is less than that of (2.9). More importantly, the model allows one to design an algorithm that uses only matrix-matrix multiplication. A Gauss-Newton algorithm for calculating low rank decomposition is developed in [65]. When the matrix to be decomposed is of low rank, this algorithm can be more effective while its complexity is similar to the gradient method but with Q linear convergence. Because the bottleneck of many current iterative algorithms is the RR procedure of the eigenvalue decomposition of smaller dense matrices, the authors of [96] proposed a unified augmented subspace algorithmic framework. Each step iteratively solves a linear eigenvalue problem: where S := span{X, AX, A 2 X, . . . , A k X} with a small k (which can be far less than p). By combining with the polynomial acceleration technique and deflation in classical eigenvalue calculations, it needs only one RR procedure theoretically to reach a high accuracy. When the problem dimension reaches the magnitude of O(10 42 ), the scale of data storage far exceeds the extent that traditional algorithms can handle. In [108], the authors consider to use a low-rank tensor format to express data matrices and eigenvectors. Let N = n 1 n 2 . . . n d with positive integer n 1 , . . . , n d . A vector u ∈ R N can be reshaped as a tensor u ∈ R n1×n2×···×n d , whose entries u i1i2...i d are aligned in reverse lexicographical order, 1 ≤ i µ ≤ n µ , µ = 1, 2, . . . , d. A tensor u can be written as the TT format if its entries can be represented by where U µ (i µ ) ∈ R rµ−1×rµ , i µ = 1, 2, . . . , n µ and fixed dimensions r µ , µ = 0, 1, . . . , d with r 0 = r d = 1. In fact, the components r µ , µ = 1, . . . , d − 1 are often equal to a value r (r is then called the TT-rank). Hence, a vector u of dimension O(n d ) can be stored with O(dnr 2 ) entries if the corresponding tensor u has a TT format. A graphical representation of u can be seen in Figure 3. The eigenvalue problem can be solved based on the subspace algorithm. By utilizing the alternating direction method with suitable truncations, the performance of the algorithm can be further improved. Fig. 3. Graphical representation of a TT tensor of order d with cores Uµ, µ = 1, 2, . . . , d. The first row is u, the second row are its entries u i 1 i 2 ...i d .
The online singular value/eigenvalue decomposition appears in principal component analysis (PCA). The traditional PCA first reads the data and then performs eigenvalue decompositions on the sample covariance matrices. If the data is updated, the principal component vectors need be investigated again based on the new data. Unlike traditional PCA, the online PCA reads the samples one by one and updates the principal component vector in an iterative way, which is essentially a random iterative algorithm of the maximal trace optimization problem. As the sample grows, the online PCA algorithm leads to more accurate main components. An online PCA is proposed and analyzed in [70]. It is proved that the convergence rate is O(1/n) with high probability. A linear convergent VR-PCA algorithm is investigated in [79]. In [59], the scheme in [70] is further proved that under the assumption of Subgaussian's stochastic model, the convergence speed of the algorithm can reach the minimal bound of the information, and the convergence speed is near-global.
2.8. Nonlinear eigenvalue problem. The nonlinear eigenvalue problems from electronic structure calculations are another important source of problems with orthogonality constraints, such as Kohn-Sham (KS) and Hartree-Fock (HF) energy minimization problems. By properly discretizing, KS energy functional can be expressed as where X ∈ C n×p satisfies X * X = I p , ρ = diag(XX * ) is the charge density and µ xc (ρ) := ∂ xc(ρ) ∂ρ and e are vectors in R n with elements all of ones. More specifically, L is a finite dimensional representation of the Laplacian operator, V ion is a constant example, w l represents a discrete reference projection function, ζ l is a constant of ±1, xc is used to characterize exchange-correlation energy. With the KS energy functional, the KS energy minimization problem is defined as Compared to the KS density functional theory, the HF theory can provide a more accurate model. Specifically, it introduces a Fock exchange operator, which is a fourth-order tensor by some discretization, V(·) : C n×n → C n×n . The corresponding Fock energy can be expressed as The HF energy minimization problem is then The first-order optimality conditions of KS and HF energy minimization problems correspond to two different nonlinear eigenvalue problems. Taking KS energy minimization as an example, the first-order optimality condition is where H ks (X) := 1 2 L + V ion + l ζ l w l w * l + diag(( L † )ρ) + diag(µ xc (ρ) * e) and Λ are diagonal matrices. The equation (2.11) is also called the KS equation. The nonlinear eigenvalue problem aims to find some orthogonal eigenvectors satisfying (2.11), while the optimization problem with orthogonality constraints minimizes the objective function under the same constraints. These two problems are connected by the optimality condition and both describe the steady state of the physical system.
The most widely used algorithm for solving the KS equation is the so-called self-consistent field iteration (SCF), which is to solve the following linear eigenvalue problems repeatedly where ρ k = diag(X k X * k ). In practice, to accelerate the convergence, we often replace the charge density ρ k by a linear combination of the previously existing m charge densities In the above expression, α = (α 0 , α 1 , . . . , α m−1 ) is the solution to the following minimization problem: min α e=1 where R = (∆ρ k , ∆ρ k−1 , . . . , ∆ρ k−m+1 ), δ j = ρ j − ρ j−1 and e is a n-dimensional vector of all entries ones. After obtaining ρ mix , we replace H ks (ρ k ) in (2.12) with H ks (ρ mix ) and execute the iteration (2.12). This technique is called charge mixing. For more details, one can refer to [72,73,84]. Since SCF may not converge, many researchers have recently developed optimization algorithms for the electronic structure calculation that can guarantee convergence. In [111], the manifold gradient method is directly extended to solve the KS minimum problem. The algorithm complexity is mainly from the calculation of the total energy and its gradient calculation, and the projection on the Stiefel manifold. Its complexity at each step is much lower than the linear eigenvalue problem and it is easy to be parallelized. Extensive numerical experiments based on the software packages Octopus and RealSPACES show that the algorithm is often more efficient than SCF. In fact, the iteration (2.12) of SCF can be understood as an approximate Newton algorithm in the sense that the complicated part of the Hessian of the total energy is not considered: Since q(X) is only a local approximation model of E ks (X), there is no guarantee that the above model ensures a sufficient decrease of E ks (X). An explicit expression of the complicated part of the Hessian matrix is derived in [92]. Although this part is not suitable for an explicit storage, its operation with a vector is simple and feasible. Hence, the full Hessian matrix can be used to improve the reliability of Newton's method. By adding regularization terms, the global convergence is also guaranteed. A few other related works include [27,113,110,33,55].
The ensemble-based density functional theory is especially important when the spectrum of the Hamiltonian matrix has no significant gaps. The KS energy minimization model is modified by allowing the charge density to contain more wave functions. Specifically, ρ(r) = p i=1 f i |ψ i (r)| 2 where p ≥ p e and the fraction occupation 0 ≤ f i ≤ 1 is to ensure that the total charge density of the total orbit is p, i.e., To calculate the fractional occupancy, the energy functional in the ensemble model introduces a temperature T associated with an entropy αR(f ), where This method is often referred as the KS energy minimization model with temperature or the ensemble KS energy minimization model (EDFT). Similar to the KS energy minimization model, by using the appropriate discretization, the wavefunction can be represented with X = [x 1 , . . . , x p ] ∈ C n×p . The discretized charge density in EDFT can be written as Obviously, ρ(X, f ) is real. The corresponding discretized energy functional is The discretized EDFT model is Although SCF can be generalized to this model, its convergence is still not guaranteed. An equivalent simple model with only one-ball constraint is proposed in [86]. It is solved by a proximal gradient method such that the terms other than the entropy function term are linearized. An explicit solution of the subproblem is then derived and the convergence of the algorithm is established.
2.9. Approximation models for integer programming. Many optimization problems arising from data analysis are NP-hard integer programmings. Spherical constraints and orthogonal constraints are often used to obtain approximate solutions with high quality. Consider optimization problem over the permutation matrices: where f (X) : R n×n → R n is differentiable, and Π n is a collection of n-order permutation matrices This constraint is equivalent to Π n := {X ∈ R n×n : X X = I n , X ≥ 0}.
It is proved in [49] that it is equivalent to an L p -regularized optimization problem over the doubly stochastic matrices, which is much simpler than the original problem. An estimation of the lower bound of the non-zero elements at the stationary points are presented. Combining with the cut plane method, a novel gradient-type algorithm with negative proximal terms is also proposed.
Given k communities S 1 , S 2 , . . . , S k and the set of partition matrix P k n , where the partition matrix X ∈ P k n means X ij = 1, i, j ∈ S t , t ∈ {1, . . . , k}, otherwise X ij = 0. Let A be the adjacency matrix of the network, d i = j A ij , i ∈ {1, . . . , n} and λ = 1/ d 2 . Define the matrix C := −(A − λdd ). The community detection problem in social networks is to find a partition matrix to maximize the modularity function under the stochastic block model: A SDP relaxation of (2.14) is A sparse and low-rank completely positive relaxation technique is further investigated in [106] to transform the model into an optimization problem over multiple nonnegative spheres: where 1 ≤ p ≤ r is usually taken as a small number so that U can be stored for large-scale data sets. The equivalence to the original problem is proved theoretically and an efficient row-by-row type block coordinate descent method is proposed. In order to quickly solve network problems whose dimension is more than 10 million, an asynchronous parallel algorithm is further developed.

Deep learning.
Batch normalization is a very popular technique in deep neural networks. It avoids internal covariance translation by normalizing the input of each neuron. The space formed by its corresponding coefficient matrix can be regarded as a Riemannian manifold. For a deep neural network, batch normalization usually involves input processing before the nonlinear activation function. Define x and w as the outputs of the previous layer and the parameter vector for the current neuron, the batch normalization of z := w x can be written as where u := w/|w|, E(z) is the expectation of random variable z and R xx are the covariance matrices of x. From the definition, we have BN(w x) = BN(u x) and Therefore, the use of the batch standardization ensures that the model does not explode with large learning rates and that the gradient is invariant to linear scaling during propagation.
Since BN(cw x) = BN(w x) holds for any constant c , the optimization problem for deep neural networks using batch normalization can be written as where L(X) is the loss function, S n−1 is a sphere in R n (can also be viewed as a Grassmann manifold), n 1 , . . . , n m are the dimensions of the weight vectors, m is the number of weight vectors, and l is the number of remaining parameters to be decided, including deviations and other weight parameters. For more information, we refer to [26].
2.11. Sparse PCA. In the traditional PCA, the obtained principle eigenvectors are usually not sparse, which leads to high computational cost for computing the principle components. Spare PCA [51] wants to find principle eigenvectors with few non-zero elements. The mathematical formulation is where X 1 = ij |X ij | and ρ > 0 is a trade-off parameter. When ρ = 0, this reduces to the traditional PCA problem. For ρ > 0, the term X 1 plays a role to promote sparsity. Problem (2.16) is a non-smooth optimization problem on the Stiefel manifold.
2.12. Low-rank matrix completion. The low-rank matrix completion problem has important applications in computer vision, pattern recognitions, statistics, etc. It can be formulated as where X is the matrix that we want to recover (some of its entries are known) and Ω is the index set of observed entries. Due to the difficulty of the rank, a popular approach is to relax it into a convex model using the nuclear norm. The equivalence between this convex problem and the nonconvex problem (2.17) is ensured under certain conditions. Another way is to use a low rank decomposition on X and then solve the corresponding unconstrained optimization problem [95]. If the rank of the ground-truth matrix A is known, an alternative model for a fixed-rank matrix completion is where P Ω is a projection with P Ω (X) ij = X ij , (i, j) ∈ Ω and 0 otherwise, and r = rank (A). The set Fr(m, n, r) := {X ∈ R m×n : rank (X) = r} is a matrix manifold, called fixed-rank manifold. The related geometry is analyzed in [87]. Consequently, problem (2.18) can be solved by optimization algorithms on manifold. Problem (2.18) can deal with Gaussian noise properly. For data sets with a few outliers, the robust low-rank matrix completion problem (with the prior knowledge r) considers: ) is a non-smooth optimization problem on the fixed-rank matrix manifold. For some related algorithms for (2.18) and (2.19), the readers can refer to [91,23].

Sparse blind deconvolution.
Blind deconvolution is to recover a convolution kernel a 0 ∈ R k and signal x 0 ∈ R m from their convolution where y ∈ R m . Since there are infinitely many pairs (a 0 , x 0 ) satisfying this condition, this problem is often ill-conditioned. To overcome this issue, some regularization terms and extra constraints are necessary. The sphere-constrained sparse blind deconvolution reformulate the problem as where µ is a parameter to control the sparsity of the signal x. This is a non-smooth optimization problem on the product manifold of a sphere and R m . Some related background and the corresponding algorithms can be found in [112].
2.14. Non-negative PCA. Since the principle eigenvectors obtained by the traditional PCA may not be sparse, one can enforce the sparsity by adding nonnegativity constraints. The problem is formulated as where A = [a 1 , . . . , a k ] ∈ R n×k are given data points. Under the constraints, the variable X has at most one non-zero element in each row. This actually helps to guarantee the sparsity of the principle eigenvectors. Problem (2.20) is an optimization problem with manifold and non-negative constraints. Some related information can be found in [102,68].
2.15. K-means clustering. K-means clustering is a fundamental problem in data mining. Given n data points (x 1 , x 2 , . . . , x n ) where each data point is a ddimensional vector, k-means is to partition them into k clusters S := {S 1 , S 2 , . . . , S k } such that the within-cluster sum of squares is minimized. Each data point belongs to the cluster with the nearest mean. The mathematical form is x∈Si x is the center of i-th cluster and card(S i ) is the cardinality of S i . Equivalently, problem (2.21) can be written as [24,61,99]: 2 is the squared Euclidean distance matrix. Problem (2.22) is a minimization over the Stiefel manifold with linear constraints and non-negative constraints.
3. Algorithms for manifold optimization. In this section, we introduce a few state-of-the-art algorithms for optimization problems on Riemannian manifold. Let us start from the concepts of manifold optimization.

Preliminaries on Riemannian manifold.
A d-dimensional manifold M is a Hausdorff and second-countable topological space, which is homeomorphic to the d-dimensional Euclidean space locally via a family of charts. When the transition maps of intersecting charts are smooth, the manifold M is called a smooth manifold. Intuitively, the tangent space T x M at a point x of a manifold M is the set of the tangent vectors of all the curves at x. Mathematically, a tangent vector ξ x to M at x is a mapping such that there exists a curve γ on M with γ(0) = x, satisfying where x (M) is the set of all real-valued functions f defined in a neighborhood of x in M. Then, the tangent space T x M to M is defined as the set of all tangent vectors to M at x. If M is equipped with a smoothly varied inner product g(·, ·) := ·, · x on the tangent space, then (M, g) is a Riemannian manifold. In practice, different Riemannian metrics may be investigated to design efficient algorithms. The is any curve on the manifold that satisfies γ(0) = x andγ(0) = ξ. The Riemannian Hessian Hess f (x) is a mapping from the tangent space T x M to the tangent space T x M: where∇ is the Riemannian connection [2]. For a function f defined on the manifold, if it can be extended to the ambient Euclidean space R n×p , we have its Riemannian gradient grad f and Riemannian Hessian Hess f : where D is the Euclidean derivative. More detailed information on the related backgrounds can be found in [2]. We next briefly introduce some typical manifolds.
At t = 0, we haveẋ(0)x + x ẋ(0) = 0. Hence, the tangent space is The projection operator is defined as For a function defined on Sp(n − 1) with respect to the Euclidean metric , its Riemannian gradient and Hessian at x can be represented by • Stiefel manifold St(n, p) := {X ∈ R n×p : X X = I p }. By a similar calculation as the spherical case, we have its tangent space: The projection operator onto T X St(n, p) is where sym(Z) := (Z + Z )/2. Given a function defined on St(n, p) with respect to the Euclidean metric g X (U, V ) = tr(U V ), U, V ∈ T X St(n, p), its Riemannian gradient and Hessian at X can be represented by The projection operator onto T X Ob(n, p) is Given a function defined on Ob(n, p) with respect to the Euclidean metric, its Riemannian gradient and Hessian at X can be represented by It denotes the set of all p-dimensional subspaces of R n . This manifold is different from other manifolds mentioned above. It is a quotient manifold since each element is an equivalent class of n×p matrices. From the definition of Grass(p, n), the equivalence relation ∼ is defined as Its element is of the form Then Grass(n, p) is a quotient manifold of St(n, p), i.e., St(n, p)/ ∼. Due to this equivalence, a tangent vector ξ of T X Grass(n, p) may have many different representations in its equivalence class. To find the unique representation, a horizontal space [2, Section 3.5.8] is introduced. For a given X ∈ R n×p with X X = I p , the horizontal space is Here, a function of the horizontal space is similar to the tangent space when computing the Riemannian gradient and Hessian. We have the projection onto the horizontal space Given a function defined on Grass(n, p) with respect to the Euclidean metric g X = tr(U V ), U, V ∈ H X Grass(n, p), its Riemannian gradient and Hessian at X can be represented by grad f (X) = P H X Grass(n,p) (∇f (X)), • Fixed-rank manifold Fr(n, p, r) := {X ∈ R n×p : rank (X) = r} is a set of all n × p matrices of rank r. Using the singular value decomposition (SVD), this manifold can be represented equivalently by where U ⊥ and V ⊥ are the orthogonal complements of U and V , respectively. The projection operator onto the tangent space is Given a function defined on Fr(n, p, r) with respect to the Euclidean metric g X (U, V ) = tr(U V ), its Riemannian gradient and Hessian at X = U ΣV can be represented by grad f (X) = P T X Fr(n,p,r) (∇f (X)), The set of symmetric positive definite matrices, i.e., SPD(n) = {X ∈ R n×n : X = X, X 0} is a manifold. Its tangent space at X is We have the projection onto T X SPD(n): Given a function defined on SPD(n, p) with respect to the Euclidean metric g X (U, V ) = tr(U V ), U, V ∈ T X SPD(n), its Riemannian gradient and Hessian at X can be represented by grad f (X) = P T X SPD(n) (∇f (X)), • The set of rank-r symmetric positive semidefinite matrices, i.e., FrPSD(n, r) = {X ∈ R n×n : X = X , X 0, rank (X) = r}. This manifold can be reformulated as which is a quotient manifold. The horizontal space at Y is We have the projection operator onto T Y H FrPSD(n,r) where the skew-symmetric matrix Ω is the unique solution of the Sylvester Given a function f with respect to the Euclidean metric g Y (U, V ) = tr(U V ), U, V ∈ T Y H FrPSD(n,r) , its Riemannian gradient and Hessian can be represented by 3.2. Optimality conditions. We next present the optimality conditions for manifold optimization problem in the following form where E and I denote the index sets of equality constraints and inequality constraints, respectively, and c i : M → R, i ∈ E ∪I are smooth functions on M. We mainly adopt the notions in [100]. By keeping the manifold constraint, the Lagrangian function of Then the first-order necessary conditions can be described as follows.
Suppose that x * is a local minima of (3.3) and that the LICQ holds at x * , then there exist Lagrangian multipliers λ * i , i ∈ E ∪I such that the following KKT conditions hold: Let x * and λ * i , i ∈ E ∪I be one of the solution of the KKT conditions (3.4). Similar to the case without the manifold constraint, we define a critical cone C(x * , λ * ) as Then we have the following second-order necessary and sufficient conditions.
• Second-order necessary conditions: Suppose that x * is a local minima of (3.3) and the LICQ holds at x * . Let λ * be the multipliers such that the KKT conditions (3.4) hold. Then we have where Hess x L(x * , λ * ) is the Riemannian Hessian of L with respect to x at (x * , λ * ). • Second-order sufficient conditions: Suppose that x * and λ * satisfy the KKT conditions (3.4). If we further have then x * is a strict local minima of (3.4).
Suppose that we have only the manifold constraint, i.e., E ∪ I is empty. For a smooth function f on the manifold M, the optimality conditions take a similar form to the Euclidean unconstrained case. Specifically, if x * is a first-order stationary point, then it holds that grad f (x * ) = 0.
If x * is a second-order stationary point, then then x * is a strict local minimum. For more details, we refer the reader to [100].
3.3. First-order type algorithms. From the perspective of Euclidean constrained optimization problems, there are many standard algorithms which can solve this optimization problem on manifold. However, since the intrinsic structure of manifolds is not considered, these algorithms may not be effective in practice. By doing curvilinear search along the geodesic, a globally convergent gradient descent method is proposed in [32]. For Riemannian conjugate gradient (CG) methods [82], the parallel translation is used to construct the conjugate directions. Due to the difficulty of calculating geodesics (exponential maps) and parallel translations, computable retraction and vector transport operators are proposed to approximate the exponential map and the parallel translation [2]. Therefore, more general Riemannian gradient descent methods and CG methods together with convergence analysis are obtained in [2]. These algorithms have been successfully applied to various applications [87,54]. Numerical experiments exhibit the advantage of using geometry of the manifold. A proximal Riemannian gradient method is proposed in [41]. Specifically, the objective function is linearized using the first-order Taylor expansion on manifold and a proximal term is added. The original problem is then transformed into a series of projection problems on the manifold. For general manifolds, the existence and uniqueness of the projection operator can not be guaranteed. But when the given manifold satisfies certain differentiable properties, the projection operator is always locally well-defined and is also a specific retraction operator [3]. Therefore, in this case, the proximal Riemannian gradient method coincides with the Riemannian gradient method. By generalizing the adaptive gradient method in [30], an adaptive gradient method on manifold is also presented in [41]. In particular, optimization over Stiefel manifold is an important special case of Riemannian optimization. Various efficient retraction operators, vector transport operators and Riemannian metric have been investigated to construct more practical gradient descent and CG methods [94,48,114]. Nonretraction based first-order methods are also developed in [33].
We next present a brief introduction of first-order algorithms for manifold optimization. Let us start with the retraction operator R. It is a smooth mapping from the tangent bundle T M := ∪ x∈M T x M to M, and satisfies where R x is the retraction operator R at x. The well-posedness of the retraction operator is shown in Section 4.1.3 of [2]. The retraction operator provides an efficient way to pull the points from the tangent space back onto the manifold. Let ξ k ∈ T x M be a descent direction, i.e., grad f (x k ), ξ k < 0. Another important concept on manifold is the vector transport operator T . It is a smooth mapping from the product of tangent bundles T M T M to the tangent bundle T M, and satisfies the following properties.
• There exists a retraction R associated with T , i.e., The vector transport is a generalization of the parallel translation [2, Section 5.4]. The general feasible algorithm framework on the manifold can be expressed as where t k is a well-chosen step size. Similar to the line search method in Euclidean space, the step size t k can be obtained by the curvilinear search on the manifold. Here, we take the Armijo search as an example. Given ρ, δ ∈ (0, 1), the monotone and nonmonotone search try to find the smallest integer h to such that respectively, where grad f (x k ), ξ k x k := g x k (grad f (x k ), ξ k ), t k = γ k δ h and γ k is an initial step size. The reference value C k+1 is a convex combination of C k and f (x k+1 ) and is calculated via C k+1 = ( Q k C k + f (x k+1 ))/Q k+1 , where C 0 = f (x 0 ), Q k+1 = Q k + 1 and Q 0 = 1. From the Euclidean optimization, we know that the Barzilai-Borwein (BB) step size often accelerates the convergence. The BB step size can be generalized to Riemannian manifold [41] as and T x k−1 →x k : T x k−1 M → T x k M denotes an appropriate vector transport mapping connecting x k−1 and x k ; see [2,47]. When M is a submanifold of an Euclidean space, the Euclidean differences are an alternative choice if the Euclidean inner product is used in (3.8). This choice is often attractive since the vector transport is not needed [94,41]. We note that the differences between first-and second-order algorithms are mainly due to their specific ways of acquiring ξ k . In practice, the computational cost and convergence behavior of different retraction operators differ a lot. Similarly, the vector transport plays an important role in CG methods and quasi-Newton methods (we will introduce them later). There are many studies on the retraction operators and vector transports. Here, we take the Stiefel manifold St(n, p) as an example to introduce several different retraction operators at the current point X for a given step size τ and descent direction −D.
• Exponential map [31] R geo This scheme needs to calculate an exponent of a 2p-by-2p matrix and an QR decomposition of an n-by-p matrix. From [31], an explicit form of parallel translation is unknown. • Cayley transform [92] (3.9) R wy where U = [P X D, X], V = [X, −P X D] ∈ R n×(2p) with P X := (I − 1 2 XX ). When p < n/2, this scheme is much cheaper than the exponential map. The associated vector transport is [114] T wy η X (ξ X ) = I − The computational cost is lower than the Cayley transform, but the Cayley transform gives a better approximation to the exponential map. The associated vector transport is then defined as [42] T pd where Y = R X η X and vec(Ω) = (Y (X +η X ))⊕(Y (X +η X )) −1 vec(Y ξ X − ξ X Y ) and ⊕ is the Kronecker sum, i.e., A ⊕ B = A ⊗ I + I ⊗ B with Kronecker product ⊗. Numerical experiments show that more iterations may be required compared to the Cayley transform.
. It can be seen as an approximation of the polar decomposition. The main cost is the QR decomposition of a n-by-p matrix. The associated vector transport is defined as [2, Example 8.1.5] where Y = R X (η X ). The vector transport above requires an associated retraction. Removing the dependence of the retraction, a new class of vector transports is introduced in [46]. Specifically, a jointly smooth operator L(x, y) : T x M → T y M is defined. In addition, L(x, x) is required to be an identity for all x. For a d-dimensional submanifold M of n-dimensional Euclidean space, two popular vector transports are defined by the projection [2, Section 8.1.3] L pj (x, y)ξ x = P TyM (ξ x ), and by parallelization [46] L pl (x, y)ξ x = B y B † x ξ x , where B : V → R n×d : z → B z is a smooth tangent basis field defined on an open neighborhood V of M and B † z is the pseudo-inverse of B z . With the tangent basis B z , we can also represent the vector transport mentioned above intrinsically, which sometimes reduces computational cost significantly [45].
To better understand Riemannian first-order algorithms, we present a Riemannian gradient method [41] in Algorithm 1. One can easily see that the difference to the Euclidean case is an extra retraction step.
The convergence of Algorithm 1 [40, Theorem 1] is given as follows.
Then, compute C k , Q k and find a step size t k satisfying (3.7).
Theorem 3. Let {x k } be a sequence generated by Algorithm 1 using the nonmonotone line search (3.7). Suppose that f is continuously differentiable on the manifold M. Then, every accumulation point x * of the sequence {x k } is a stationary point of problem (1.1), i.e., it holds grad f (x * ) = 0.
Proof. At first, by using grad f (x k ), η k x k = − grad f (x k ) 2 x k < 0 and applying [103, Lemma 1.1], we have f (x k ) ≤ C k and x k ∈ L for all k ∈ N. Next, due to there always exists a positive step size t k ∈ (0, γ k ] satisfying the monotone and nonmonotone Armijo conditions (3.6) and (3.7), respectively. Now, let x * ∈ M be an arbitrary acccumulation point of {x k } and let {x k } K be a corresponding subsequence that converges to x * . By the definition of C k+1 and (3.6), we have Hence, {C k } is monotonically decreasing and converges to some limitC ∈ R ∪ {−∞}.
Using f (x k ) → f (x * ) for K k → ∞, we can inferC ∈ R and thus, we obtain x k } → 0. Let us now assume grad f (x * ) = 0. In this case, we have {t k } K → 0 and consequently, by the construction of Algorithm 1, the step size δ −1 t k does not satisfy (3.7), i.e., it holds for all k ∈ K sufficiently large. Since the sequence {η k } K is bounded, the rest of the proof is now identical to the proof of [2, Theorem 4.3.1]. In particular, applying the mean value theorem in (3.10) and using the continuity of the Riemannian metric, we can easily derive a contradiction. We refer to [2] for more details.

Second-order type algorithms.
A gradient-type algorithm usually is fast in the early iterations, but it often slows down or even stagnates when the generated iterations are close to an optimal solution. When a high accuracy is required, secondorder type algorithms may have its advantage.
By utilizing the exact Riemannian Hessian and different retraction operators, Riemannian Newton methods, trust-region methods, adaptive regularized Newton method have been proposed in [85,1,2,41]. When the second-order information is not available, the quasi-Newton type method becomes necessary. As in the Riemannian CG method, we need the vector transport operator to compare different tangent vectors from different tangent spaces. In addition, extra restrictions on the vector transport and the retraction are required for better convergence property or even convergence [74,75,78,42,44,46,43]. Non-vector-transport based quasi-Newton method is also explored in [38].

Riemannian trust-region method.
One of the popular second-order algorithms is a Riemannian trust-region (RTR) algorithm [1,2]. At the k-th iteration x k , by utilizing the Taylor expansion on manifold, RTR constructs the following subproblem on the Tangent space: where ∆ k is the trust-region radius. In [69], extensive methods for solving (3.11) are summarized. Among them, the Steihaug CG method, also named as truncated CG method, is most popular due to its good properties and relatively cheap computational cost. By solving this trust-region subproblem, we obtain a direction ξ k ∈ T x k M satisfying the so-called Cauchy decrease. Then a trial point is computed as z k = R x k (ξ k ), where the step size is chosen as 1. To determine the acceptance of z k , we compute the ratio between the actual reduction and the predicted reduction When ρ k is greater than some given parameter 0 < η 1 < 1, z k is accepted. Otherwise, z k is rejected. To avoid the algorithm stagnating at some feasible point and promote the efficiency as well, the trust-region radius is also updated based on ρ k . The full algorithm is presented in Algorithm 2.
For the global convergence, the following assumptions are necessary for secondorder type algorithms on manifold. Use the truncated CG method to obtain ξ k by solving (3.11). 5 Compute the ratio ρ k in (3.12).
Assumption 5. There exists two constants β RL > 0 and δ RL > 0 such that for all x ∈ M and ξ ∈ T x M with ξ = 1, Then the global convergence to a stationary point [2,Theorem 7.4.2] is presented as follows.

Adaptive regularized Newton method.
From the perspective of Euclidean approximation, an adaptive regularized Newton algorithm (ARNT) is proposed for specific and general manifold optimization problems [92,97,41]. In the subproblem, the objective function is constructed by the second-order Taylor expansion in the Euclidean space and an extra regularization term, while the manifold constraint is kept. Specifically, the mathematical formulation is where H k is the Euclidean Hessian or its approximation. From the definition of Riemannian gradient and Hessian, we have 15 Update ξ k according to (3.15).
where U ∈ T x k M, P ⊥ Tx k M := I − P Tx k M is the projection onto the normal space and the Weingarten map x k M is a symmetric linear operator which is related to the second fundamental form of M. To solve (3.13), a modified CG method is proposed in [41] to solve the Riemannian Newton equation at x k , Since Hessm k (x k ) may not be positive definite, CG may be terminated if a direction with negative curvature, says d k , is encountered. Different from the truncated CG method used in RTR, a linear combination of s k (the output of the truncated CG method) and the negative curvature direction d k is used to construct a descent direction A detailed description on the modified CG method is presented in Algorithm 3. Then, Armijo search along ξ k is adopted to obtain a trial point z k . After obtaining z k , we compute the following ratio between the actual reduction and the predicted reduction, If ρ k ≥ η 1 > 0, then the iteration is successful and we set x k+1 = z k ; otherwise, the iteration is not successful and we set x k+1 = x k , i.e., x k , otherwise.
2 while stopping conditions not met do 3 Compute a new trial point z k by doing Armijo search along ξ k obatined by Algorithm 3.
5 Update x k+1 from the trial point z k based on (3.17). 6 Update σ k according to (3.18).
The regularization parameter σ k+1 is updated as follows where 0 < η 1 ≤ η 2 < 1 and 0 < γ 0 < 1 < γ 1 ≤ γ 2 . These parameters determine how aggressively the regularization parameter is adjusted when an iteration is successful or unsuccessful. Putting these features together, we obtain Algorithm 4, which is dubbed as ARNT. We next present the convergence property of Algorithm 4 with the exact Euclidean Hessian (i.e., H k = ∇ 2 f (x k )) starting from a few assumptions.
We note that the assumptions (A. For the local convergence rate, we make the following assumptions.
Assumption 9. Let {x k } be generated by Algorithm 4.
for all x ∈ M, all ξ ∈ T x M with ξ x = 1 and all t < δ R .
Then we have the following results on the local convergence rate.
The detailed convergence analysis can be found in [41].

Quasi-Newton type methods. When the Riemannian Hessian Hess f (x)
is computationally expensive or even not available, quasi-Newton-type methods turn out to be an attractive approach. In literatures [74,75,78,42,44,46,43], extensive variants of quasi-Newton methods are proposed. Here, we take the Riemannian Broyden-Fletcher-Goldfarb-Shanno (BFGS) as an example to show the general idea of quasi-Newton methods on Riemannian manifold. Similar to the quasi-Newton method in the Euclidean space, an approximation B k+1 should satisfy the following secant equation where s k = T S α k ξ k α k ξ k and y k = β −1 k grad f (x k+1 ) − T S α k ξ k grad f (x k ) with parameter β k . Here, α k and ξ k is the step size and the direction used in the k-th iteration. T S α k ξ k is an isometric vector transport operator associated with the retraction R, i.e., Additionally, T S should satisfy the following locking condition, Then, the scheme of the Riemannian BFGS is With this choice of β k and the isometric property of T S , we can guarantee the positive Algorithm 5: Riemannian BFGS method 1 Input: Initial guess x 0 ∈ M, isometric vector transport T S associated with the retraction R, initial Riemannian Hessian approximation B 0 : T x0 M → T x0 M, which is symmetric positive definite, Wolfe condition parameters 0 < c 1 < 1 2 < c 2 < 1.
Obtain x k+1 by doing a Wolfe search along ξ k , i.e., finding α k > 0 such that the following two conidtions are satisfied Update B k+1 by (3.19).
definiteness of B k+1 . After obtaining the new approximation B k+1 , the Riemannian BFGS method solves the following linear system to get ξ k . The detailed algorithm is presented in Algorithm 5. The choice of β k = 1 can also guarantee the convergence but with more strict assumptions. One can refer to [46] for the convergence analysis. Since the computation of differentiated retraction may be costly, authors in [43] investigate another way to preserve the positive definiteness of the BFGS scheme. Meanwhile, the Wolfe search is replaced by the Armijo search. As a result, the differentiated retraction can be avoided and the convergence analysis is presented as well.
The aforementioned quasi-Newton methods rely on the vector transport operator. When the vector transport operation is computationally costly, these methods may be less competitive. Noticing the structure of the Riemannian Hessian Hess f (x k ), i.e., where the second term W x k (U, P ⊥ Tx k M (∇f (x k ))) is often much cheaper than the first term P Tx k M (∇ 2 f (x k )[U ]). Similar to the quasi-Newton methods in unconstrained nonlinear least square problems [52] [83, Chapter 7], we can focus on the construction of an approximation of the Euclidean Hessian ∇ 2 f (x k ) and use exact formulations of remaining parts. Furthermore, if the Euclidean Hessian itself consists of cheap and expensive parts, i.e., (3.20) where the computational cost of H e (x k ) is much more expensive than H c (x k ), an approximation of ∇ 2 f (x k ) is constructed as where C k is an approximation of H e (x k ) obtained by a quasi-Newton method in the ambient Euclidean space. If an objective function f is not equipped with the structure (3.20), H k is a quasi-Newton approximation of ∇ 2 f (x k ). In the construction of the quasi-Newton approximation, a Nyström approximation technique [38, Section 2.3] is explored, which turns to be a better choice than the BB type initialization [69,Chapter 6]. Since the quasi-Newton approximation is constructed in the ambient Euclidean space, the vector transport is not necessary. Then, subproblem (3.13) is constructed with H k . From the expression of the Riemannian Hessian Hessm k in (3.14), we see that subproblem (3.13) gives us a way to approximate the Riemannian Hessian when an approximation H k to the Euclidean Hessian is available. The same procedures of ARNT can be utilized for (3.13) with the approximate Euclidean Hessian H k . An adaptive structured quasi-Newton method given in [38] is presented in Algorithm 6.
2 while stopping conditions not met do 3 Check the structure of ∇ 2 f (x k ) to see if it can be written in a form as (3.20). 4 Construct an approximation H k by utilizing a quasi-Newton method. 5 Construct and solve the subproblem (3.13) (by using the modified CG method or the Riemannian gradient type method) to obtain a new trial point z k . 6 Compute the ratio ρ k via (3.12).
7 Update x k+1 from the trial point z k based on (3.17). 8 Update τ k according to (3.18).
To explain the differences between the two quasi-Newton algorithms more straightforwardly, we take the HF total energy minimization problem (2.10) as an example. From the calculation in [38], we have the Euclidean gradients where H ks (X) := 1 2 L+V ion + l ζ l w l w * l +Diag(( L † )ρ)+Diag(µ xc (ρ) * e) and H hf (X) = H ks (X) + V(XX * ). The Euclidean Hessian of E ks and E f along a matrix U ∈ C n×p are ∂ρ 2 e (X U + X Ū )e X, Since ∇ 2 E f (X) is significantly more expensive than ∇ 2 E ks (X), we only need to ap- Then, a quasi-Newton approximation C k of ∇ 2 E f is obtained without requiring vector transport. By adding the exact formulation of ∇ 2 E ks (X k ), we have an approximation H k , i.e., H k = ∇ 2 E ks + C k .
A Nyström approximation for C k is also investigated. Note that the spectrum of ∇ 2 E ks (X) dominates the spectrum of ∇ 2 E f (X). The structured approximation H k is more reliable than a direct quasi-Newton approximate to ∇ 2 E hf (X) because the spectrum of ∇ 2 E ks is inherited from the exact form. The Remaining procedure is to solve subproblem (3.13) to update X k .
3.5. Stochastic algorithms. For problems arising from machine learning, the objective function f is often a summation of a finite number of functions For unconstrained situations, there are many efficient algorithms, such as Adam, Adagrad, RMSProp, Adelta, SVRG, etc. One can refer to [58]. For the case with manifold constraints, combining with retraction operators and vector transport operator, these algorithms can be well generalized. However, in the implementation, due to the considerations of the computational costs of different parts, they may have different versions. Riemannian stochastic gradient method is first developed in [14]. Later, a class of first-order methods and their accelerations are investigated for geodesically convex optimization in [105,66]. With the help of parallel translation or vector transport, Riemannian SVRG methods are generalized in [104,76]. In consideration of the computational cost of the vector transport, non-vector transport based Riemannian SVRG is proposed in [50]. Since an intrinsic coordinate system is absent, the coordinate-wise update on manifold should be further investigated. A compromised approach for Riemannian adaptive optimization methods on product manifolds are presented in [10].
Here, the SVRG algorithm [50] is taken as an example. At the current point X s,k , we first calculate the full gradient G(X s,k ), then randomly sample a subscript from 1 to m and use this to construct a stochastic gradient with reduced variance as G(X s,k , ξ s,k ) = ∇f (X s,0 ) + ∇f i s,k (X s,k ) − ∇f i s,k (X s,0 ) , finally move along this direction with a given step size to next iteration point For Riemannian SVRG [50], after obtaining the stochastic gradient with reduced Euclidean variance, it first projects this gradient to the tangent space ξ s,k = P T X s,k M (G(X s,k )).
Then, the following retraction step X s,k+1 = R X s,k (−τ s ξ s,k ), Calculate a random Euclidean gradient G(X s,k ) Calculate a random Riemann gradient ξ s,k = P T X s,k M (G(X s,k )).
Update X s,k+1 in the following format is executed to get the next feasible point. The detailed version is outlined in Algorithm 7.
3.6. Algorithms for Riemannian non-smooth optimization. As shown in subsections 2.11 to 2.15, many practical problems are with non-smooth objective function and manifold constraints, i.e., where g is smooth and h is non-smooth. Riemannian subgradient methods [29,15] are firstly investigated to solve this kind of problems and their convergence analysis is exhibited in [36] with the help of Kurdyka-Lojasiewicz (K L) inequalities. For locally Lipschitz functions on Riemannian manifolds, a gradient sampling method and a nonsmooth Riemannian trust-region method are proposed in [35,37]. Proximal gradient methods on manifold are presented in [5,28], where the inner subproblem is solved inexactly by subgradient type methods. The corresponding complexity analysis is given in [11,12]. Different from the constructions of the subproblem in [5,28], a more tractable subproblem without manifold constraints is investigated in [25] for convex h(x). By utilizing the semi-smooth Newton method [98], the proposed proximal gradient method on manifold enjoys a faster convergence. Another class of methods is based on operator-splitting techniques. Some variants of the alternating direction method of multipliers (ADMM) are studied in [56,53,90,107,13,60].
We briefly introduce the proximal gradient method on manifold [25] here. Assume that the convex function h is Lipschitz continuous. At each iteration x k , the following subproblem is constructed where t > 0 is a step size. Given a retraction R, problem (3.22) can be seen as a firstorder approximation of f (R x k (d)) near the zero element 0 x k on T x k M. Specifically, it follows from the definition of the Riemannian gradient that grad g(x k ) = ∇g(R x k (0)). From the Lipschitz continuous property of h and the definition of R, we have where L h is the Lipschitz constant of h. Therefore, we conclude Then the next step is to solve (3.22). Since (3.22) is convex and with linear constraints, the KKT conditions are sufficient and necessary for the global optimality. Specifically, we have , it is proved in [25] that E is monotone and then the semi-smooth Newton method in [98] is utilized to solve the nonlinear equation E(λ) = 0 to obtain a direction d k . Combining with a curvilinear search along d k with R x k , the decrease on f is guaranteed and the global convergence is established.
3.7. Complexity Analysis. The complexity analysis of the Riemannian gradient method and the Riemannian trust region method has been studied in [16]. Similar to the Euclidean unconstrained optimization, the Riemannian gradient method (using a fixed step size or Armijo curvilinear search) converges to grad f (x) ≤ ε up to O(1/ε 2 ) steps. Under mild assumptions, a modified Riemannian trust region method converges to grad f (x) ≤ ε, Hess f (x) − √ εI at most O(max{1/ε 1.5 , 1/ε 2.5 }) iterations. For objective functions with multi-block convex but non-smooth terms, an ADMM of complexity of O(1/ε 4 ) is proposed in [107]. For the cubic regularization methods on the Riemannian manifold, recent studies [109,4] show a convergence to grad f (x) ≤ ε, Hess f (x) − √ εI with complexity of O(1/ε 1.5 ).

Analysis for manifold optimization.
4.1. Geodesic convexity. For a convex function in the Euclidean space, any local minima is also a global minima. An interesting extension is the geodesic convexity of functions. Specifically, a function defined on manifold is said to be geodesically convex if it is convex along any geodesic. Similarly, a local minima of a geodesically convex function on manifold is also a global minima. Naturally, a question is how to distinguish the geodesically convex function.
Definition 11. Given a Riemannian manifold (M, g), a set K ⊂ M is called g-fully geodesic, if for any p, q ∈ K, any geodesic γ pq is located entirely in K.
For example, D c := {P ∈ S n ++ | det(P ) = c} with a positive constant c is not a convex set in R n×n , but is a fully geodesic set of Riemannian manifolds (S n ++ , g), where the Riemannian metric g at P is g P (U, V ) := tr(P −1 U P −1 V ). Now we present the definition of the g-geodesically convex function.
Definition 12. Given a Riemannian manifold (M, g) and a g-fully geodesic set K ⊂ M, a function f : K → R is g-geodesically convex if for any p, q ∈ K and any geodesic γ pq : [0, 1] → K connecting p, q, it holds: A g-fully geodesically convex function may not be convex. For example, f (x) := (log x) 2 , x ∈ R + is not convex in the Euclidean space, but is convex with respect to the manifold (R + , g), where g x (u, v) := ux −1 v.
Therefore, for a specific function, it is of significant importance to define a proper Riemannian metric to recognize the geodesic convexity. A natural problem is, given a manifold M and a smooth function f : M → R, whether there is a metric g such that f is geodesic convex with respective to g? It is generally not easy to prove the existence of such a metric. From the definition of the geodesic convexity, we know that if a function has a non-global local minimum, then this function is not geodesically convex for any metric. For more information on geodesic convexity, we refer to [88].

4.2.
Convergence of self-consistent field iterations. In [62,63], several classical theoretical problems from KSDFT are studied. Under certain conditions, the equivalence between KS energy minimization problems and KS equations are established. In addition, a lower bound of non-zero elements of the charge density is also analyzed. By treating the KS equation as a fixed point equation with respect to a potential function, the Jacobian matrix is explicitly derived using the spectral operator theory and the theoretical properties of the SCF method are analyzed. It is proved that the second-order derivatives of the exchange-correlation energy are uniformly bounded if the Hamiltonian has a sufficiently large eigenvalue gap. Moreover, SCF converges from any initial point and enjoys a local linear convergence rate. Related results can be found in [27,114,22,113,6,110].
Specifically, for the KS equation (2.11), we define the potential function Then, we have H ks (ρ) = H(V (ρ)). From (2.11), X are the eigenvectors corresponding to the p-smallest eigenvalues of H(V ), which is dependent on V . Then, a fixed point mapping for V can be written as where F φ (V ) = diag(X(V )X(V ) ). Therefore, each iteration of SCF is to update V k as For SCF with a simple charge-mixing strategy, the update scheme can be written as where α is an appropriate step size. Under some mild assumptions, SCF converges with a local linear convergence rate.

Pursuing global optimality.
In the Euclidean space, a common way to escape the local minimum is to add white noise to the gradient flow, which leads to a stochastic differential equation where B(t) is the standard n-by-p Brownian motion. A generalized noisy gradient flow on the Stiefel manifold is investigated in [101] where B M (t) is the Brownian motion on the manifold M := St(n, p). The construction of a Brownian motion is then given in an extrinsic form. Theoretically, it can converge to the global minima by assuming second-order continuity.

Community detection.
For community detection problems, a commonly used model is called the degree-correlated stochastic block model (DCSBM). It assumes that there are no overlaps between nodes in different communities. Specifically, the hypothesis node set [n] = {1, ..., n} contains k communities, {C * 1 , . . . , C * k } satisfying C * a ∩ C * b = ∅, ∀a = b and ∪ k a=1 C * a = [n]. In DCSBM, the network is a random graph, which can be represented by a matrix with all elements 0 to 1 represented by B ∈ S k . Let A ∈ {0, 1} n×n be the adjacency matrix of this network and A ii = 0, ∀i ∈ [n]. Then for i ∈ C * a , j ∈ C * b , i = j, A ij = 1, with probability B ab θ i θ j , 0, with probability 1 − B ab θ i θ j , where the heterogeneity of nodes is characterized by the vector θ. More specifically, larger θ i corresponds to i with more edges connecting other nodes. For DCSBM, the aforementioned relaxation model (2.15) is proposed in [106]. By solving (2.15), an approximation of the global optimal solution can be obtained with high probability.
The following theorem tells the approximation quality of an ε-approximate concave point of (2.3). Another problem with similar applications is the Z 2 synchronization problem [7]. Specifically, given noisy observations Y ij = z i z j + σW ij , where W i>j ∼ N (0, 1) and W ij = W ji , W ii = 0, we want to estimate the unknown labels z i ∈ {±1}. It can be seen as a special case of the max cut problem with p = 2. The following results are presented in [7].
Theorem 16. If σ < 1 8 √ n, then, with a high probability, all second-order stable points Q of problem (2.3) (p = 2) have the following non-trivial relationship with the true label z, i.e., for each such σ, there is ε such that 1 n Q z 2 ≥ ε.  where C ∈ S n×n is a cost matrix, A : S n×n → R m is a linear operator and A(X) = b leads to m equality constraints on X, i.e, tr(A i X) = b i with A i ∈ S n×n , b ∈ R m , i = 1, . . . , m. Define C as the constraint set C = {X ∈ S n×n : A(X) = b, X 0}.
If C is compact, it is proved in [9,71] that (4.7) has a global minimum of rank r with r(r+1) 2 ≤ m. This allows to use the Burer-Monteiro factorizations [19] (i.e., let X = Y Y with Y ∈ R n×p , p(p+1) 2 ≥ m) to solve the following non-convex optimization problem Here, we define the constraint set Since M is non-convex, there may exist many non-global local minima of (4.8). It is claimed in [20] that each local minimum of (4.8) maps to a global minimum of (4.7) if p(p+1) 2 > m. By utilizing the optimality theory of manifold optimization, any second-order stationary point can be mapped to a global minimum of (4.7) under mild assumptions [18]. Note that (4.9) is generally not a manifold. When the dimension of the space spanned by {A 1 Y, . . . , A m Y }, denoted by rank A, is fixed for all Y , M p defines a Riemannian manifold. Hence, we need the following assumptions. By comparing the optimality conditions of (4.8) and the KKT conditions of (4.7), the following equivalence between (4.7) and (4.8) is established in [18,Theorem 1.4].
Theorem 18. Let p be such that p(p+1) 2 > rank A. Suppose that Assumption 17 holds. For almost any cost matrix C ∈ S n×n , if Y ∈ M p satisfies first-and second-order necessary optimality conditions for (4.8), then Y is globally optimal and X = Y Y is globally optimal for (4.7). 4.7. Little Grothendieck problem with orthogonality constraints. Given a positive semidefinite matrix C ∈ R dn×dn , the little Grothendieck problem with orthogonality constraints can be expressed as A SDP relaxation of (4.10) is as follows [8] (4.11) max G∈R dn×dn Gii=I d×d , G 0 tr(CG).
For the original problem (4.10), a randomized approximation algorithm is presented in [8]. Specifically, it consists of the following two procedures.
• Let G be a solution to problem (4.11). Denote by the Cholesky decomposition G = LL . Let X i be a d × (nd) matrix such that L = (X 1 , X 2 , . . . , X n ) . • Let R ∈ R (nd)×d be a real-valued Gaussian random matrix whose entries are i.i.d.N (0, 1 d ). The approximate solution of the problem (4.10) can be calculated as follows where P(Y ) = arg min Z∈O d Z − Y F with Y ∈ R d×d . For the solution obtained in the above way, a constant approximation ratio on the objective function value is shown, which recovers the known 2 π approximation guarantee for the classical little Grothendieck problem.
Theorem 19. Given a symmetric matrix C 0. Let V 1 , . . . , V n ∈ O d be obtained as above. Then Z ∈ R d×d is a Gaussian random matrix whose components i.i.d. N (0, 1 d ) and σ j (Z) is the j-th singular value of Z.

5.
Conclusions. Manifold optimization has been extensively studied in literatures. We review the definition of manifold optimization, a few related applications, algorithms and analysis. However, there are still many issues and challenges. Many manifold optimization problems that can be effectively solved are still limited to relatively simple structures such as orthogonal constraints, rank constraints and etc. For other manifolds with complicated structures, what are the most efficient choices of Riemannian metrics and retraction operators are not obvious. Another interesting topic is to combine the manifold structure with the characteristics of specific problems and applications, such as graph-based data analysis, real-time data flow analysis, biomedical image analysis, etc. Nonsmooth problems appear to be more and more attractive.