# Effective Algorithms for Solving Trace Minimization Problem in Multivariate Statistics

Authors: Jiao-fen Li; Ya-qiong Wen; Xue-lin Zhou; Kai Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper develops two novel and fast Riemannian second-order approaches for solving a class of matrix trace minimization problems with orthogonality constraints, which is widely applied in multivariate statistical analysis. The existing majorization method is guaranteed to converge but its convergence rate is at best linear. A hybrid Riemannian Newton-type algorithm with both global and quadratic convergence is proposed firstly. A Riemannian trust-region method based on the proposed Newton method is further provided. Some numerical tests and application to the least squares fitting of the DEDICOM model and the orthonormal INDSCAL model are given to demonstrate the efficiency of the proposed methods. Comparisons with some latest Riemannian gradient-type methods and some existing Riemannian second-order algorithms in the MATLAB toolbox Manopt are also presented.


---

## Body

## 1. Introduction

Many multivariate statistical techniques are based on fitting a model to observed data. In most cases, these models are fitted by least squares, and as a consequence, one has to minimize a least squares loss function possibly subject to certain constraints on the parameters in the model. Least squares loss functions can often be expressed in terms of matrix trace functions. Sometimes, these have straightforward closed-form solutions, but also often such solutions are not available. In such cases, one may resort to iterative algorithms. The main purpose of the present paper is to provide several efficient algorithms that can be used for optimizing such matrix trace functions.A class of functions that would include all those of interest in multivariate statistical analysis can probably not be given, but in practice many statistical problems can be reformulated as special cases of the optimization of the general function [1–3]:(1)c+trAX+∑j=1mtrBjXCjXT,where A∈ℝp×np≤n, Bj∈ℝn×n, Cj∈ℝp×p, j=1,…,m; X is an unknown n×p matrix, c is a constant that does not depend on X, and trM means the sum of the diagonal elements of M. The matrices A, Bj, and Cj at this point are arbitrary but with the required orders. In fact, the problems of least squares fitting of a model MX (involving a matrix X of model parameters) to a data matrix Y, that is, the minimization of σX=Y−MXF2 over X (with ⋅F denoting the Frobenius norm), can often be written in a form like that of (1). To give an example, consider the case where the data matrix Y has submatrices Yk and the model MX has submatrices of the form FkXG; then,(2)σX=Y−MXF2=∑k=1mYk−FkXGF2=∑k=1mtrYkTYk−2∑k=1mtrGYkTFkX+∑k=1mtrGGTXTFkTFkX.Then, minimization of (2) is equivalent to minimization of (1) by letting c=∑k=1mtrYkTYk, A=−2∑k=1mGYkTFk, Bk=FkTFk, and Ck=GGT, k=1,…,m.The functionfX can be minimized over arbitrary X, but in many practical applications of multivariate statistics, constraints are imposed, the most frequent one being columnwise orthonormality: XTX=Ip [1–10], where Ip denotes the identity matrix of size p. Another common constraint is that of reduced rank; that is, the rank of X is forced to be less than or equal to r (r≤minn,p). The present paper discusses only the case of columnwise orthonormality. Therefore, the considered matrix trace minimization problem with constraints can be mathematically formulated by the following problem.Problem 1.
(3)minimizefX=c+trAX+∑j=1mtrBjXCjXTsubject toX∈ℝn×p,XTX=Ip,where A∈ℝp×n, Bj∈ℝn×n, and Cj∈ℝp×p, j=1,…,m are given, c is a constant, and 1≤p≤n.
Next, we present several practical applications of Problem1 in multivariate statistics. The first example was proposed in [1, 4] for a least squares fit of a model for linear dynamical systems. Following their notation, let B be a fixed n×n matrix, D, F, and H be fixed p×p matrices, X and Y be fixed n×p matrices, and Z be a variable n×p. One of their subproblems is to minimize a function(4)f1Z=Z−BZFT−XDTF2+Y−ZHTF2,over Z, subject to the constraint ZTZ=Ip. This function can be expanded as(5)f1Z=c1−2trDXTZ+2trFTDXTBZ−2trHTYTZ+trZHTHZT+trBTBZFTFZT−2trBZFTZT,where c1 denotes a constant term not depending on Z. It is readily verified that f1Z is the special case of fX with X replaced by Z, A=−2DXT+2FTDXTB−2HTYT,B1=In, B2=BTB, B3=−2B, C1=HTH, C2=FTF, and C3=FT.
The second example is the least squares fitting of the “DEDICOM model” [1, 5], which can be written as minimizing(6)f2X,R=F−XRXTF2,where F is an arbitrary fixed n×n matrix and X and R are variable matrices of orders n×r and r×r, respectively, with X columnwise orthonormal. More detailed description of this model will be given in Numerical Experiments section. Kiers et al. [6] have given an alternating least squares algorithm for fitting this model. The DEDICOM function, with R considered fixed, can be rewritten as(7)f2X=c2−2trFTXRXT,where c2 is a constant not depending on X. This is a special case of fX with A=0, B1=−2FT, C1=R.
The third example is to fit the “orthonormal INDSCAL model” (see [7], pp. 110, 118 and [1, 10]). Let Sj be a fixed symmetric matrix of order n×n, X be a variable matrix of order n×p, and Dj be a variable diagonal matrix of order p, j=1,…,m. Fitting the orthonormal INDSCAL model is equivalent to minimizing the function(8)f3X,Dj=∑j=1mSj−XDjXTF2,subject to XTX=Ip. As in Berge [8], f3 can be minimized by alternatingly updating X and Dj, j=1,…,m. Function f3, with Dj,j=1,…,m, considered fixed, can be rewritten as(9)f3X=c3−2∑j=1mtrSjXDjXT,where c3 is a constant with respect to X. This is the particular case of fX with A=0, Bj=−2Sj, Cj=Dj, j=1,…,m.
As a final example, we consider the simultaneous components analysis problem [3, 9] of minimizing(10)f4S,Ql,…,Qm=∑j=1mZj−ZjRj−1SQjTF2,over Qj and S, subject to STS=Ip. By means of an alternating least squares algorithm, they updated all parameter matrices in turn. Function f4, with Qj,j=1,…,m, considered fixed, can be expanded as(11)f4S=c4+tr−2∑jQjTS+∑jtrRj−1SQjTQjST,which, ignoring the constant c4, is a special case of fX with X replaced by S, A=−2∑jQjT, Bj=Rj−1 and Cj=QjTQjj=1,…,m.
Concerning numerical approach for solving Problem1, actually, Kiers et al. [1–3] have proposed the well-known majorization algorithm, which is based on majorizing fX by a different function gX whose minimum is easily obtained. The majorizing function is chosen such that the parameter set for which the majorizing function attains its minimum will also decrease the function to be minimized. For a more detailed description of such strategies, see [1–3]. However, since each iteration of the majorization method requires a full eigenvalue decomposition and the rate of convergence is at best linear, the method can potentially be very slow. An average of computational results of 10 random tests of the majorization method on Problem 1 is reported in Table 1, where the constant c=0 and m=2 and the given matrices A, B1, C1, B2, and C2 are all generated randomly by means of (in MATLAB style) A=randnp,n, B1=randnn,n, C1=randnp,p, B2=randnn,n, and C2=randnp,p. In Table 1, “n, p” means the dimension indexes and “CT.” and “IT.” mean the averaged total computing time in seconds and the averaged number of iterations to achieve the stopping criteria, respectively. “Grad.” means the averaged norm of Riemannian gradient gradfXk, “fdiff=fXk+1−fXk” means the averaged difference of the objective function values, and “Xfeasi” means the averaged value of feasibility XkTXk−IpF. We can observe from Table 1 and the convergence histories of the objective function value that are unreported here that the existing majorization algorithm is a monotonically convergent algorithm for minimizing the matrix trace function iteratively. However, it requires a large number of iterations to converge to within a given tolerance because the rate of convergence of this method is at best linear.
To overcome the slow linear rate of convergence, in this work we reconsider Problem1 under the framework of Riemannian optimization, by noting that the feasible set, Stn,p=X∈ℝn×pXTX=Ip, with p≤n is referred to the well-known Stiefel manifold, an np−1/2pp+1-dimensional embedded submanifold of the vector space ℝn×p. Riemannian optimization refers to optimization on Riemannian manifolds and has been extensively studied over decades. Optimization over the Stiefel manifold is an important special case of Riemannian optimization, which has recently aroused considerable research interests due to the wide applications in different fields such as the linear eigenvalue problem, the orthogonal Procrustes problem, the nearest low-rank correlation matrix problem, the Kohn–Sham total energy minimization, and singular value decomposition. Since optimization over the Stiefel manifold can be viewed as a general nonlinear optimization problem with constraints, many standard algorithms [11] in the Euclidean space can be generalized to manifold setting directly and have been explored and successfully applied to various applications, e.g., Riemannian steepest descent method [12], Riemannian curvilinear search method with Barzilai–Borwein (BB) steps [13], Riemannian Dai’s nonmonotone-based conjugate gradient method [14, 15], Riemannian Polak–Ribière–Polyak-based nonlinear conjugate gradient method [16, 17], and Riemannian Fletcher–Reeves-based conjugate gradient method [18]. However, as we know, gradient-type algorithms often perform reasonably well but might converge slowly when the generated iterates are close to an optimal solution. Usually, fast local convergence cannot be expected if only the gradient information is used. In the Euclidean space ℝn, it is well known that higher rates of convergence can be achieved by using second-order information on the cost function. The classical choice is Newton’s method; it plays a central role in the development of numerical techniques for optimization because of its simple formulation and its quadratic convergence properties. The history of Newton’s method on manifolds can be traced back to Gabay [19] who proposed a formulation for the method on embedded submanifolds of ℝn. Edelman et al. [20] proposed a formulation of Newton’s method on Stiefel manifold. Recently, Aihara and Sato [21] presented a matrix-free implementation of Riemannian Newton’s method on the Stiefel manifold. Sato [22] developed the Riemannian Newton’s method for the joint diagonalization problem on the Stiefel manifold, where the dimension of Newton’s equation is reduced and can be effectively solved by means of the Kronecker product and the vec and veck operators. Very recently, Hu et al. [23] proposed a regularized Newton method for optimization problems on Riemannian manifolds. They used a second-order approximation of the objective function in the Euclidean space to form a sequence of quadratic subproblems while keeping the manifold constraints. Numerical experiments and comparisons with other state-of-the-art methods indicated that their proposed algorithm is very promising. Another popular second-order algorithm on Stiefel manifold is the Riemannian trust-region (RTR) algorithm [24–28], which has been successfully applied to various applications. In particular, Yang et al. [29] presented the RTR algorithm for H2 model reduction of bilinear systems, where the H2 error norm is treated as a cost function on the Stiefel manifold. Sato [30] developed the RTR algorithm to the join singular value decomposition of multiple rectangular matrices, which is formulated as a Riemannian optimization problem on the product of two Stiefel manifolds.
Motivated by these works, in this paper, we are interested in extending Riemannian Newton’s method and the Riemannian trust-region method to the underlying matrix trace minimization problem. We first derive the specific expression of Riemannian gradient and Hessian of the objective function of Problem1. Specifically, we use the Kronecker product and the vectorization operators to reduce the dimension of the involved Newton’s equation and to transform the equation into a standard symmetric linear system Hx=b. The resultant system can be efficiently solved by means of direct inversion or some well-known Krylov subspace methods, such as the conjugate residual method (page 182 of [31]). Because Newton’s method is not guaranteed to have global convergence, we need to prepare a suitable starting point that is sufficiently close to an optimal solution. Here, we applied OptStiefelGBB, a state-of-the-art algorithm proposed by Wen and Yin [13] which has the benefit of requiring very little memory and has been proved to be globally convergent, to obtain a suitable initial point. The resulting hybrid Riemannian Newton-type algorithm is globally and quadratically convergent. In the Euclidean space ℝn, it is well known that the pure Newton method converges only locally, and it cannot distinguish between local minima, local maxima, and saddle points. Compared to the pure Newton-type algorithm, the advantage of the trust-region algorithm is its more stable behavior [25]. For the RTR algorithm, convergence to stationary points is guaranteed for all initial points. By utilizing the Taylor expansion on Stiefel manifold, RTR algorithm constructs a trust-region subproblem on the tangent space. Here, we proposed a new trust-region subproblem, which has a lower computational cost than the original subproblem, based on our expression of the Riemannian Hessian of the objective function. The new trust-region subproblem can be solved by the classical truncated conjugate gradient, which is most popular due to its good properties and relatively cheap computational cost.
This paper is organized as follows. In Section2, some basic geometric properties of the Stiefel manifold are given, and the representation matrix formula of Riemannian gradient and Hessian of the objective function are also derived. A hybrid Newton-type algorithm with globally and quadratically convergent is provided in Section 3. A Riemannian trust-region-based algorithm for Problem 1 is described in Section 4. Numerical experiments and application to the least squares fitting of the DEDICOM model and the orthonormal INDSCAL model and comparisons with some latest Riemannian gradient-type algorithms mentioned above and some existing second-order Riemannian algorithms in the MATALB toolbox Manopt are reported in Section 5. The conclusion is presented in Section 6.Table 1
Numerical results for the existing majorization method.
n, pCT.IT.Grad.fdiffXfeasi30, 100.59917931.005 ×10−41.561 ×10−124.108 ×10−1550, 305.32535951.001 ×10−42.436 ×10−131.221 ×10−14100, 5051.80438861.001 ×10−48.030 ×10−141.795 ×10−14100, 70131.74067921.000 ×10−44.213 ×10−142.231 ×10−14200, 50333.11041041.001 ×10−44.545 ×10−141.637 ×10−14200, 70929.85253181.000 ×10−42.913 ×10−142.232 ×10−14


## 2. Preliminaries

In this section, we first recall some notations, definitions, and basic properties of Riemannian manifolds used throughout the paper. The tangent space atx on a manifold ℳ is denoted by Txℳ. For manifolds ℳ and N and a mapping f:ℳ⟶N, the differential of f at x∈ℳ is denoted by Dfx, which is a mapping from Txℳ to Tfxℳ. Given a smooth function f on a manifold ℳ⊂ℝn1×n2, the symbol f¯ is the extension of f to the ambient Euclidean space ℝn1×n2. The symbols ∇ and grad denote the Euclidean and Riemannian gradients, respectively; i.e., given a smooth function f on a manifold ℳ⊂ℝn1×n2, ∇ and grad act on f¯ and f, respectively. The symbol Hess denotes the Riemannian Hessian. The concept of a retraction, which is a smooth map from the tangent bundle of ℳ into ℳ that approximates the exponential map to the first order, is given as follows.Definition 1.
(see [24], Definition 4.1.1). A retraction on a differentiable manifold ℳ is a smooth mapping ℛ from the tangent bundle Tℳ onto ℳ satisfying the following two conditions (here ℛx denotes the restriction of ℛ to Txℳ):(1)
ℛx0x=x, where 0x denotes the zero element of Txℳ.(2)
For anyX∈ℳ, it holds that(12)limTxℳ∈ξ⟶0ℛxξ−x+ξξ=0.
Given a retractionℛ and a smooth manifold ℳ, the general feasible algorithm framework on the manifold can be expressed as(13)xk+1=ℛxtkξk,where tk is the step size at the k-th iterate xk and ξk∈Txkℳ is a tangent vector.
We next introduce some basic geometric properties of the involved Stiefel manifoldStn,p, the reader is referred to [20, 24] for more details. The tangent space TXStn,p at X∈Stn,p can be expressed as(14)TXStn,p=ξ∈ℝn×pξTX+XTξ=0.
Since the manifoldStn,p is an embedded submanifold of the matrix Euclidean space ℝn×p, a natural metric for TXStn,p is the Euclidean metric Y,ZX=trYTZ, which induces a norm ZX=ZF, where ⋅F is the Frobenius norm of a matrix. In what follows, we denote by ⋅,⋅ and ⋅F the Riemannian metric and its induced norm on Stn,p, respectively. Under the Riemannian metric, the orthogonal projection PX at X∈Stn,p onto TXStn,p is expressed as(15)PXV=V−XsymXTV,∀V∈ℝn×p,where symM denotes the symmetric part of M, i.e., symM:=M+MT/2. We next introduce several different retraction operators on the Stiefel manifold at the current point X for a given step size t and descent direction ξ (one can refer to [27] for more details).(1)
Exponential map [20]:(16)ℛXexptξ=X,Qexpt−XTξ−RTℝ0Ip0,whereQR=−In−XXTξ is the QR decomposition of −In−XXTξ. This scheme needs to calculate an exponent of a 2p-by-2p matrix and an QR decomposition of an n-by-p matrix.(2)
Cayley transformation [13]:(17)ℛXcayleytξ=X+tPξI2p−t2QξTPξ−1QξTX,wherePξ=ΠXξ,X∈ℝn×2p and Qξ=X,−ΠXξ∈ℝn×2p with ΠX=In−1/2XXT. When p<n/2, this scheme is much cheaper than the exponential map.(3)
Polar decomposition (see [24], p. 58):(18)ℛXpolartξ=X+tξIp+t2ξTξ−1/2.The computational cost is lower than the Cayley transform but the Cayley transformation gives a better approximation to the exponential map.(4)
QR decomposition (see [24], p. 59):(19)ℛXQRtξ=qfX+tξ,whereqf⋅ denotes the Q factor of the QR decomposition of the matrix in parentheses. 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 Riemannian gradient and Hessian of an objective function are basic concepts in Riemannian optimization; we next derive the representation matrix formulas of Riemannian gradientgradf and Riemannian Hessian Hessf of the objective function in Problem 1. The Riemannian gradient, gradfX, of an objective function f at X∈Stn,p is defined to be a unique tangent vector which satisfies(20)gradfX,ξ=DfXξ,ξ∈TXStn,p.
The Hessian,HessfX, of f at X is defined to be a linear transformation of the tangent space TXStn,p through the covariant derivative ∇ξgradf of gradf evaluated at X:(21)HessfXξ≔∇ξgradf,ξ∈TXStn,p,where the covariant derivative is defined through the Levi-Civita connection ∇ on Stn,p (see [24], Definition 5.5.1).
In what follows, we definef¯X to be a function with the same form as fX defined in ℝn×p; that is,(22)f¯X=c+trAX+∑j=1mtrBjXCjXT,X∈ℝn×p.
The following two lemmas derive the explicit expressions ofgradf and Hessf.Lemma 1.
The Riemannian gradient off at X∈Stn,p can be expressed as(23)gradfX=PXAT+∑j=1mBjTXCjT+BjXCj.Proof.
SinceStn,p is a Riemannian submanifold of ℝn×p endowed with the induced metric, gradfX is equal to the orthogonal projection of the Euclidean gradient ∇f¯ at X onto TXStn,p. On the other hand, according to the matrix trace function differentiation, the Euclidean gradient ∇f¯ at X can be computed as(24)∇f¯X=AT+∑j=1mBjTXCjT+BjXCj.
Hence, by using the projectionPX given in (15), we obtain (23).Lemma 2.
Letξ be a tangent vector at X∈Stn,p, the Riemannian Hessian of f at X is expressed as a linear map on TXStn,p and given by(25)HessfXξ=PX∑j=1mBjTξCjT+BjξCj−ξsymXTAT+XT∑j=1mBjTXCjT+BjXCj.Proof.
ForX∈ℝn×p, by simple calculation, the Euclidean Hessian ∇2f¯ can be written as(26)∇2f¯Xξ=∑j=1mBjTξCjT+BjξCj.
By using the classical expression of the Riemannian connection on a Riemannian submanifold of a Euclidean space (see [32], §5.3.3) and choosing ∇ξζ≔PXDζXξ, the Riemannian Hessian hessf of f on Stn,p can also be expressed by using ∇2f¯ and the projection map (15). That is,(27)HessfXξ=PXDgradfXξ.
Note thatgradfX=PX∇f¯X. We view the right-hand side of this relation as a product of two matrix functions of X, rather than a composition of a map and a function. Then, (27) can be written as [32](28)HessfXξ=PXDgradfXξ=PXDPX∇f¯Xξ=PXDPXξ∇f¯X+PXDf¯Xξ=PXDPXξ∇f¯X+∇2f¯Xξ,in which we used the relation PXPX=PX in the last equality of (28), and for any M∈ℝn×p,(29)DPXξM=DM−XsymXTMξ=−ξsymXTM−XsymξTM.
Furthermore, ifS is p×p symmetric matrix, we have(30)PXXS=XS−XsymXTXS=XS−XS=0.
This relation can reduce the computational cost ofPXDPXξ∇f¯X in the last equality of (28). Indeed, it follows from (29) and (30) that(31)PXDPXξ∇f¯X=−PXξsymXT∇f¯X.
Consequently, by using (29) and (31), we obtain a more concrete expression for HessfXξ acting on ξ as(32)HessfXξ=PX∇2f¯Xξ−ξsymXT∇f¯X.
By substituting (24) and (26) into (32), we obtain (25).
At the end of this section, we introduce the vec and veck operators, which are useful for rewriting a matrix equation by transforming the matrix into an unknown column vector, and some useful properties will be used in the sequel. For anyM=aij∈ℝm×n, the vec operator is defined as(33)vecM=a11,…,am1,a12,…,am2,a1n,…,amnT∈ℝmn.
For anyN=bij∈ASn×n, where ASn×n means the set of all n-by-n skew-symmetric matrices, the veck operator is defined as(34)veckN=2b21,b31,…,bn1,b32,b42,…,bn2,…,bn−1n−2,bnn−2,bnn−1T∈ℝnn−1/2.
The above definitions yield the following inner product equivalence:(35)trMTM′=vecMTvecM′,∀M,M′∈ℝn×n,(36)trNTN′=veckNTveckN,∀N,N′∈ASn×n.(i)
ForM∈ℝm×n, X∈ℝn×s, and M′∈ℝs×t, one can get vecMXM′=M′T⊗MvecX, where ⊗ is the Kronecker product.(ii)
There exists anmn×mn vec-permutation matrix [33] Tm,n=∑i=1n∑j=1mEij⊗EijT such that(37)vecM=Tm,nvecMT,M∈ℝm×n,whereEij is the n×m matrix with a 1 in position i,j and zeros elsewhere. From [33], we have Tm,nTn,m=Imn. Because Tm,n is a permutation matrix and so is orthogonal, then we have Tm,n−1=Tm,nT=Tn,m, and then(38)vecMT=Tn,mvecM.(iii)
ForM∈ℝn×n, let skewM be the skew-symmetric parts of M, i.e., skewM:=M−MT/2, and we have(39)vecskewM=12In2−Tn,nvecM.(iv)
ForN∈ASn×n, we have(40)vecN=KnveckN,where the matrixKn∈ℝn2×nn−1/2 is defined as the following form:(41)Kn=12e2e3⋯en−1en0⋯00⋯0−e10⋯00e3⋯en−1en⋯00−e1⋯00−e2⋯00⋯0⋮⋮⋮⋮⋮⋮⋮⋮00⋯−e100⋯−e20⋯en00⋯0−e10⋯0−e2⋯−en−1,where ei is the ith column of In. Obviously, Kn is standard column orthogonal, that is, KnTKn=Inn−1/2. Furthermore, for N∈ASn×n, we have(42)veckN=KnTvecN.


## 3. Riemannian Newton’s Method for Problem1

Since we have already obtained the matrix expressions ofgradf and Hessf and some other requisites for Riemannian optimization algorithms, in this section, we develop Riemannian Newton’s method for Problem 1. In Newton’s method, the search direction ξk∈TXkStn,p at Xk∈Stn,p is determined to be the solution to Newton’s equation(43)HessfXkξ=−gradfXk.Afterξk is obtained, the candidate for the new iterate is then given by(44)Xk+1=ℛXkξk,where the step size is set to one for simplicity. By substituting (23) and (25) into (43), complete Riemannian Newton’s equation for Problem 1 at X∈Stn,p can be written as(45)PX∑j=1mBjTξCjT+BjξCj−ξsymXTAT+XT∑j=1mBjTXCjT+BjXCj=−PXAT+∑j=1mBjTXCjT+BjXCj.Now, the problem is how to solve the above Newton’s equation. Noting that it is complicated and difficult to solve because it must be solved forξ∈TXStn,p for given X∈Stn,p, i.e., ξ must satisfy ξTX+XTξ=0. To overcome these difficulties, we wish to obtain the representation matrix of Hess f as a linear transformation on TXStn,p for arbitrarily fixed X so that we can rewrite (45) into a standard linear equation, which permits many quite effective numerical approaches. To this end, we introduce the following lemma.Lemma 3 (see [24], p. 42).
An equivalent form of the tangent space toStn,p at X∈Stn,p is given by(46)TXStn,p=XU+X⊥VU∈ASp×p,V∈ℝn−p×p,where X⊥ is an arbitrary n×n−p matrix that satisfies XTX⊥=0 and X⊥TX⊥=In−p.From Lemma3, one can easily check that the vector space of all tangent vectors has a dimension of np−1/2pp+1. Let ξ∈TXStn,p be the search direction at X∈Stn,p. From Lemma 3, ξ can be expressed by(47)ξ=XU+X⊥V,U∈ASp×p,V∈ℝn−p×p.SinceHessfXξ is also determined to be a unique tangent vector in TXStn,p, there exist unique matrices UH∈ASp×p and VH∈ℝn−p×p which satisfy(48)HessfXξ=XUH+X⊥VH.The following proposition shows that we can writeUH and VH by using U and V.Proposition 1.
LetX∈Stn,p and X⊥∈Stn−p,p satisfy XTX⊥=0. If a tangent vector ξ∈TXStn,p is expressed as (47), then the Hessian HessfX of the objective function acts on ξ as HessfXξ=XUH+X⊥VH with(49)UH=skew∑j=1mXTBjTXUCjT+XTBjXUCj+XTBjTX⊥VCjT+XTBjX⊥VCj−UsymXTAT+XT∑j=1mBjTXCjT+BjXCj,(50)VH=∑j=1mX⊥TBjTXUCjT+X⊥TBjXUCj+X⊥TBjTX⊥VCjT+X⊥TBjX⊥VCj−VsymXTAT+XT∑j=1mBjTXCjT+BjXCj.Proof.
From (25) and (48), we have the following equality:(51)XUH+X⊥VH=PX∑j=1mBjTξCjT+BjξCj−ξsymXTAT+XT∑j=1mBjTXCjT+BjXCj.
Note thatXTPXW=skewXTW and X⊥TPXW=X⊥TW for any W∈ℝn×p. Multiplying equation (51) by XT from the left and using the relations XTX=IpandXTX⊥=0 yields(52)UH=XTPX∑j=1mBjTξCjT+BjξCj−ξsymXTAT+XT∑j=1mBjTXCjT+BjXCj=skewXT∑j=1mBjTξCjT+BjξCj−XTξsymXTAT+XT∑j=1mBjTXCjT+BjXCj.
By substituting (47) into (52) and using the relation XTξ=XTXU+X⊥V=U, we obtain the expression of UH in (49) by using U and V. Similarly, we multiply (51) by X⊥T from the left to obtain(53)VH=X⊥TPX∑j=1mBjTξCjT+BjξCj−ξsymXTAT+XT∑j=1mBjTXCjT+BjXCj=X⊥T∑j=1mBjTξCjT+BjξCj−X⊥TξsymXTAT+XT∑j=1mBjTXCjT+BjXCj.
By substituting (47) into (53) and using the relation X⊥Tξ=X⊥TXU+X⊥V=V, we obtain the expression of VH in (50).
Since the HessianHessfXξ is a linear transformation with respect to ξ=XU+X⊥V, and, therefore, with respect to the set of matrices U,V with U∈ASp×p and V∈ℝn−p×p. Meanwhile, HessfXξ is also in the tangent space TXStn,p and can be uniquely expressed as HessfXξ=XUH+X⊥VH with UH∈ASp×p,VH∈ℝn−p×p. Note that the dimension of the tangent space TXStn,p is K:=pp−1/2+pn−p; then, we can regard that the Hessian HessfXξ is a linear transformation on ℝK that transforms a K-dimension vector into another K-dimension vector. The following lemma obtains the representation matrix of that linear transformation, in which we use the vec and veck operators taking matrices and skew-symmetric matrices into vectors.Lemma 4.
LetH be a linear transformation on ℝK that acts on veckUT,vecVTT with U∈ASp×p,V∈ℝpn−p as(54)HveckUvecV=veckUHvecVH,where UH and VH are given in equations (49) and (50). Then, the representation matrix HA of H is given by(55)HA=H11H12H21H22,where(56)H11=12KpTIp2−Tp,p∑j=1mCj⊗XTBjTX+CjT⊗XTBjX−symXTAT+XT∑j=1mBjTXCjT+BjXCj⊗IpKp,H12=12KpTIp2−Tp,p∑j=1mCj⊗XTBjTX⊥+CjT⊗XTBjX⊥,H21=∑j=1mCj⊗X⊥TBjTX+CjT⊗X⊥TBjXKp,H22=∑j=1mCj⊗X⊥TBjTX⊥+CjT⊗X⊥TBjX⊥−symXTAT+XT∑j=1mBjTXCjT+BjXCj⊗In−p.Proof.
From equations (49) and (50) together with (39)–(42), veckUH and vecVH are calculated as follows:(57)veckUH=KpTvecUH=KpTvecskew∑j=1mXTBjTXUCjT+XTBjXUCj+XTBjTX⊥VCjT+XTBjX⊥VCj−UsymXTAT+XT∑j=1mBjTXCjT+BjXCj=12KpTIp2−Tp,p∑j=1mCj⊗XTBjTX+CjT⊗XTBjXvecU+Cj⊗XTBjTX⊥+CjT⊗XTBjX⊥vecV−symXTAT+XT∑j=1mBjTXCjT+BjXCj⊗IpvecU=H11veckU+H12vecV,vecVH=vec∑j=1mX⊥TBjTXUCjT+X⊥TBjXUCj+X⊥TBjTX⊥VCjT+X⊥TBjX⊥VCj−VsymXTAT+XT∑j=1mBjTXCjT+BjXCj=∑j=1mCj⊗X⊥TBjTX+CjT⊗X⊥TBjXvecU+Cj⊗X⊥TBjTX⊥+CjT⊗X⊥TBjX⊥vecV−symXTAT+XT∑j=1mBjTXCjT+BjXCj⊗In−pvecV=H21veckU+H22vecV.
This completes the proof.
The following proposition considers the structure of the representation matrixHA.Proposition 2.
The representation matrixHA is always symmetric.Proof.
Letξ and η be two tangent vectors at X which are expressed by(58)ξ=XU1+X⊥V1,η=XU2+X⊥V2,where U1,U2∈ASp×p and V1,V2∈ℝn−p×p. From the induced metric for TXStn,p, together with (35) and (36), one can get(59)ξ,η=trU1TU2+trV1TV2=veckU1TveckU2+vecV1TvecV2.
Since the HessianHessfX is symmetric with respect to the induced inner product for TXStn,p [24], we have(60)HessfXξ,η=HessfXη,ξ,ξ,η∈TXStn,p.
This is equivalent to(61)veckU1vecV1THAveckU2vecV2=HAveckU1vecV1TveckU2vecV2.
Hence, we obtain(62)HA=HATbecause ξ and η can be arbitrary tangent vectors at X.
We are now in the position to solve Newton’s equation (45). From (48), Newton’s equation (45) can be rewritten as(63)XUH+X⊥VH=−PXAT+∑j=1mBjTXCjT+BjXCj.
By multiplying (63) from the left by XT and X⊥T, (63) is equivalent to(64)UH=−XTPXAT+∑j=1mBjTXCjT+BjXCj,VH=−X⊥TPXAT+∑j=1mBjTXCjT+BjXCj.
Note again thatXTPXW=skewXTW and X⊥TPXW=X⊥TW for any W∈ℝn×p. Applying the vec operation to the two equations of (64), we get(65)veckUH=−KpTIp2−Tp,pvecXTAT+XT∑j=1mBjTXCjT+BjXCj,vecVH=−vecX⊥TAT+X⊥T∑j=1mBjTXCjT+BjXCj.
Together with Lemma4, Newton’s equation (45) can be represented in the following form:(66)HAx=g,where HA is given in (55), which is the representation symmetric matrix of Hess HessfX:(67)x=veckUvecV,(68)g=−KpTIp2−Tp,pvecXTAT+XT∑j=1mBjTXCjT+BjXCjvecX⊥TAT+X⊥T∑j=1mBjTXCjT+BjXCj.
We have thus derived an equivalent form, which is a standard symmetric linear equation, of Newton’s equation (45). After we solve this symmetric linear equation (66), we can obtain veckU and vecV; then, we can easily reshape U∈ASp×p and V∈ℝn−p×p. Therefore, we can calculate the solution ξ=XU+X⊥V of Newton’s equation (45).
Based on the above discussion, the corresponding Riemannian Newton’s method for solving Problem1 can be described as follows.
If we know a good approximate solution of the problem in advance, Newton’s method works effectively. This is because Newton’s method generates locally but quadratically convergent sequences in general (Theorem 6.3.2 of [24, 34]).Theorem 1.
(local convergence). Consider Algorithm1 with retraction ℛ as described in Section 2. Let X∗∈Stn,p be a critical point of f; gradfX∗=0. Assume that HessfX∗ is nondegenerate at X∗∈Stn,p. Then, there exists a neighborhood U of X∗ in Stn,p such that for all X0∈U, Algorithm 1 generates an infinite sequence Xk converging quadratically to X∗.From [17, 21, 22], we have several remarks for Algorithm 1 as follows.Algorithm 1: Riemannian Newton’s method for Problem 1.
(1)
Choose an initial pointX0∈Stn,p.(2)
fork=0,1,2,…do(3)
ComputeX⊥k that satisfies XkTX⊥k=0 and X⊥kTX⊥k=In−p.(4)
Compute the partition matrixHAk by Lemma 4 and the vector gk by (68). Compute xk∈ℝnp−1/2pp+1 that satisfies the symmetric linear equation HAkxk=gk(69).(5)
Partitionxk∈ℝnp−1/2pp+1 as xkT=x1kT,x2kT, where x1k∈ℝpp−1/2,x2k∈ℝn−pp.(6)
ConstructUk∈ASp×p and Vk∈ℝn−p×p that satisfy veckUk=x1k and vecVk=x2k.(7)
Computeξk=XkUk+X⊥kVk.(8)
Compute the next iterateXk+1=ℛXkξk, where ℛ is a retraction on Stn,p.(9)
end forRemark 1.
In the Step 3 of Algorithm1, we can arbitrarily fix X⊥k to satisfy XkTX⊥k=0 and X⊥kTX⊥k=In−p. Such X⊥k can be computed by the full QR decompositions. In practice, using the MATLAB’s qr function with input Xk∈Stn,p returns an orthogonal matrix Q∈ℝn×n such that Xk=QR10 where R1 is a p×p upper triangular matrix. If we partition Q=Q1,Q2 with Q1∈Stn,p and Q2∈Stn,n−p, then Q2TXk=0,In−pR10=0. Therefore, we can choose Q2 as X⊥k.Remark 2.
The symmetric linear equation (69) can be solved by direct inversion if HAk is invertible. However, if the dimension of the problem is large, it is difficult to equate by direct inversion. In this case, some Krylov subspace methods such as the conjugate residual (CR) method (page 182 of [31]) can be used to solve (69). The CR method is used to solve linear equations of the form Ax=b where A is an invertible and symmetric matrix and b is nonzero. Note that the system matrix is only required to be symmetric, not symmetric positive definite. However, we should point out that other solvers can also be applied to solve (69) because there are no specific structural requirements for the variable x. When HessfX is positive definite, the conjugate gradient (CG) method can be used. Moreover, if we consider (69) as just a linear equation of standard form, then the representation matrix HA is not necessarily symmetric. In this case, more Krylov subspace methods can be used, such as the generalized minimal residual (GMRES) method and biconjugate gradient stabilized (BiCGSTAB) method. Specifically, we use the CR method for all the numerical experiments in Section 5.Remark 3.
In general, the Newton vectorξk, solution of (43), is not necessarily a descent direction of fX. Indeed, we have(69)DfXξ=gradfX,ξ=−gradfX,HessfX−1gradfX,which is not guaranteed to be negative without additional assumptions on the operator HessfX. A sufficient condition for ξk to be a descent direction is that HessfXk be positive definite, i.e., ξ,HessfXξ>0 for all ξ≠0. The positive definiteness of the representation matrix HA can be used to check the positive definiteness of the Hessian HessfX, which can be used to check whether a critical point obtained by Algorithm 1 is a local minimum. We note that for any ξ=XU1+X⊥V1∈TXStn,p with U1∈ASp×p, V1∈ℝn−p×p,(70)HessfXξ,ξ=veckU1vecV1THAveckU1vecV1.
Therefore, if the symmetric matrixHA is positive definite, then HessfX is positive definite, and X is a local minimum. By using the partition form (55) of HA, we can easily derive a necessary condition for the positive definiteness of HA as H11 and H22 are positive definite. We can obtain from this fact an easy way to check the necessity for the positive definiteness as all the diagonal components of H11 and H22 are positive. We next consider the case where Algorithm 1 arrives at a critical point of f and assume that HA is semipositive definite. Then, HA is positive definite if and only if 0≠detHA=detH11detH22−H21H11−1H12. Therefore, under the assumption of semipositive definiteness of HA, the condition of positive definiteness of HA is necessary and sufficient for detH11≠0 and detH22−H21H11−1H12≠0.Remark 4.
Similar to unconstrained optimization in the Euclidean space,gradfXkF<ε for some constant ε>0 is a reasonable stopping criterion for Problem 1. In fact, the Lagrangian of Problem 1 is(71)ℒX,Λ=f¯X−12trΛXTX−Ip,where Λ is a symmetric matrix representing the Lagrange multiply. Then, the first-order optimality conditions in the Euclidean sense are(72)XTX=Ip,∂∂XℒX,Λ=∇f¯X−XΛ=0.
From the second equation, we haveΛ=XT∇f¯X. Since Λ must be symmetric, we further obtain Λ=∇f¯XTX. Under the conditions XTX=Ip, ∂/∂XℒX,Λ=0 is equivalent to gradfX=0, since it follows from (15), (23), and (24) that(73)PX∇f¯X=In−12XXT∇f¯X−X∇f¯XTX=In−12XXT∂∂XℒX,Λ.
Thus, first-order critical points in the Euclidean sense can be interpreted as stationary points in the Riemannian sense.
Because Newton’s method is not guaranteed to have global convergence, we need to prepare a suitable starting point that is sufficiently close to an optimal solution. Here, we apply a state-of-the-art Riemannian gradient-type algorithm, OptStiefelGBB, to obtain a suitable initial point. OptStiefelGBB, proposed by Wen and Yin [13], is simple and has the benefit of requiring very little memory. The retraction adapted in OptStiefelGBB is the Cayley transformation (17), and an initial step size computed by the Barzilai–Borwein (BB) method [35] is used to speed up the convergence. Since the BB step size does not necessarily decrease the objective value at every iteration, it may invalidate convergence, and a nonmonotone line search method based on a strategy in [36] is also used to guarantee global convergence.
To speed up the convergence, we propose a hybrid method for Problem1. The strategy is as follows: we start to perform the OptStiefelGBB until gradfXkF is sufficiently small. Then, we switch the method to Newton’s one. One possibility of the switching criterion can be set as(74)gradfXkF≤ε1,where ε1>0. The algorithm is stated as follows.
Since the cost functionfX in Problem 1 is smooth and the Stiefel manifold Stn,p is compact, OptStiefelGBB applied to Problem 1 generates a sequence converging globally to a critical point of f [13]. On the other hand, according to Theorem 1, Newton’s method generates a quadratically convergent sequence if the initial point is sufficiently close to an optimal solution. Therefore, we know that the proposed hybrid algorithm is globally and quadratically convergent if ε1 is sufficiently small.Proposition 3.
Ifε1>0 is sufficiently small, Algorithm 2 with the switching criterion (74) generates a globally and quadratically convergent sequence.Algorithm 2: Hybrid method for Problem 1.
(1)
Choose an initial pointX0∈Stn,p and a parameter ε1>0. Set k:=0.(2)
whilegradfXkF>ε1do(3)
Perform OptStiefelGBB. All parameters needed in OptStiefelGBB are the same as those in [13].(4)
k≔k+1.(5)
end while(6)
SetX0≔Xk and k≔0.(7)
Perform Steps 2–9 in Algorithm1.


## 4. Riemannian Trust-Region Method for Problem1

With the preliminary results discussed in Section2, we now state the Riemannian trust-region method for solving Problem 1. The trust-region method is an iterative method for minimizing a cost function. At each iteration step, a quadratic model of the cost function is obtained. This model is assumed to be suitable in a region (the trust region) around the current iterate. Then, an update is computed as the minimizer of the model in the trust region. The quality of the trial update is evaluated; it is consequently accepted or rejected, and the trust-region radius is adjusted.In a Euclidean spaceE, if f:E⟶ℝ and ⋅,⋅ is the inner product in E, the trust-region subproblem for finding the update ξ∈E for the current iterate xk∈E is given by(75)minξ∈Emξ=fxk+Dfxkξ+12D2fxkξ,ξs.t.ξ,ξ≤Δk2,where DFxkz denotes the directional derivative of the function F at xk in the direction of z and Δk is the trust-region radius. The quality of the model m is evaluated by means of the quotient(76)ρk=fxk−fxk+ξm0−mξ.Ifρ is close to 0 or negative, then the model is very inaccurate; i.e., the step must be rejected, and the trust-region radius must be reduced. If ρ is larger but still small, the step is accepted, and the trust-region radius is reduced. Finally, if ρ is close to 1, then there is a good correspondence between the model and the cost function; the step is accepted, and the trust-region radius can be increased.In the Riemannian trust-region method, at thek-th iteration Xk, by utilizing the Taylor expansion on manifold, we consider the following trust-region subproblem on the tangent space:(77)minξ∈TXkStn,pmXkξ=fXk+gradfXk,ξ+12HessfXkξ,ξs.t.ξ,ξ≤Δk.From (47) and (48), we have(78)HessfXkξ,ξ=XkUHk+X⊥kVHk,XkU+X⊥kV=trUHkTU+trVHkTV.On the other hand, sincegradfXk∈TXkStn,p, then from Lemma 3, gradfX can also be uniquely expressed as gradfXk=XkUgk+X⊥kVgk with Ugk=XkTgradfXk and Vgk=X⊥kTgradfXk. And(79)gradfXkξ,ξ=XkUgk+X⊥kVgk,XkU+X⊥kV=trUgkTU+trVgkTV.Then, the Riemannian trust-region subproblem (77) can be written as(80)minU,VmXkU,V=fXk+trUgkTU+VgkTV+12UHkTU+12VHkTVs.t.trUTU+trVTV≤Δk,U∈ASp×p,V∈ℝn−p×p.In order to provide a guidance for selecting the new trust-region radiusΔk , we introduce the quotient(81)ρk=fXk−fℛXkξkmXk0Xk−mXkξk,which is also used to judge the acceptance or rejection of the candidate ℛXkξk=ℛXkXkUk+X⊥kVk. Due to the fact that mXk0Xk=fXk, ρk can be given by(82)ρk=2fℛXkξk−2fXktr2UgkTUk+2VgkTVk+UHkTUk+VHkTVk.As in the Euclidean space setting, the constants1/4 and 3/4 are compared with the ratio ρk and the result determines the trust-region radius in the next iteration. Except that, the constant ρ′∈0,1/4 is used to measure ρk. The trust-region step will be taken as the next iteration if ρk>ρ′, and rejected, otherwise. Specifically, we present the Riemannian trust-region method for Problem 1 as follows.The trust-region subproblem (80) is solved iteratively and forms the inner iteration of Algorithm 3. Note that (80) is a minimization problem in ASp×p×ℝn−p×p and could then be solved by many classical approaches. A widely used approach is the truncated conjugate gradient algorithm (tCG) due to the following reasons:(1)
The tCG is a Krylov subspace-based solver in which if the initial guessU0∈ASp×p and V0∈ℝn−p×p, then the sequences Uj and Vj generated by tCG always satisfy Uj∈ASp×p and Vj∈ℝn−p×p, and hence the approximation ξj=XkUj+X⊥kVj∈TXkStn,p.(2)
The sequencemXkUj,Vj is strictly decreasing while trUjTUj+trVjTVj is strictly increasing, which makes the RTR a descent approach and guarantees the global convergence (see [[24], Proposition 7.3.2.]).Algorithm 3: Riemannian trust-region method for Problem 1.
(1)
Choose parametersΔ¯>0, Δ0∈0,Δ¯, ρ′∈0,1/4, and an initial point X0∈Stn,p.(2)
fork=0,1,2,…do(3)
Obtainξk=XkUk+X⊥kVk by solving the trust-region subproblem (80).(4)
Evaluateρk from (82).(5)
ifρk<1/4then(6)
Δk+1=1/4Δk.(7)
else ifρk>3/4 and ξkF=Δkthen(8)
Δk+1=min2Δk,Δ¯.(9)
else(10)
Δk+1=Δk.(11)
end if(12)
ifρk>ρ′then(13)
Xk+1=ℛXkξk.(14)
else(15)
Xk+1=Xk.(16)
end if(17)
end forThe pesudocode of tCG is presented in Algorithm3 for completeness. Note that we use indices in superscript without round brackets to denote the evolution of U and V within the inner iteration, while superscripts with round brackets are used in the outer iteration. As suggested by Absil et al. [24], the algorithm can stop in either after a fixed number of iterations or by the criterion:(83)trRUjTRUj+trRVjTRVj≤maxR0FminR0Fθ,κ,where κ,θ>0 are real parameters and R0F=trRU0TRU0+trRV0TRV0. In our numerical testing, we set κ=0.1 and θ=1.Remark 5.
Fixk, and let ξj∈TXkStn,p be the jth iterate of the tCG method for the subproblem expressed by (77). In the existing RTR method, which directly solves the subproblem (77), ξj is updated as a tangent vector in TXkStn,p, while in Algorithm 4, we update Uj∈ASp×p and Vj∈ℝn−p×p, where ξj=XkUj+X⊥kVj. That is, if the tCG method for (77) terminates with Uj and Vj, we only have to compute ξj=XkUj+X⊥kVj. We do not have to compute ξi for i<j. Furthermore, in the existing RTR method, hessfXkδj for some δj=XkMj+X⊥kNj∈TXkStn,p must to computed at each iteration of the tCG method. However, the proposed Algorithm 4 needs only MHj∈ASp×p and NHj∈ℝn−p×p, where hessfXkδj=XkMHj+X⊥kNHj. Therefore, if the number of iterations in the inner tCG method needed for solving the trust-region subproblems is sufficiently large, the proposed Algorithm 4 may have a shorter total computational time than the existing RTR method which directly solves the subproblem (77). These facts imply that our proposed method can reduce the computational cost.
The convergence of the general RTR-tCG on a smooth Riemannian manifold has been extensively analyzed in [24, 37] where both the global convergence and local convergence rate have been established under appropriate assumptions. Therefore, to understand the performance of the RTR-tCG method for solving Problem 1, we only need to check these assumptions, which leads to the following conclusions.Algorithm 4: The truncated CG (tCG) method for the trust-region subproblem (80).
(1)
Initialization: setU0=0∈ASp×p and V0=0∈ℝn−p×p; RU0=Ugk and RV0=Vgk; and M0=−RU0 and N0=−RV0.(2)
forj=0,1,2,…do(3)
iftrMHjTMj+trNHjTNj≤0 (negative curvature) then(4)
Computeτ such that U=Uj+τMj and V=Vj+τNj minimize mXkU,V in (80) and satisfy trUTU+trVTV=Δk.(5)
ReturnU and V.(6)
else(7)
Setαj=trRUjTRUj+trRVjTRVj/trMHjTMj+trNHjTNj.(8)
SetUj+1=Uj+αjMj and Vj+1=Vj+αjNj.(9)
iftrUj+1TUj+1+trVj+1TVj+1≥Δk (exceeded trust-region) then(10)
Computeτ≥0 such that U=Uj+τMj and V=Vj+τNj satisfy trUTU+trVTV=Δk.(11)
ReturnU and V.(12)
else(13)
SetRUj+1=RUj+αjMHj and RVj+1=RVj+αjNHj.(14)
Setβj+1=trRUj+1TRUj+1+trRVj+1TRVj+1/trRUjTRUj+trRVjTRVj.(15)
SetMj+1=−RUj+1+βj+1Mj and Nj+1=−RVj+1+βj+1Nj.(16)
end if(17)
end if(18)
end forTheorem 2 (global convergence).
LetXk be a sequence of iterates generated by Algorithms 3 and 4, and let an iterate be accepted if ρk>ρ′ with ρ′∈0,1/4. Then,(84)limk⟶∞gradfXk=0.In fact, this conclusion is a straightforward result of Theorem 4.4 and Corollary 4.6 of [37], based on the fact that the cost function fX and the adopted retraction as described in Section 2 are smooth, and the involved manifold Stn,p is a smooth and compact Riemannian manifold.Theorem 3 (local convergence speed).
Consider Algorithms3 and 4 with retraction ℛ as in Section 2 and stopping criterion in Algorithm 4 as in (83). Let X∗∈Stn,p be a (nondegenerate) local minimum of f. Then, there exists c>0 such that for all sequences Xk generated by the algorithm converging to X∗, there exists K>0 such that for all k>K,(85)distXk+1,X∗≤cdistXk,X∗minθ+1,2,with θ>0 as in (83), and “dist” defines the Riemannian distance on Stn,p (see [24], p. 46).


## 5. Numerical Experiments

In this section, we report the numerical performance of the proposed algorithms for Problem1. All the numerical experiments were completed on a personal computer with a Intel(R) Core(TM)2 Quad of 2.33 GHz CPU and 3.00 GB of RAM equipped with MATLAB R2019b. We report numerical experiments and application to the least squares fitting of the DEDICOM model and the orthonormal INDSCAL model. To illustrate the efficiency of our algorithms, we also compare Algorithms 2 and 3 with some latest Riemannian conjugate gradient- (CG-) type methods which are all applicable to Problem 1 with necessary modifications and some existing Riemannian second-order algorithms in the MATALB toolbox Manopt.


### 5.1. Algorithmic Issues

Throughout our experiments, in the implementation of Algorithm2, the parameter ε1 in (74) is set to be ε1=0.1. In the implementation of the inner iteration of Algorithm 2, the CR method to the symmetric linear equation (69), the initial matrix is all set to be 0 with suitable size at each repetition of Newton’s method, and were stopped when the relative residual 2-norms rLk2/gk2 became less than 10−6. Here, the residual is rLk=gk−HAkxLk, where xLk denotes the Lth approximation of CR at the kth repetition of Newton’s method. The largest number of CR method is set to be 1000. In the implementation of Algorithm 3, OptStiefel is also applied to obtain a suitable initial point, and the parameters are set to be Δ0=3, ρ′=0.1 and Δ¯=∞.The stopping criteria for Algorithms2 and 3 are set to be(86)gradfXkF<ε,with ε=10−6. For all the compared Riemannian gradient-type algorithms, we use the same stopping criterion as that in OptStiefelGBB [13]: we let algorithms run up to K iterations and stop it at iteration k<K if gradfXkF≤ε or relXk≤εX and relfk≤εf or(87)meanrelXk−mink,T+1,…,relXk≤10εX,meanrelfk−mink,T+1,…,relfk≤10εf,for some constants ε, εX, εf∈0,1, and T,K∈ℕ+, where(88)relXk=Xk+1−XkFp,relfk=fXk+1−fXkfXk+1.The default values ofεX, εf, T, and K are 10−6, 10−8, 5, and 2000, respectively. Throughout, “IT.,” “CT.,” “Grad.,” and “Obj.” mean the number of iterations, the total computing time in seconds, the norm of Riemannian gradient gradfXk, and the objective function value at the final iterate by implementing the proposed algorithms, respectively.


### 5.2. Numerical Experiments and Application to the Least Squares Fitting of the DEDICOM Model

DEDICOM is a model proposed by Harshman [5] for the analysis of asymmetric data. According to the DEDICOM model, a square data matrix F, containing entries Fij representing the (asymmetric) relation of object i to object j, is decomposed as(89)F=XRXT+N,where X is an n by pp<n matrix of weights for the n objects on p dimensions or aspects, R is a square matrix of order p, representing (asymmetric) relations among the p dimensions, and N is an error matrix with entries n/j representing the part of the relation of object i to object j that is not explained by the model. The objective of fitting this model to the data is to explain the data by means of relations among as small a number of dimensions as possible. These dimensions can be considered as “aspects” of the objects. The “loadings” of the objects on these aspects are given by matrix X, which is constrained to be columnwise orthonormal. The entries in matrix X indicate the importance of the aspects for the objects. The dimensionality of R and X and hence the number of aspects to be determined are to be based on some external criterion, defined by the user.The DEDICOM model has to be fit in the least squares sense over matricesX and R of order n by p and p by p, respectively. The loss function that is to be minimized can be written as(90)fX,R=F−XRXTF2,over matrix X and R, subject to the constraint XTX=Ip. Kiers and ten Berge [1, 2] have given an algorithm for this problem based on alternating least squares method that can be interpreted as alternately updating R and X. Because XTX=Ip, the minimum of f over R for fixed X is given by R=XTFX. Minimizing f over X, for fixed R, is equivalent to minimizing(91)fX=c−2trFTXRXT,subject toXTX=Ip,which is a particular case of Problem 1 with A=0, B1=−2FT, and C1=R. An algorithm based on majorization was provided in [2] for updating this subproblem on X. Because majorization algorithm requires a large number of iterations to converge within a given tolerance, we here test numerical performance of our proposed algorithms for updating X.Without loss of generality, the square data matrixF in the DEDICOM model is given by F=randnn,n and the fixed R is given by R=randp,p. The dimensions n and p were set to 500 and 5, respectively. The iterations of Algorithms 2 and 3 were started with X0=qfrandn,p or X0=qfX∗+ϑ∗randn,p, where the parameter ϑ was set to 10, 102, and 103, and X∗ is a solution of problem generated by majorization algorithm.Figure1 shows the convergence histories of the proposed Algorithm 2 after switching to Newton’s method (Algorithm 1), and Algorithm 3 after OptStiefel is performed to generate a suitable initial point for different choices of initial matrix X0. The number of iterations is plotted on the horizontal axis, and the log10 of norm of the Riemannian gradient gradfXk is plotted on the vertical axis in these figures. More numerical results are reported in Table 2, in which the term “IT.” in Algorithms 2 and 3 means the number of iterations required by OptStiefelGBB to satisfy the switching criteria and the number of extra iterations required by Algorithm 2 or Algorithm 3 to achieve the stopping criteria and “Obj.” means the objective function value −2trFTXkRXkT at the final iterate. From Figure 1 and Table 2, we observe that once Xk arrives in the appropriate region determined by ε1, Algorithms 2 and 3 generate sequences which quickly converge. We can also see that the effectiveness of Algorithms 2 and 3 was not affected by different choices of initial matrix X0, and the convergence speed of these two algorithm is not very sensitive to the choice of X0. This is in accordance with the global convergence of these two algorithms established in Proposition 3 and Theorem 2, respectively. For this reason, we use X0=qfrandn,p in the subsequent examples.Figure 1
Convergence histories of Algorithm2 after switching to Newton’s method and Algorithm 3 after OptStiefel is performed to generate a suitable initial point, with different initial matrix X0. (a) X0=qfrandn,p. (b) X0=qfX∗+10randn,p. (c) X0=qfX∗+102randn,p. (d) X0=qfX∗+103randn,p.
(a)(b)(c)(d)Table 2
Numerical results of Algorithms2 and 3 with different initial matrix X0.
CT.IT.Grad.Obj.X0=randn,pAlgorithm211.85230/45.10 ×10−10−309.19Algorithm31.15230/57.77 ×10−8−309.19X0=qfrandn,pAlgorithm29.95149/41.52 ×10−12−310.24Algorithm30.67149/52.74 ×10−7−310.24X0=qfX∗+102randn,pAlgorithm211.49189/46.58 ×10−11−309.19Algorithm31.50189/86.77 ×10−10−309.19X0=qfX∗+103randn,pAlgorithm26.41264/33.24 ×10−9−310.24Algorithm30.62264/43.36 ×10−7−310.24Next, we also show the detailed numerical results when applying the CR method to the symmetric linear equation (69) directly in the implementation of step 7 of Algorithm 2. Figure 2 shows the convergence histories of the CR method for (69). The number of iterations of the CR method is plotted on the horizontal axis, and the log10 of the relative residual norm, rLk2/gk2, is plotted on the vertical axis in the figures. More numerical results are reported in Table 3. We can observe from Figure 2 and Table 3 that the CR method is effective for solving the Newton equation (69) as the symmetric linear system, and that Newton’s method, (69), generates locally quadratically convergent sequences.Figure 2
Convergence histories of the CR method to the Newton equation (69), with different initial matrix X0 in Algorithm 2. (a) X0=qfrandn,p. (b) X0=qfX∗+10randn,p. (c) X0=qfX∗+102randn,p. (d) X0=qfX∗+103randn,p.
(a)(b)(c)(d)Table 3
Numerical results of the CR method to the symmetric linear equation (69) with different initial matrix X0 in Algorithm 2.
LCT.IT.rLk2/gk2X0=qfrandn,p11.981949.38 ×10−722.051898.86 ×10−731.991969.74 ×10−742.141939.39 ×10−7X0=qfX∗+10randn,p11.781599.08 ×10−721.631609.56 ×10−731.831699.13 ×10−741.631639.38 ×10−7X0=qfX∗+102randn,p12.022029.56 ×10−722.021899.61 ×10−732.032009.97 ×10−742.121939.69 ×10−7X0=qfX∗+103randn,p11.661629.45 ×10−721.791609.61 ×10−731.711689.16 ×10−7Finally, we numerically compare Algorithm3 in which the modified trust-region subproblems (80) are solved by Algorithm 4 with the existing trust-region method (denoted by Existing RTR) in MATLAB toolbox Manopt [38]. We fix n=500. For each p∈5,10,…,50, we compute the total computing time needed for convergence. Table 4 and Figure 3 show that the performance of the proposed Algorithm 3 is slightly better than that of the existing method in terms of computational time. This phenomenon becomes more pronounced as the matrix dimension increases.Table 4
Computational time (seconds) between Algorithm3 and the existing trust-region method.
p5101520253035404550Algorithm32.803.664.989.7320.6312.6628.2922.4635.3038.47Existing RTR3.806.3910.8322.1727.6626.0143.2549.1567.5472.38Figure 3
Computational time (seconds) of Algorithm3 and the existing trust-region method.


### 5.3. Numerical Comparison with Latest Riemannian CG Methods for the Least Squares Fitting of the Orthonormal INDSCAL Model

The orthonormal INDSCAL model is used to consider a three-way array consisting ofm symmetric n×n slices Sj [7, 8, 10]. In the model, these slices Sj are composed by the doubly centered dissimilarity matrices. The orthonormal INDSCAL model decomposes each slice as(92)Sj=XDjXT+Ej,where X is a n×r matrix (assumed to be columnwise orthonormal), Dj is a diagonal matrix, and Ej is a n×n matrix, containing the errors of the model fit. That means all slices share a common loading matrix X and differ from each other only by the diagonal elements of Dj called idiosyncratic saliences. The orthonormal INDSCAL model, appeared initially in psychometric literature, is now widely used in social sciences, marketing research, etc.The orthonormal INDSCAL problem seeks forX,D1,D2,…,Dm such that the model fits the data in least squares sense, i.e., for given and fixed n×n symmetric matrices Sj, j=1,2,…,m, minimizes the function:(93)minimizefX,Dj=∑i=1mSi−XDiXTF2subject toX,D1,D2,…,Dm∈Stn,r×Drm,where Drm=Dr×⋯×Dr and Dr denotes the set of all r×r diagonal matrices. This optimization problem has no direct analytical solution. The standard numerical solution of this problem is given by an alternating least squares algorithm [1, 8, 10]. Minimizing f over X, for fixed Dj, j=1,2,…,m, is equivalent to minimizing(94)fX=c−2∑j=1mtrSjXDjXT,subject toXTX=Ir,which is also a particular case of Problem 1 with A=0, Bj=−2Sj, Cj=Dj, j=1,…,m. A majorization-based algorithm was provided in [1] for updating this subproblem on X. Here, we test our proposed algorithms for this X subproblem.To illustrate the efficiency of our algorithms, we also compare them with the existing Riemannian trust-region algorithm (denoted by existing RTR) and the existing Riemannian BFGS algorithm (denoted by existing RBFGS) in MATLAB toolbox Manopt [38] and some latest Riemannian conjugate gradient- (CG-) type methods which are all applicable with necessary modifications. To this end, we first introduce the iterative framework of each compared Riemannian CG method. For solving Problem 1, update the iterated using the following scheme starting in X0∈Stn,p with ξ0=−gradfX0,(95)Xk+1=ℛtkξk,where tk>0 is the step size and ξk∈TXStn,p is(96)ξk+1=−gradfXk+βk+1Ttkξkξk,where T is a vector transport on Stn,p [24], which is a smooth mapping from the product of tangent bundles TXStn,p⊕TXStn,p to the tangent bundle TXStn,p, where ⊕ is the Whitney sum. There are several expressions to update the parameter βk+1 of equation (1), and some of the most popular expressions are β of the Fletcher–Reeves method:(97)βk+1FR=gradfXk+1F2gradfXkF2,and β of the Polak–Ribière method:(98)βk+1PR=gradfXk+1,gradfXk+1−TtkξkgradfXkgradfXkF2.In this experiment, we compare Algorithms2 and 3 with the following Riemannian CG methods. We note that the differences between these Riemannian CG methods are mainly due to their specific ways of updating the parameters tk and βk+1 and constructing retractions and vector transports.(i)
“RGPRP-CG”: the Riemannian geometric Polak–Ribière–Polay CG method used in [16] designed to solve Problem 1. The retraction adapted is the QR-based retraction (19), and the vector transport is constructed by using the orthogonal projection operator onto the tangent space, which has the form(99)Ttkξkξk=PℛXktkξkξk=ξk−qfXk+tkξksymqfXk+tkξkTξk.The linear search is the Armijo-type line search and the CG directionξk+1 is given by(100)ξk+1=−gradfXk+1+βk+1Tαkξkξk−θk+1Yk,where(101)Yk=gradfXk+1−TαkξkgradfXk,withβk+1 as in (3) and θk+1=gradfXk+1,Tαkξkξk/gradfXkF2.(ii)
“RFR-CG”: the Riemannian Fletcher–Reeves CG method used in [18] designed to solve Problem 1, with the same QR-based retraction and vector transport as “GPRP.” The step size is determined by tk=maxt¯ρj,j=0,1,2,… such that(102)fℛXktkξk−fXk≤δ1tkgradfXk,ξk−δ2tk2ξkF2,wheret¯ is the initial step length, δ1∈0,1, and δ2>0. The CG direction ξk+1 is given as(103)ξk+1=−gradfXk+1+βk+1Tαkξkξk−θk+1gradfXk+1,withβk+1 as in (2) and θk+1=gradfXk+1,ξk+1/gradfXkF2.(iii)
“RCG-Cayley”: the Riemannian case of Dai’s nonlinear CG method used in [14] designed to solve Problem 1 by using the Cayley transformation (17) to construct the retraction and the differentiated retraction to construct the vector transport(104)Ttkξkξk=In−tk2Wξk−2ξk,whereWξk=PξkTQξk, in which Pξk=ΠXkξk,Xk and Qξk=Xk,−ΠXkξk with ΠXk=Ip−1/2XkXkT. The step size tk is determined by the following nonmonotone line search condition instead of the Wolfe conditions:(105)fℛXktkξk≤maxfXk,…,fXk−hk+δtkgradfXk,ξk,wherehk=minh−1,k and h is some positive integer. The parameter βk+1 is computed as βk+1∈0,βk+1D with(106)βk+1D=gradfXk+1F2maxgradfXk+1,Ttkξkξk−gradfXk,ξk,−gradfXk,ξk.(iv)
“RCG-QR”: a modification of “RCG-Cayley,” by using theQR-based retraction and the corresponding differentiated retraction as a vector transport (see p. 173 of [24]).(v)
“RCG-Pol”: another modification of “RCG-Cayley” by using the polar decomposition (18) to construct a retraction and naturally using the corresponding differentiated retraction as a vector transport (see p. 173 of [24]).In our experiments, the indexm=2. We consider Si to be correlation matrices of the form(107)S1=sij1n×n=1,i=j;1i+j−1,i≠j,S2=sij2n×n=1,i=j;0.3,i−j=1;0,else.In this case, the INDSCAL can be seen as simultaneous principal component analysis of several correlation matrices [10]. We use the diagonal matrices D1=diagp,p−1,…,1 and D2=0.5∗diag1,2,…,p. All starting points X0 were generated randomly by X0=qfrandn,p. In particular, similar to [16, 18], a reasonable initial step length in the “RGPRR-CG” and “RFR-CG” methods can be set to(108)t¯=gradfXk,ξkHessfXk,ξk.All the parameters in the compared Riemannian CG methods are chosen in a standard way as in their corresponding literature.Figures4 and 5 show the convergence histories of the considered algorithms with different problem sizes n,p. The number of iterations is plotted on the horizontal axis, and the log10 of norm of the Riemannian gradient gradfXk is plotted on the vertical axis in these figures. The convergence curves generated by Algorithm 2 or Algorithm 3 are composed of two pieces, one of which is an initial part of the curve generated by OptStiefelGBB, and the second piece is the vertical line segment generated by Newton’s method. These curves show that the switching of OptStiefel to Newton’s one or trust region’s one accelerates drastically the speed of convergence. The convergence curve of OptStiefel shows that the convergence of the sequence generated by the OptStiefel is very slow. If the sequence generating method is switched to Newton’s method or trust-region method after the iterate Xk gets into the ε1-neighborhood of the target, the generated sequence reaches, in a few steps, the final point within the required accuracy. More numerical results are reported in Table 5, where the term “IT.” in Algorithms 2 and 3 is the same as that in Table 3 and “Obj.,” means the objective function value −2∑j=1mtrSjXkDjXkT at the final iterate. We can observe from Table 5 that Algorithm 3, with the implementation of OptStiefel to generate a suitable initial point, works very effectively for the orthonormal INDSCAL fitting problem and always performs much better than the existing Riemannian BFGS algorithm and all the compared Riemannian gradient-types methods in terms of total computing time and the accuracy of the solution. On the other hand, we note that the OptStiefel part of Algorithm 3 seems to be like preconditioning. In most cases, the preconditioning scheme can bring a significant reduction in terms of iteration steps than the existing Riemannian RTR method.Figure 4
Convergence histories of the considered algorithms with different problem sizesn,p.
(a)(b)(c)(d)Figure 5
Convergence histories of the considered algorithms with different problem sizesn,p.
(a)(b)(c)(d)Table 5
Comparison results of the considered algorithms.
CT.IT.Grad.Obj.n=50,p=5Algorithm20.04127/26.49 ×10−7−73.76Algorithm30.35127/36.41 ×10−11−73.76Existing-RTR0.44161.28 ×10−9−73.76Existing-RBFGS0.872818.93 ×10−7−73.76OptStiefel0.042701.27 ×10−4−73.76RCG-Cayley0.114461.83 ×10−4−73.76RCG-Pol0.8420007.87 ×10−1−73.75RCG-QR0.073691.20 ×10−4−73.76RGPRP-CG0.042331.45 ×10−4−73.76RFR-CG0.2920001.32 ×10−3−73.76n=100,p=5Algorithm20.19140/49.86 ×10−9−76.97Algorithm30.63140/55.89 ×10−11−76.97Existing-RTR0.57152.95e − 10−76.97Existing-RBFGS1.313878.63 ×10−7−76.97OptStiefel0.195011.26 ×10−4−76.97RCG-Cayley0.609852.05 ×10−4−76.97RCG-Pol2.9920001.64 ×10−1−76.95RCG-QR0.436966.43 ×10−4−76.97RGPRP-CG0.153531.32 ×10−4−76.97RFR-CG0.9320003.93 ×10−1−76.87n=200,p=5Algorithm21.74258/45.72 ×10−8−79.90Algorithm31.15258/57.54 ×10−11−79.90Existing-RTR0.82183.29 ×10−7−79.90Existing-RBFGS2.966768.33 ×10−7−79.90OptStiefel0.576481.11 ×10−4−79.90RCG-Cayley1.4510387.97 ×10−5−79.90RCG-Pol5.7016047.51 ×10−2−79.87RCG-QR1.339653.60 ×10−4−79.90RGPRP-CG0.495111.46 ×10−4−79.90RFR-CG1.9920005.32 ×10−2−79.90n=300,p=5Algorithm24.60157/45.25 ×10−8−81.55Algorithm31.55157/59.16 ×10−8−81.55Existing-RTR1.06154.85 ×10−7−81.55Existing-RBFGS3.586508.95 ×10−7−81.55OptStiefel1.779921.36 ×10−4−81.55RCG-Cayley2.097492.03 ×10−4−81.55RCG-Pol15.1420001.55 ×10−2−81.54RCG-QR1.957031.34 ×10−4−81.55RGPRP-CG0.834261.51 ×10−4−81.55RFR-CG3.8720001.74 ×10−1−81.48n=400,p=5Algorithm221.52135/135.83 ×10−8−82.68Algorithm32.59135/54.53 ×10−9−82.69Existing-RTR1.41173.07 ×10−7−82.69Existing-RBFGS8.277039.69 ×10−7−82.69OptStiefel3.415292.76 ×10−4−82.69RCG-Cayley8.217871.61 ×10−4−82.69RCG-Pol58.0120006.95 ×10−2−82.67RCG-QR14.5814201.15 ×10−3−82.69RGPRP-CG4.196261.64 ×10−4−82.69RFR-CG15.7220001.63 ×10−1−82.65n=500,p=5Algorithm228.07138/125.12 ×10−8−83.55Algorithm35.40138/82.87 ×10−8−83.55Existing-RTR4.09193.37 ×10−9−83.55Existing-RBFGS15.399171.00 ×10−6−83.55OptStiefel13.0013171.66 ×10−4−83.55RCG-Cayley29.7218292.70 ×10−4−83.55RCG-Pol82.8518383.66 ×10−2−83.54RCG-QR24.9815694.00 ×10−4−83.55RGPRP-CG11.2610511.22 ×10−4−83.55RFR-CG21.0920001.88 ×10−1−83.48n=600,p=5Algorithm256.49147/148.44 ×10−7−84.25Algorithm36.28147/94.36 ×10−8−84.25Existing-RTR2.97209.91 ×10−9−84.25Existing-RBFGS16.568009.75 ×10−7−84.25OptStiefel20.0314361.28 ×10−4−84.25RCG-Cayley32.2114272.60 ×10−4−84.25RCG-Pol125.5320003.98 ×10−1−84.24RCG-QR24.0810672.53 ×10−4−84.25RGPRP-CG17.4711501.47 ×10−4−84.25RFR-CG30.4620004.96 ×10−1−84.12n=700,p=5Algorithm267.50143/114.94 ×10−7−84.83Algorithm35.08143/69.72 ×10−7−84.84Existing-RTR3.36173.83 ×10−7−84.84Existing-RBFGS16.686091.04 ×10−6−84.84OptStiefel13.637201.22 ×10−4−84.84RCG-Cayley25.078391.57 ×10−4−84.84RCG-Pol171.9620001.53 ×10−1−84.82RCG-QR22.827934.87 ×10−5−84.84RGPRP-CG14.487171.35 ×10−4−84.84RFR-CG41.2120009.98 ×10−2−84.81We can also see from Table5 that when the problem size is large, Algorithm 2 uses much more computing time than that of other methods. This is because Newton’s method involve an inner iteration, the CR method to the symmetric linear equation (69), it consumes a lot of time to get an approximate solution to the equation, especially when the system dimension is large. Even though Newton’s method takes the longest time per iteration among the considered algorithms, the convergence is very quick as a whole.


## 5.1. Algorithmic Issues

Throughout our experiments, in the implementation of Algorithm2, the parameter ε1 in (74) is set to be ε1=0.1. In the implementation of the inner iteration of Algorithm 2, the CR method to the symmetric linear equation (69), the initial matrix is all set to be 0 with suitable size at each repetition of Newton’s method, and were stopped when the relative residual 2-norms rLk2/gk2 became less than 10−6. Here, the residual is rLk=gk−HAkxLk, where xLk denotes the Lth approximation of CR at the kth repetition of Newton’s method. The largest number of CR method is set to be 1000. In the implementation of Algorithm 3, OptStiefel is also applied to obtain a suitable initial point, and the parameters are set to be Δ0=3, ρ′=0.1 and Δ¯=∞.The stopping criteria for Algorithms2 and 3 are set to be(86)gradfXkF<ε,with ε=10−6. For all the compared Riemannian gradient-type algorithms, we use the same stopping criterion as that in OptStiefelGBB [13]: we let algorithms run up to K iterations and stop it at iteration k<K if gradfXkF≤ε or relXk≤εX and relfk≤εf or(87)meanrelXk−mink,T+1,…,relXk≤10εX,meanrelfk−mink,T+1,…,relfk≤10εf,for some constants ε, εX, εf∈0,1, and T,K∈ℕ+, where(88)relXk=Xk+1−XkFp,relfk=fXk+1−fXkfXk+1.The default values ofεX, εf, T, and K are 10−6, 10−8, 5, and 2000, respectively. Throughout, “IT.,” “CT.,” “Grad.,” and “Obj.” mean the number of iterations, the total computing time in seconds, the norm of Riemannian gradient gradfXk, and the objective function value at the final iterate by implementing the proposed algorithms, respectively.


## 5.2. Numerical Experiments and Application to the Least Squares Fitting of the DEDICOM Model

DEDICOM is a model proposed by Harshman [5] for the analysis of asymmetric data. According to the DEDICOM model, a square data matrix F, containing entries Fij representing the (asymmetric) relation of object i to object j, is decomposed as(89)F=XRXT+N,where X is an n by pp<n matrix of weights for the n objects on p dimensions or aspects, R is a square matrix of order p, representing (asymmetric) relations among the p dimensions, and N is an error matrix with entries n/j representing the part of the relation of object i to object j that is not explained by the model. The objective of fitting this model to the data is to explain the data by means of relations among as small a number of dimensions as possible. These dimensions can be considered as “aspects” of the objects. The “loadings” of the objects on these aspects are given by matrix X, which is constrained to be columnwise orthonormal. The entries in matrix X indicate the importance of the aspects for the objects. The dimensionality of R and X and hence the number of aspects to be determined are to be based on some external criterion, defined by the user.The DEDICOM model has to be fit in the least squares sense over matricesX and R of order n by p and p by p, respectively. The loss function that is to be minimized can be written as(90)fX,R=F−XRXTF2,over matrix X and R, subject to the constraint XTX=Ip. Kiers and ten Berge [1, 2] have given an algorithm for this problem based on alternating least squares method that can be interpreted as alternately updating R and X. Because XTX=Ip, the minimum of f over R for fixed X is given by R=XTFX. Minimizing f over X, for fixed R, is equivalent to minimizing(91)fX=c−2trFTXRXT,subject toXTX=Ip,which is a particular case of Problem 1 with A=0, B1=−2FT, and C1=R. An algorithm based on majorization was provided in [2] for updating this subproblem on X. Because majorization algorithm requires a large number of iterations to converge within a given tolerance, we here test numerical performance of our proposed algorithms for updating X.Without loss of generality, the square data matrixF in the DEDICOM model is given by F=randnn,n and the fixed R is given by R=randp,p. The dimensions n and p were set to 500 and 5, respectively. The iterations of Algorithms 2 and 3 were started with X0=qfrandn,p or X0=qfX∗+ϑ∗randn,p, where the parameter ϑ was set to 10, 102, and 103, and X∗ is a solution of problem generated by majorization algorithm.Figure1 shows the convergence histories of the proposed Algorithm 2 after switching to Newton’s method (Algorithm 1), and Algorithm 3 after OptStiefel is performed to generate a suitable initial point for different choices of initial matrix X0. The number of iterations is plotted on the horizontal axis, and the log10 of norm of the Riemannian gradient gradfXk is plotted on the vertical axis in these figures. More numerical results are reported in Table 2, in which the term “IT.” in Algorithms 2 and 3 means the number of iterations required by OptStiefelGBB to satisfy the switching criteria and the number of extra iterations required by Algorithm 2 or Algorithm 3 to achieve the stopping criteria and “Obj.” means the objective function value −2trFTXkRXkT at the final iterate. From Figure 1 and Table 2, we observe that once Xk arrives in the appropriate region determined by ε1, Algorithms 2 and 3 generate sequences which quickly converge. We can also see that the effectiveness of Algorithms 2 and 3 was not affected by different choices of initial matrix X0, and the convergence speed of these two algorithm is not very sensitive to the choice of X0. This is in accordance with the global convergence of these two algorithms established in Proposition 3 and Theorem 2, respectively. For this reason, we use X0=qfrandn,p in the subsequent examples.Figure 1
Convergence histories of Algorithm2 after switching to Newton’s method and Algorithm 3 after OptStiefel is performed to generate a suitable initial point, with different initial matrix X0. (a) X0=qfrandn,p. (b) X0=qfX∗+10randn,p. (c) X0=qfX∗+102randn,p. (d) X0=qfX∗+103randn,p.
(a)(b)(c)(d)Table 2
Numerical results of Algorithms2 and 3 with different initial matrix X0.
CT.IT.Grad.Obj.X0=randn,pAlgorithm211.85230/45.10 ×10−10−309.19Algorithm31.15230/57.77 ×10−8−309.19X0=qfrandn,pAlgorithm29.95149/41.52 ×10−12−310.24Algorithm30.67149/52.74 ×10−7−310.24X0=qfX∗+102randn,pAlgorithm211.49189/46.58 ×10−11−309.19Algorithm31.50189/86.77 ×10−10−309.19X0=qfX∗+103randn,pAlgorithm26.41264/33.24 ×10−9−310.24Algorithm30.62264/43.36 ×10−7−310.24Next, we also show the detailed numerical results when applying the CR method to the symmetric linear equation (69) directly in the implementation of step 7 of Algorithm 2. Figure 2 shows the convergence histories of the CR method for (69). The number of iterations of the CR method is plotted on the horizontal axis, and the log10 of the relative residual norm, rLk2/gk2, is plotted on the vertical axis in the figures. More numerical results are reported in Table 3. We can observe from Figure 2 and Table 3 that the CR method is effective for solving the Newton equation (69) as the symmetric linear system, and that Newton’s method, (69), generates locally quadratically convergent sequences.Figure 2
Convergence histories of the CR method to the Newton equation (69), with different initial matrix X0 in Algorithm 2. (a) X0=qfrandn,p. (b) X0=qfX∗+10randn,p. (c) X0=qfX∗+102randn,p. (d) X0=qfX∗+103randn,p.
(a)(b)(c)(d)Table 3
Numerical results of the CR method to the symmetric linear equation (69) with different initial matrix X0 in Algorithm 2.
LCT.IT.rLk2/gk2X0=qfrandn,p11.981949.38 ×10−722.051898.86 ×10−731.991969.74 ×10−742.141939.39 ×10−7X0=qfX∗+10randn,p11.781599.08 ×10−721.631609.56 ×10−731.831699.13 ×10−741.631639.38 ×10−7X0=qfX∗+102randn,p12.022029.56 ×10−722.021899.61 ×10−732.032009.97 ×10−742.121939.69 ×10−7X0=qfX∗+103randn,p11.661629.45 ×10−721.791609.61 ×10−731.711689.16 ×10−7Finally, we numerically compare Algorithm3 in which the modified trust-region subproblems (80) are solved by Algorithm 4 with the existing trust-region method (denoted by Existing RTR) in MATLAB toolbox Manopt [38]. We fix n=500. For each p∈5,10,…,50, we compute the total computing time needed for convergence. Table 4 and Figure 3 show that the performance of the proposed Algorithm 3 is slightly better than that of the existing method in terms of computational time. This phenomenon becomes more pronounced as the matrix dimension increases.Table 4
Computational time (seconds) between Algorithm3 and the existing trust-region method.
p5101520253035404550Algorithm32.803.664.989.7320.6312.6628.2922.4635.3038.47Existing RTR3.806.3910.8322.1727.6626.0143.2549.1567.5472.38Figure 3
Computational time (seconds) of Algorithm3 and the existing trust-region method.


## 5.3. Numerical Comparison with Latest Riemannian CG Methods for the Least Squares Fitting of the Orthonormal INDSCAL Model

The orthonormal INDSCAL model is used to consider a three-way array consisting ofm symmetric n×n slices Sj [7, 8, 10]. In the model, these slices Sj are composed by the doubly centered dissimilarity matrices. The orthonormal INDSCAL model decomposes each slice as(92)Sj=XDjXT+Ej,where X is a n×r matrix (assumed to be columnwise orthonormal), Dj is a diagonal matrix, and Ej is a n×n matrix, containing the errors of the model fit. That means all slices share a common loading matrix X and differ from each other only by the diagonal elements of Dj called idiosyncratic saliences. The orthonormal INDSCAL model, appeared initially in psychometric literature, is now widely used in social sciences, marketing research, etc.The orthonormal INDSCAL problem seeks forX,D1,D2,…,Dm such that the model fits the data in least squares sense, i.e., for given and fixed n×n symmetric matrices Sj, j=1,2,…,m, minimizes the function:(93)minimizefX,Dj=∑i=1mSi−XDiXTF2subject toX,D1,D2,…,Dm∈Stn,r×Drm,where Drm=Dr×⋯×Dr and Dr denotes the set of all r×r diagonal matrices. This optimization problem has no direct analytical solution. The standard numerical solution of this problem is given by an alternating least squares algorithm [1, 8, 10]. Minimizing f over X, for fixed Dj, j=1,2,…,m, is equivalent to minimizing(94)fX=c−2∑j=1mtrSjXDjXT,subject toXTX=Ir,which is also a particular case of Problem 1 with A=0, Bj=−2Sj, Cj=Dj, j=1,…,m. A majorization-based algorithm was provided in [1] for updating this subproblem on X. Here, we test our proposed algorithms for this X subproblem.To illustrate the efficiency of our algorithms, we also compare them with the existing Riemannian trust-region algorithm (denoted by existing RTR) and the existing Riemannian BFGS algorithm (denoted by existing RBFGS) in MATLAB toolbox Manopt [38] and some latest Riemannian conjugate gradient- (CG-) type methods which are all applicable with necessary modifications. To this end, we first introduce the iterative framework of each compared Riemannian CG method. For solving Problem 1, update the iterated using the following scheme starting in X0∈Stn,p with ξ0=−gradfX0,(95)Xk+1=ℛtkξk,where tk>0 is the step size and ξk∈TXStn,p is(96)ξk+1=−gradfXk+βk+1Ttkξkξk,where T is a vector transport on Stn,p [24], which is a smooth mapping from the product of tangent bundles TXStn,p⊕TXStn,p to the tangent bundle TXStn,p, where ⊕ is the Whitney sum. There are several expressions to update the parameter βk+1 of equation (1), and some of the most popular expressions are β of the Fletcher–Reeves method:(97)βk+1FR=gradfXk+1F2gradfXkF2,and β of the Polak–Ribière method:(98)βk+1PR=gradfXk+1,gradfXk+1−TtkξkgradfXkgradfXkF2.In this experiment, we compare Algorithms2 and 3 with the following Riemannian CG methods. We note that the differences between these Riemannian CG methods are mainly due to their specific ways of updating the parameters tk and βk+1 and constructing retractions and vector transports.(i)
“RGPRP-CG”: the Riemannian geometric Polak–Ribière–Polay CG method used in [16] designed to solve Problem 1. The retraction adapted is the QR-based retraction (19), and the vector transport is constructed by using the orthogonal projection operator onto the tangent space, which has the form(99)Ttkξkξk=PℛXktkξkξk=ξk−qfXk+tkξksymqfXk+tkξkTξk.The linear search is the Armijo-type line search and the CG directionξk+1 is given by(100)ξk+1=−gradfXk+1+βk+1Tαkξkξk−θk+1Yk,where(101)Yk=gradfXk+1−TαkξkgradfXk,withβk+1 as in (3) and θk+1=gradfXk+1,Tαkξkξk/gradfXkF2.(ii)
“RFR-CG”: the Riemannian Fletcher–Reeves CG method used in [18] designed to solve Problem 1, with the same QR-based retraction and vector transport as “GPRP.” The step size is determined by tk=maxt¯ρj,j=0,1,2,… such that(102)fℛXktkξk−fXk≤δ1tkgradfXk,ξk−δ2tk2ξkF2,wheret¯ is the initial step length, δ1∈0,1, and δ2>0. The CG direction ξk+1 is given as(103)ξk+1=−gradfXk+1+βk+1Tαkξkξk−θk+1gradfXk+1,withβk+1 as in (2) and θk+1=gradfXk+1,ξk+1/gradfXkF2.(iii)
“RCG-Cayley”: the Riemannian case of Dai’s nonlinear CG method used in [14] designed to solve Problem 1 by using the Cayley transformation (17) to construct the retraction and the differentiated retraction to construct the vector transport(104)Ttkξkξk=In−tk2Wξk−2ξk,whereWξk=PξkTQξk, in which Pξk=ΠXkξk,Xk and Qξk=Xk,−ΠXkξk with ΠXk=Ip−1/2XkXkT. The step size tk is determined by the following nonmonotone line search condition instead of the Wolfe conditions:(105)fℛXktkξk≤maxfXk,…,fXk−hk+δtkgradfXk,ξk,wherehk=minh−1,k and h is some positive integer. The parameter βk+1 is computed as βk+1∈0,βk+1D with(106)βk+1D=gradfXk+1F2maxgradfXk+1,Ttkξkξk−gradfXk,ξk,−gradfXk,ξk.(iv)
“RCG-QR”: a modification of “RCG-Cayley,” by using theQR-based retraction and the corresponding differentiated retraction as a vector transport (see p. 173 of [24]).(v)
“RCG-Pol”: another modification of “RCG-Cayley” by using the polar decomposition (18) to construct a retraction and naturally using the corresponding differentiated retraction as a vector transport (see p. 173 of [24]).In our experiments, the indexm=2. We consider Si to be correlation matrices of the form(107)S1=sij1n×n=1,i=j;1i+j−1,i≠j,S2=sij2n×n=1,i=j;0.3,i−j=1;0,else.In this case, the INDSCAL can be seen as simultaneous principal component analysis of several correlation matrices [10]. We use the diagonal matrices D1=diagp,p−1,…,1 and D2=0.5∗diag1,2,…,p. All starting points X0 were generated randomly by X0=qfrandn,p. In particular, similar to [16, 18], a reasonable initial step length in the “RGPRR-CG” and “RFR-CG” methods can be set to(108)t¯=gradfXk,ξkHessfXk,ξk.All the parameters in the compared Riemannian CG methods are chosen in a standard way as in their corresponding literature.Figures4 and 5 show the convergence histories of the considered algorithms with different problem sizes n,p. The number of iterations is plotted on the horizontal axis, and the log10 of norm of the Riemannian gradient gradfXk is plotted on the vertical axis in these figures. The convergence curves generated by Algorithm 2 or Algorithm 3 are composed of two pieces, one of which is an initial part of the curve generated by OptStiefelGBB, and the second piece is the vertical line segment generated by Newton’s method. These curves show that the switching of OptStiefel to Newton’s one or trust region’s one accelerates drastically the speed of convergence. The convergence curve of OptStiefel shows that the convergence of the sequence generated by the OptStiefel is very slow. If the sequence generating method is switched to Newton’s method or trust-region method after the iterate Xk gets into the ε1-neighborhood of the target, the generated sequence reaches, in a few steps, the final point within the required accuracy. More numerical results are reported in Table 5, where the term “IT.” in Algorithms 2 and 3 is the same as that in Table 3 and “Obj.,” means the objective function value −2∑j=1mtrSjXkDjXkT at the final iterate. We can observe from Table 5 that Algorithm 3, with the implementation of OptStiefel to generate a suitable initial point, works very effectively for the orthonormal INDSCAL fitting problem and always performs much better than the existing Riemannian BFGS algorithm and all the compared Riemannian gradient-types methods in terms of total computing time and the accuracy of the solution. On the other hand, we note that the OptStiefel part of Algorithm 3 seems to be like preconditioning. In most cases, the preconditioning scheme can bring a significant reduction in terms of iteration steps than the existing Riemannian RTR method.Figure 4
Convergence histories of the considered algorithms with different problem sizesn,p.
(a)(b)(c)(d)Figure 5
Convergence histories of the considered algorithms with different problem sizesn,p.
(a)(b)(c)(d)Table 5
Comparison results of the considered algorithms.
CT.IT.Grad.Obj.n=50,p=5Algorithm20.04127/26.49 ×10−7−73.76Algorithm30.35127/36.41 ×10−11−73.76Existing-RTR0.44161.28 ×10−9−73.76Existing-RBFGS0.872818.93 ×10−7−73.76OptStiefel0.042701.27 ×10−4−73.76RCG-Cayley0.114461.83 ×10−4−73.76RCG-Pol0.8420007.87 ×10−1−73.75RCG-QR0.073691.20 ×10−4−73.76RGPRP-CG0.042331.45 ×10−4−73.76RFR-CG0.2920001.32 ×10−3−73.76n=100,p=5Algorithm20.19140/49.86 ×10−9−76.97Algorithm30.63140/55.89 ×10−11−76.97Existing-RTR0.57152.95e − 10−76.97Existing-RBFGS1.313878.63 ×10−7−76.97OptStiefel0.195011.26 ×10−4−76.97RCG-Cayley0.609852.05 ×10−4−76.97RCG-Pol2.9920001.64 ×10−1−76.95RCG-QR0.436966.43 ×10−4−76.97RGPRP-CG0.153531.32 ×10−4−76.97RFR-CG0.9320003.93 ×10−1−76.87n=200,p=5Algorithm21.74258/45.72 ×10−8−79.90Algorithm31.15258/57.54 ×10−11−79.90Existing-RTR0.82183.29 ×10−7−79.90Existing-RBFGS2.966768.33 ×10−7−79.90OptStiefel0.576481.11 ×10−4−79.90RCG-Cayley1.4510387.97 ×10−5−79.90RCG-Pol5.7016047.51 ×10−2−79.87RCG-QR1.339653.60 ×10−4−79.90RGPRP-CG0.495111.46 ×10−4−79.90RFR-CG1.9920005.32 ×10−2−79.90n=300,p=5Algorithm24.60157/45.25 ×10−8−81.55Algorithm31.55157/59.16 ×10−8−81.55Existing-RTR1.06154.85 ×10−7−81.55Existing-RBFGS3.586508.95 ×10−7−81.55OptStiefel1.779921.36 ×10−4−81.55RCG-Cayley2.097492.03 ×10−4−81.55RCG-Pol15.1420001.55 ×10−2−81.54RCG-QR1.957031.34 ×10−4−81.55RGPRP-CG0.834261.51 ×10−4−81.55RFR-CG3.8720001.74 ×10−1−81.48n=400,p=5Algorithm221.52135/135.83 ×10−8−82.68Algorithm32.59135/54.53 ×10−9−82.69Existing-RTR1.41173.07 ×10−7−82.69Existing-RBFGS8.277039.69 ×10−7−82.69OptStiefel3.415292.76 ×10−4−82.69RCG-Cayley8.217871.61 ×10−4−82.69RCG-Pol58.0120006.95 ×10−2−82.67RCG-QR14.5814201.15 ×10−3−82.69RGPRP-CG4.196261.64 ×10−4−82.69RFR-CG15.7220001.63 ×10−1−82.65n=500,p=5Algorithm228.07138/125.12 ×10−8−83.55Algorithm35.40138/82.87 ×10−8−83.55Existing-RTR4.09193.37 ×10−9−83.55Existing-RBFGS15.399171.00 ×10−6−83.55OptStiefel13.0013171.66 ×10−4−83.55RCG-Cayley29.7218292.70 ×10−4−83.55RCG-Pol82.8518383.66 ×10−2−83.54RCG-QR24.9815694.00 ×10−4−83.55RGPRP-CG11.2610511.22 ×10−4−83.55RFR-CG21.0920001.88 ×10−1−83.48n=600,p=5Algorithm256.49147/148.44 ×10−7−84.25Algorithm36.28147/94.36 ×10−8−84.25Existing-RTR2.97209.91 ×10−9−84.25Existing-RBFGS16.568009.75 ×10−7−84.25OptStiefel20.0314361.28 ×10−4−84.25RCG-Cayley32.2114272.60 ×10−4−84.25RCG-Pol125.5320003.98 ×10−1−84.24RCG-QR24.0810672.53 ×10−4−84.25RGPRP-CG17.4711501.47 ×10−4−84.25RFR-CG30.4620004.96 ×10−1−84.12n=700,p=5Algorithm267.50143/114.94 ×10−7−84.83Algorithm35.08143/69.72 ×10−7−84.84Existing-RTR3.36173.83 ×10−7−84.84Existing-RBFGS16.686091.04 ×10−6−84.84OptStiefel13.637201.22 ×10−4−84.84RCG-Cayley25.078391.57 ×10−4−84.84RCG-Pol171.9620001.53 ×10−1−84.82RCG-QR22.827934.87 ×10−5−84.84RGPRP-CG14.487171.35 ×10−4−84.84RFR-CG41.2120009.98 ×10−2−84.81We can also see from Table5 that when the problem size is large, Algorithm 2 uses much more computing time than that of other methods. This is because Newton’s method involve an inner iteration, the CR method to the symmetric linear equation (69), it consumes a lot of time to get an approximate solution to the equation, especially when the system dimension is large. Even though Newton’s method takes the longest time per iteration among the considered algorithms, the convergence is very quick as a whole.


## 6. Conclusion

We have dealt with the trace function minimization problem with orthogonal constraints as an optimization problem on the manifoldStn,p and have developed Riemannian Newton’s method and Riemannian trust-region method for the problem. For diminishing the difficulty in solving Newton’s equation in Newton’s method, we have computed the representation matrix of the Hessian of the objective function using the vectorization operators and have reduced Newton’s equation into a standard symmetric linear equation with dimension reduction. In addition, we have proposed a new trust-region method in which trust-region subproblems are solved by the truncated conjugate gradient method based on our expressions of the Hessian of the objective function. Furthermore, we have performed numerical experiments to verify that our present algorithms almost always perform much better than some Riemannian gradient-type methods.



---

*Source: 3054764-2020-08-10.xml*