{'fraction_non_alphanumeric': 0.04764056224899598, 'fraction_numerical': 0.016039156626506023, 'mean_word_length': 4.3665140086206895, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parallelizable Algorithms for Optimization Problems with Orthogonality Constraints

To construct a parallel approach for solving optimization problems with orthogonality constraints is usually regarded as an extremely difficult mission, due to the low scalability of the orthonormalization procedure. However, such demand is particularly huge in some application areas such as materials computation. In this paper, we propose a proximal linearized augmented Lagrangian algorithm (PLAM) for solving optimization problems with orthogonality constraints. Unlike the classical augmented Lagrangian methods, in our algorithm, the prime variables are updated by minimizing a proximal linearized approximation of the augmented Lagrangian function, meanwhile the dual variables are updated by a closed-form expression which holds at any first-order stationary point. The orthonormalization procedure is only invoked once at the last step of the above mentioned algorithm if high-precision feasibility is needed. Consequently, the main parts of the proposed algorithm can be parallelized naturally. We establish global subsequence convergence, worst-case complexity and local convergence rate for PLAM under some mild assumptions. To reduce the sensitivity of the penalty parameter, we put forward a modification of PLAM, which is called parallelizable column-wise block minimization of PLAM (PCAL). Numerical experiments in serial illustrate that the novel updating rule for the Lagrangian multipliers significantly accelerates the convergence of PLAM and makes it comparable with the existent feasible solvers for optimization problems with orthogonality constraints, and the performance of PCAL does not highly rely on the choice of the penalty parameter. Numerical experiments under parallel environment demonstrate that PCAL attains good performance and high scalability in solving discretized Kohn-Sham total energy minimization problems.

1. Introduction. In this paper, we consider the following matrix variable optimization problem with orthogonality constraints. where I p is the p-by-p identity matrix with 2p ≤ n, and f : R n×p −→ R is a continuously differentiable function. The feasible set of the orthogonality constraints is also known as Stiefel manifold, S n,p = {X ∈ R n×p | X ⊤ X = I p }.
Throughout this paper, we assume ASSUMPTION 1.1. (Blanket Assumption) f is continuously differentiable.
The twice differentiability of f will be particularly mentioned once it is required in some theoretical analyses.
The numerical experiments under serial computing demonstrate the way how to choose default settings for our algorithms, and show that the infeasible algorithms are at least as efficient as the existent feasible algorithms in solving a bunch of test problems. The numerical experiments under parallel computing illustrate the computational complexity of PCAL and expose its high scalability.
1.3. Organization and Notations. The motivation of new approaches will be introduced in the next section. In Section 3, we will present the algorithm frameworks. We investigated the theoretical behaviors of the new proposed algorithms in Section 4. Numerical experiments will be demonstrated in Section 5. In the last section, we will draw a brief conclusion and discuss possible future works.
Notations. S p := {X ∈ R p×p | X ⊤ = X} refers to the p-by-p real symmetric matrices set. λ max (A) and λ min (A) stand for the largest and smallest eigenvalues of given symmetric real matrix A, respectively. σ max (A) and σ min (A) denote the largest and smallest singular values of given real matrix A, respectively. X † := (X ⊤ X) −1 X ⊤ refers to the pseudo inverse of X. Diag(v) ∈ S n denotes a diagonal matrix with all entries of v ∈ R n in its diagonal, and diag(A) ∈ R n extracts the diagonal entries of matrix A ∈ R n×n . For convenience, Φ(M ) := Diag(diag(M )) represents the diagonal matrix with the diagonal entries of square matrix M in its diagonal. Ψ (A) := 1 2 (A + A ⊤ ) stands for the average of a square matrix and its transpose.

Motivation.
As mentioned in the previous section, almost all the existing practically useful methods require feasible iterates all the time. To realize feasibility, either explicit or implicit orthonormalization requires to be invoked. Such kind of calculation lacks of scalability and hence becomes the bottleneck computation in the corresponding algorithms. For example, we consider the discretized Kohn-Sham total energy minimization (1.2). In each iteration, the function value and first-order derivative evaluation cost O(n log n + np) or O(np) flops per iteration, depending on whether plane wave or finite difference, respectively, is used in the discretization scheme. For the main iteration of any algorithm for solving (1.2) developed in recent decade, the computational cost per iteration is O(np 2 ) for BLAS3 calculation, plus O(p 3 ) for orthonormalization which can hardly be parallelized.
To break through this bottleneck, we suggest to use infeasible methods to take the place of feasible methods.
There is no existent infeasible approach for general purpose reported to be efficient for optimization problems with orthogonality constraints. Previous infeasible approaches designed for (1.1) either work specially for Rayleigh-Ritz trace minimization [23,33], or adopt ADMM framework after introducing auxiliary variables to split the objective and orthogonality constraints [17]. The previous ones can hardly be extended to general objective, while the latter ones does not have good performance in general.
In the following subsections, we introduce how we come up with a new idea on constructing an efficient infeasible algorithm for problem (1.1).
2.1. The Optimality Condition. We start from the optimality condition of the optimization problem with orthogonality constraints (1.1). The first-order optimality condition of problem (1.1) can be written as ∇f (X) = XΛ; X ⊤ X = I p , (2.1) where Λ ∈ S p consists of the Lagrangian multipliers of the orthogonality constraints. Condition (2.1) has the following equivalent form where Λ is eliminated.
∇f (X) − X∇f (X) ⊤ X = 0; X ⊤ X = I p . (2.2) 3 DEFINITION 2.1. We call X a first-order stationary point, if condition (2.2) holds. We call X a second-order stationary point, if it is a first-order stationary point and satisfies where T (X) := {Y | Y ⊤ X + X ⊤ Y = 0} is the tangent space of the orthogonality constraints at X.
The following proposition can be easily verified and hence its proof is omitted here. PROPOSITION 2.2. If X is a local minimizer of (1.1), it has to be a second-order stationary point. If X is a strict local minimizer 3 , if and only if X is a first-order stationary point and satisfies

Augmented Lagrangian Method.
A straightforward idea to solve (1.1) without requiring feasibility in each iteration is to employ the Augmented Lagrangian Method (ALM) [29,26,3], which is described in Algorithm 1.

Algorithm 1: Augmented Lagrangian Method (ALM)
1 Input: choose initial guess Λ 0 for the dual variables, and set k := 0; 2 while certain stopping criterion is not reached do 3 Minimize the augmented Lagrangian function with respect to the prime variables X: where the augmented Lagrangian function fo problem (1.1) is defined as Update the Lagrangian multipliers Update the penalty parameter β if necessary. Set k := k + 1.
It is well-known that the augmented Lagrangian function (2.5) is an exact penalty if the parameter β is sufficiently large. Algorithm 1 works very well for problem with linear constraints. For optimization problems with nonlinear constraints, it is not clear how to choose the parameter β in practice, which is very sensitive to the numerical performance.
The purpose of this work is to find an infeasible algorithm for solving (1.1) at similar cost of the existent feasible methods. Otherwise, we can hardly gain much from the parallelization. To this end, we carefully test Algorithm 1 and try our best to tune the parameter β. Unfortunately, for solving optimization problems with orthogonality constraints (1.1), the efficiency of classical ALM is far from being satisfactory. Therefore, we need to employ a new idea to remould the classical ALM. According to the conditions (2.1) and (2.2), it is not difficult to verify that the Lagrangian multipliers Λ have the following closed-form expression at any first-order stationary point, A straightforward idea is to use the following symmetrized form of (2.7) as a new multipliers updating rule. The symmetrization is necessary because the symmetry of the expression ∇f (X) ⊤ X can not be guaranteed in each iteration.
As we will demonstrate in the following lemma and the theoretical analyses in Section 4, an explicit lower bound of the penalty parameter β can be estimated if updating rule (2.8) is applied. Hence, the update of the penalty parameter β can be waived. Moreover, the numerical experiments verify the validation of this new updating rule.
Proof. For convenience, we abuse the notation slightly by deleting the superscript * from X * . First, we have Since X is the second-order stationary point of (2.9) with Λ = ∇f (X) ⊤ X, we have Left multiplying X ⊤ into both sides of (2.14), we have Suppose X = U ΣV ⊤ is the singular value decomposition of X in economy-size, which implies X ⊤ X = V Σ 2 V ⊤ . Then, we further have Left multiplying V ⊤ and right multiplying V to both sides of the above equality, we arrive at Taking the Φ operator and using the fact that On the other hand, since n ≥ 2p, there existsŨ ∈ S n,p satisfyingŨ ⊤ U = 0. let S =Ũ DV ⊤ . If S = 0, we substitute S into (2.11) and obtain Here I stands for the identity mapping from R n×p to R n×p . Combining with the second-order optimality condition (2.13), relationship (2.18) and the assumption on β, we have which leads to contradiction. Hence, S = 0, which immediately implies that Σ = I p . Therefore, we have X ∈ S n,p . Together with (2.12) and (2.13), we can easily show that the optimality condition (2.1) and (2.2) hold. This completes the proof. Lemma 2.3 guarantees that the augmented Lagrangian function is still an exact penalty function with the Lagrangian multipliers updated by explicit formula (2.8). However, to achieve the convergence results for first-order methods, we need a first-order version of Lemma 2.3. Moreover, to obtain the global convergence rate, the feasibility should be controlled by the first-order optimality violation. LEMMA 2.4. For any X * satisfying σ min (X * ) > 0, suppose β > (||∇f (X * )|| 2 · ||X * || 2 + δ) /σ 2 min (X * ) with δ > 0. Then it holds with Λ * = Ψ (∇f (X * ) ⊤ X * ). In particular, if it happens that X * is a first-order stationary point of , then X * is also a first-order stationary point of problem (1.1).
Proof. For brevity, we denote G = ∇ X L β (X, Λ). Left multiplying X ⊤ into both sides of (2.10) and using the singular value decomposition X = U ΣV ⊤ , we have Left multiplying V ⊤ and right multiplying V to both sides of the above equality, we obtain Taking the Φ operator and using the fact (2.16), we arrive at Hence, it holds that Submitting (2.21) into (2.20), we arrive at and complete the proof.
3. Parallelizable Algorithms. In this section, we introduce a parallelizable approach and one of its variant for optimization problem with orthogonality constraints (1.1). Both of these two approaches are based on the augmented Lagrangian function (2.5) and employ the new idea of updating the multipliers by explicit expression instead of dual ascent step in Algorithm 1.
Another distinction between our algorithms and the classical ALM is that the minimization subproblem for the prime variables is replaced by a proximal linearized approximation.
3.1. The Proximal Linearized Augmented Lagrangian Algorithm. We describe our main algorithm framework in Algorithm 2.
The main calculation costs of Algorithm 2 concentrate at Step 3 and 4.
Step 3 only involves BLAS3 calculation. 7 Algorithm 2: Proximal Linearized Augmented Lagrangian Algorithm (PLAM) 1 Input: choose initial guess X 0 , and set k := 0; 2 while certain stopping criterion is not reached do 3 Compute the Lagrangian multipliers Minimize the following proximal linearized Lagrangian function Set k := k + 1.
The minimization subproblem (3.2) in Step 4 is nothing but a gradient step where the last step is due to the updating formula (3.1). Apparently, the arithmetic operations involved in (3.3) belong to BLAS3 as well.
We notice that 1/η k is nothing but the stepsize of gradient step. Hence, the proximal parameter η k can be chosen in the same manner as how we choose stepsize for gradient methods. This issue will be described in details in Section 5.

Parallelizable Column-wise Block
Minimization. An obvious demerit of PLAM is the boundedness of the iterate sequence can hardly be expected without any restriction on the penalty parameter β and the proximal parameter η k . Theoretically, to guarantee the global convergence, β should be sufficiently large. Accordingly, η k should be large as well which means sufficiently small stepsize is required and slow convergence can be expected. In fact, according to the empirical observations, the performance of PLAM is very sensitive to parameters β and η k . In other word, it is not easy to tune these two parameters to guarantee good performance of Algorithm 2 in general. Therefore, we put forward an upgraded version of PLAM. It is based on PLAM, but redundant column-wise unit sphere constraints are imposed to Step 4. Therefore, the proximal gradient takes the place of the gradient step in the Step 4 of Algorithm 2. With redundant constraints, the resulting iterate sequence will then be restricted to a compact set and hence bounded. We describe the framework of this upgraded PLAM in Algorithm 3.
Subproblem (3.5) in Algorithm 3 can be solved in a column-wisely parallel fashion. In fact, it is of closed-form solution For PCAL, we can update the Lagrangian multipliers in the same manner as PLAM, i.e. by formula (3.1). To 8 Algorithm 3: Parallelizable Column-wise Block Minimization for PLAM (PCAL) 1 Input: choose initial guess X 0 , and set k := 0; 2 while certain stopping criterion is not reached do 3 Compute the Lagrangian multipliers by (3.1) or , and set k := k + 1.
obtain a better performance, we can also use the heuristic formula (3.4). The motivation of updating formula (3.4) comes from the following observation. In the KKT condition (2.1), we impose an additional term for the redundant sphere constraints. Namely, where D is a diagonal matrix. Furthermore, D is determined by the Lagrangian multiplier of X i in the subproblem (3.5).

Computational Cost.
In this subsection, we compare the computational cost per iteration among MOp-tQR, PLAM and PCAL. The computational cost of the basic linear algebra operations and the overall costs of the aforementioned algorithms are listed in Table 3.1.
Here, those terms in red represent the corresponding operations that cannot be parallelized.
In practice, we calculate XΨ (∇f (X) ⊤ X) instead of X(∇f (X) ⊤ X) for KKT evaluation since they are very close to each other around any first-order stationary point. Consequently, it saves 2np 2 flops computational cost.

Convergence of PLAM.
In this section, we focus on the theoretical analyses of our proposed PLAM. The global convergence, worst-case complexity and Q-linear local convergence rate will be established under different mild assumptions.

Global Convergence of PLAM.
Besides blanket Assumption 1.1, to prove the convergence of Algorithm 2, we need to impose a mild condition on the initial guess, and restrictive conditions on β and η k . To facilitate the narrative, we first state all these conditions here. ASSUMPTION 4.1. For a given X 0 , we say it is a qualified initial guess, if there exists σ ∈ (0, 1) so that Assumption 4.1 is not restrictive at all. Therefore two types of points satisfying this assumption and can be 9 evaluate function Now, we list all the special notations to be used in this section.
We introduce the following merit function According to the twice continuous differentiability of f (X), ∇f (X) is Lipschitz continuous on the compact set C. Namely, there exists constant L h > 0, related to β, so that The algorithm parameters β and η k , and the constants used in the proof can be selected by the following rules.
Now we give a sketch of our proof. Suppose {X k } is the iterate sequence generated by Algorithm 2. The main steps of the proof include: (1) Any iterate X k is in C, and σ is a unified lower bound of the smallest singular values of the iterates X k ; (2) The merit function h(X) is bounded below; (3) {h(X k )} monotonically decreases, and hence {X k } has at least one convergent subsequence; (4) Any cluster point of {X k }, say X * , is a first-order stationary point of the augmented Lagrangian function (5) Any cluster point of {X k }, say X * , is a first-order stationary point of the original optimization problem with orthogonality constraints (1.1).
Next we provide five concrete lemmas or corollaries following the above-mentioned sketch.
Proof. We use mathematical induction. The argument (4.6) directly holds for X 0 resulting from Assumption 4.1. Next we investigate whether (4.6) holds at X k+1 provided that it holds for X k .

It is not difficult to verify that
holds for any X k ∈ C. By using the facts X k ∈ C, (4.1) and (4.5), we have which implies ||X k+1 ⊤ X k+1 − I p || F ≤ R. This shows (4.6) is true for k + 1.
According to the facts σ min (X k ) ≥ σ and X k ∈ C, we have By using the fact that tr(AB) = tr(AB ⊤ ) if A is symmetric, we have Hence, we have Notice that L c = 2R + 4M 2 is the Lipschitz constant of ∇c(X) over C. Due to the fact (4.5), (4.8) and (4.9), we have According to the Taylor expansion, we have According to assumption, R ≤ 1 − σ 2 , we can easily obtain that σ min (X k+1 ) ≥ σ. This completes the proof.
This lemma immediately holds by the continuous differentiability of h(X) and the compactness of C, and hence, the proof is omitted.
We keep using the notations (4.7) and investigate where the second last inequality is implied by relation (4.9). Hence, we arrive at Substituting (4.5) and (4.11) into the Taylor expansion, we have We complete the proof.
With the boundedness of h(X) at C, Lemma 4.5 immediately implies the convergent of {h(X k )}. More precisely, we have the following corollary.
Moreover, {X k } has at least one convergent subsequence. Any cluster point of {X k }, X * , is a first-order stationary point of minimizing the augmented Lagrangian function Proof. This is a direct corollary of Lemmas 4.3 and 4.5.
Finally, we give the global convergence rate of PLAM, namely, the worst case complexity.
Proof. This is a direct corollary of Lemma 2.4 and Theorem 4.7.
REMARK 4.9. The sublinear convergence rate of Corollary 4.8 actually tells us that Algorithm 2 terminates after

Local Convergence Rate of PLAM and PCAL.
In this subsection, we consider the local convergence of PLAM once the optimization problem with orthogonality constraints (1.1) has an isolated local minimizer. THEOREM 4.10. Suppose X * is an isolated minimizer of (1.1), and we denote The algorithm parameters satisfy β ≥ L+MN +τ 2 and η k ∈ [η,η], whereη ≥ η ≥ L + M N + 2β. Then, there exists ε > 0 such that starting from any X 0 satisfying ||X 0 − X * || F < ε, and the iterate sequence {X k } generated by Algorithm 2 converges to X * Q-linearly.
Subtracting the second one from the first one and using the Taylor expansion, we have 16) where δ k = X k − X * . Recall the expression of Hessian (2.13), the fact that ∇f (X * ) ⊤ X * = Ψ (∇f (X * ) ⊤ X * ) and the assumption on η, we have 14 On the other hand, δ k can be decomposed as the summation of three terms: δ k = X * S + X * W + K, (4.18) where S ∈ R p×p is symmetric, W ∈ R p×p is skew-symmetric, K ∈ R n×p is perpendicular to X * . Since X * is a strict local minimizer, and T (X) is closed, we have τ > 0. Hence, it holds that as X * W + K ∈ T (X). Moreover, it follows from the assumption on β that Combining (4.19), (4.20), the symmetry of S, the skew symmetry of W , K ⊤ X * = 0 together with the assumption on η, we arrive at Notice that (4.17) implies the positive semi-definiteness of the linear operator Together with (4.21), we can conclude that which completes the proof. 5. Numerical Experiments. In this section, we evaluate the numerical performance of our proposed algorithms PLAM and PCAL. We first introduce the implementation details and the testing problems in Subsection 5.1 and 5.2, respectively. Then, we report the numerical experiments which are mainly of three folds.
In the first part, we mainly determine the default settings of our proposed algorithms, which will be discussed in Subsection 5.3. Then, in Subsection 5.4, we compare our PLAM and PCAL with a few existing solvers by testing a bunch of instances, which are chosen from a MATLAB toolbox KSSOLV [36]. Another one is the proximal parameter η, whose reciprocal is the step size of the gradient step in Algorithm 2 and 3. Similar to β, we can not use the rigorous restriction in the theoretical analysis. In practice, we have the following strategies to choose this parameter: (iii) Barzilai-Borwein(BB) strategy [2]: (iv) Alternating BB strategy [9]: Unless specifically mentioned, the stopping criterion used for both serial and parallel experiments can be described as follows, The maximum number of iteration for all those solvers is set to 3000.

Testing Problems.
In this subsection, we introduce six types of problems which will be used in the numerical experiments. 4 More information at http://lsec.cc.ac.cn/chinese/lsec/LSSC-IVintroduction.pdf Problem 1: A simplification of discretized Kohn-Sham total energy minimization.
where the matrix L ∈ S n and ρ(X) := diag(XX ⊤ ). In the numerical experiments, we set α = 1, and L is randomly generated by Gauss distribution, i.e., L=randn(n) in MATLAB language, and set L := 1 2 (L + L ⊤ ). In this instance, s = L 2 .
Problem 2: A class of quadratic minimization with orthogonality constraints.
where the matrices A ∈ S n and G ∈ R n×p . This problem is adequately discussed in [11]. In the numerical experiments, the matrices A and G are randomly generated in the same manner as in [11]. Namely, where the matrices P =qr(rand(n,n))∈ R n×n ,Q=rand(n,p)∈ R n×p , Q ∈ R n×p and Q i =Q i /||Q i || 2 (i = 1, 2, ..., p), and matrices Λ ∈ R n×n and D ∈ R p×p are diagonal matrices with Here, parameter θ ≥ 1 determines the decay of eigenvalues of A; parameter ζ ≥ 1 refers to the growth rate of column's norm of G. The parameter κ > 0 represents the scale difference between the quadratic term and the linear term. The default setting of these parameters are κ = 1, θ = 1.01, ζ = 1.01, ξ = 1. In this instance, s = A 2 .
Problem 3: Rayleigh-Ritz trace minimization, which is a special case of Problem 2.
where the matrices A ∈ S n . In our experiments, the matrix A is generated in the same manner as in Problem 2. In this instance, s = A 2 .
Problem 4: Another class of quadratic minimization with orthogonality constraints. X ⊤ X = I p . (5.8) where the matrices A ∈ S n and B ∈ S p . This problem is out of the scope of problems discussed in [11], but can be solved by PLAM or PCAL. The matrices A and B are randomly generated by A =randn(n), A := 1 2 (A + A ⊤ ) and B =randn(p), B := 1 2 (B + B ⊤ ), respectively. In this instance, s = A 2 · B 2 Problem 5: Discretized Kohn-Sham total energy minimization instances from KSSOLV [36].
where the discretized Kohn-Sham total energy function E(X) is defined by (1.2). All the data comes from MATLAB toolbox KSSLOV.
Problem 6: A synthetic instance of discretized Kohn-Sham total energy minimization.
where the matrix L ∈ R n×n and ρ(X) := diag(XX ⊤ ). The parameter γ = 2( 3 π ) 1/3 and ρ(X) 1 3 denotes the component-wise cubic root of the vector ρ(X). This problem adopts a special exchange functional − 3 4 γρ(X) ⊤ ρ(X) 1 3 (the correlation term is ignored), which is introduced in [22]. The generation of L is in the same manner as in Problem 1.

Default Settings.
In this subsection, we determine the default settings for the proposed algorithms PLAM and PCAL.
In the first experiment, we test PLAM and PCAL with these four different choices of η k on Problem 1-4. Here, we only illustrate the results of η BB1 for strategy (iii), since its performances overwhelms those with η BB2 . The penalty parameter is fixed as β = s + 0.1. FIG. 5.1 shows the results of PLAM and PCAL with different η k . From subfigures (a)-(d), we observe that PLAM with η ABB outperforms others. Under the same setting, a comparison among PCAL with different η k is reported in subfigures (e)-(h). We notice that PCAL with η ABB is superior to the other η k choices. Then we set η ABB as the default setting for PLAM and PCAL. We next compare the performance among PLAM and PCAL variations corresponding to different β. In the comparison, we set β varying among 0, 0.01s, 0.1s, s + 0.1, 10s + 1. The proximal parameter is fixed as its default η = η ABB . We present all the numerical results in Figure 5.2. We notice from subfigures (a)-(d) that PLAM with small β might be divergent in some cases, while large β causes slow convergence. Therefore, a suitable chosen β, often unreachable in practice, is necessary for good performance of PLAM. On the other hand, the dependence on β of PCAL can be learnt from subfigures (e)-(h). The smaller β for PCAL has the better performance in some instances, and the behavior of PCAL is completely not sensitive to β in other instances. To take more distinctive look at the difference between PLAM and PCAL, we present a comparison in Figure 5.3. Therefore, in practice, we suggest an approximation of s to be the default β of PLAM and 1 for PCAL. Since it is easier to tune β for PCAL than PLAM, we choose PCAL to be the default algorithm of ours in Subsection 5.5. There are two distinctions between PLAM and ALM. Firstly, a gradient step takes the place of solving the subproblem to some given precision in the update of the prime variables. Secondly, a closed-form expression is used to update the Lagrangian multipliers in stead of dual ascend. In order to show that the new update formula for multipliers is a crucial fact of the efficiency of PLAM and PCAL, we compare PLAM and PCAL with PLAM-DA and PCAL-DA, respectively. Here PLAM-DA and PCAL-DA stand for Algorithm 2 and 3 with Step 3 using dual ascend to 19 update the multipliers, respectively. We report the numerical results in Figure 5.4. It can be observed that the closedform expression for updating Lagrangian multipliers is superior to dual ascend in solving optimization problems with orthogonality constraints.

. A comparion bewteen PLAM and PCAL on multilplier
In the end of this subsection, we show how KKT and feasibility violations decay in the iterations, when PLAM and PCAL are used to solve Problem 1. The numerical results are presented in Figure 5.5. We notice that the decay of feasibility violations is nonmonotone and has a similar variation tendency as KKT violations, which coincides our theoretical analysis Lemma 2.4. If we want a high accuracy for the feasibility but a mild one for KKT conditions, we can set a mild tolerance for KKT violation and impose an orthonormalization step as a post process when we obtain the last iterate by PLAM or PCAL. Table 5.1 illustrates that such post process does not affect the KKT violation, but do improve the feasibility. Here, "stop" and "orth" represent the relative values at the last iterate and the one after post process, respectively. Hereinafter, the orthonormalization post process, achieved by an internal function qr(·) assembled in MATLAB, is the default last step of PLAM and PCAL. Our test is based on KSSOLV 5 [36], which is a MATLAB toolbox for electronic structure calculation. It allows 5   researchers to investigate their own algorithms easily and friendly for different steps in electronic structure calculation. We choose two integrated solvers in KSSOLV. One is the self-consistent field (SCF) iteration, which minimize a quadratic surrogate of the objective of (5.9) with orthogonality constraints in each iteration [21]. SCF and its variations are the most widely used in real KSDFT calculation. The other one is called trust-region direct constrained minimization (TRDCM) [38], which combines the trust-region framework and SCF to solve the subproblem. Besides SCF and TRDCM, which are particularly for KSDFT, we also pick up two state-of-the-art solvers in solving general optimization problems with orthogonality constraints. One is OptM 6 , which is based on the algorithm proposed in [34]. OptM adopts Cayley transform to preserve the feasibility on the Stiefel manifold in each iteration. Nonmonotone line search with BB stepsize is the default setting in OptM. Another existing solver for comparison, we intend to choose MOptQR, which is based on a projection-like retraction method introduced in [1]. Its original version is MOptQR-LS (manifold QR method with line search 7 ). For fair comparison, we implement the same alternating BB stepsize strategy as PLAM and PCAL to MOptQR-LS, and form the MOptQR used in this section.
We select 18 testing problems with respect to different molecules, which are assembled in KSSOLV. For all the methods, the stopping criterion is set as (I n − XX ⊤ )H(X)X F < 10 −5 . And we set the max iteration number MaxIter = 200 for methods SCF and TRDCM, while MOptQR, OptM, PLAM and PCAL set their max iteration number with MaxIter = 1000 to get a comparable solution with other methods. The penalty parameter β PLAM for PLAM is tuned case by case to achieve a good performance. Meanwhile, β PCAL for PCAL is always set as 1. Other parameters for all these methods take their default values. For all of the testing algorithms, we set the same initial guess X 0 by using the function "getX0", which is provided by KSSOLV. The numerical results are illustrated in Tables 5.2, 5.3 and 5.4.
Here, "E tot " represents the total energy function value, and "KKT violation", "Iteration", "Feasibility violation" and "Time(s)" stand for (I n − XX ⊤ )H(X)X F , the number of iteration, X ⊤ X − I p F and the total running wall-clock time in second, respectively. From the tables, we observe that PCAL has a better performance than other algorithms, and in most cases, it obtains a comparable total energy function value and a lower KKT violation. In particular, in the large size problem "graphene30", PCAL achieves the same total energy function value and same magnitude KKT violation in much less CPU time than others. In the problem "qdot", we observe that only PLAM and PCAL can output a point satisfying the KKT violation tolerance, while all the other algorithms terminate abnormally. Therefore, we can conclude that PCAL and PLAM perform comparable with the existent feasible algorithms in solving discretized Kohn-Sham total energy minimization.

Parallel Efficiency.
In this subsection, we examine the parallel efficiency of our algorithms PLAM and PCAL. To investigate the parallel scalability, we need to test large scale problems in a single core, which consumes Both algorithms are implemented in the C++ language and parallelized by OpenMP. The linear algebra library we used in comparison is Eigen 8 (version 3.3.4), which is an open and popular C++ template library for matrix computation. We define the speedup factor for running a code on m cores as speedup-factor(m) = wall-clock time for a single core run wall-clock time for a m-core run .
BLAS3 type arithmetic operations contribute a high proportion in computational cost in both PCAL and MOptQR.
Therefore, a good parallel strategy for BLAS3 calculation is unnegligible in saving CPU time. Given this, we first determine the parallel strategy for matrix-matrix multiplication by a set of tests. We have two choices. The library Eigen provides its own multi-threading computing 9 that is the default parallel strategy for dense matrix-matrix products and row-major-sparse * dense vector/matrix products in OpenMP. Another strategy is to parallelize BLAS3 computation in the manner of column-wise product. Namely, when we calculate AB, we multiply matrix A by each column of B in parallel. To figure out which strategy is better, we test the parallel scalability of BLAS3 computation under these two schemes. We generate A=Random(1000,10000) and B=Random(10000,1000), where "Random(·,·)" is an internal generation function provided by Eigen. We run the code in parallel with 1, 2, 4, 8, 16, 32, 64 and 96 cores, respectively. The result of matrix-matrix multiplication AB is illustrated in Figure 5.6. "Eigen" and "Column-wise" represent the default parallel strategy and column-wise product strategy, respectively. We can observe that columnwise parallelization obviously outperforms the default setting of Eigen in multi-threading computing. Hence, in the following implementation, we choose column-wise parallelization strategy for BLAS3 in our experiments. Next, we investigate the parallel scalability of the new proposed PCAL and MOptQR. According to the existent numerical report of Eigen 10 , we select the class "LLT" in Eigen to compute QR factorization. The calculation of orthonormalization consists of a small size (p-by-p) Cholesky decomposition and solving a p-by-p linear system. The maximum number of iterations for MOptQR and PCAL is set to 1000. All the parameters for MOptQR and PCAL take their default values. The initial guess X 0 is generated by X 0 =random(n,p)" and X 0 =qr(X 0 ).
We first focus on the test Problems 1 and 2. For Problem 1, we set L as a block diagonal matrix, i.e., L = Diag(L 1 , . . . , L s ), where L i ∈ R 5×5 is a tridiagonal matrix with 2 on its main diagonal and −1 on subdigonal, for i = 1, . . . , s. The coefficient α is set to 1. For the generation of Problem 2, we set A as a tridiagonal matrix with 2 on its main diagonal and −1 on subdigonal, and G=Random(n,p). The advantage of such generation is to make function value and gradient calculations parallelizable. In the first group of tests, we aim to figure out how MOptQR and PCAL perform with the increasing width of variables. We set n = 10000 and p varying from a set of increasing values {500, 1000, 1500, 2000, 2500}. Both algorithms are run in parallel with 96 cores. The wall-clock time results are shown in Figure 5.7. Here, "#cores" stands for the number of cores. From Figure 5.7, we notice that PCAL always takes less amount of wall-clock time than MOptQR. As the width of the matrix variable increases, the running time of MOptQR increases much more rapidly than PCAL. In Figure 5.8, we show wall-clock time of four categories: three categories: "BLAS3" (dense-dense matrix multiplication), "Func" (function value and gradient evaluation) and "Orth" (orthonormalization including QR factorization for MOptQR and the final correction step in PCAL). These are the major computational components of both PCAL and MOptQR, albeit in different proportions. We have to clarify two issues: firstly, we categorize these categories of calculation only at the highest solver level. As such, any matrix-matrix multiplication involved in function value and gradient evaluation is not counted as in the "BLAS3" category. Secondly, although the correctness of such a classification scheme may be debatable, it does not alter the overall fact, as is clearly shown by our computational results, that the category "BLAS3" is much more scalable than the category "Orth" on our test platform. The running time of each category is measured in terms of the percentage of wall-clock time spent in that category over the total wall-clock time. We can clearly see that for PCAL the run time of "BLAS3" dominates the entire computation in almost all cases. The "BLAS3" time increases steadily as p increases from 500 to 2500, while the "Func" time decreases steadily. The run time of "Orth" is negligible. However, for MOptQR, the "BLAS3" time takes around 60% of total run time and decreases steadily with the increasing of p. Meanwhile, the "Orth" time takes around "40%" of total run time and increasing steadily. Now, we set n = 10000 and p = 1000, 2000, and run PCAL and MOptQR in parallel with 1, 2, 4, 8, 16, 32, 64 and 96 cores, respectively. Figure 5.9 and 5.10 illustrate the speedup factors associated with total running wall-clock time, "BLAS3", "Func" and "Orth", respectively. From these two figures, we can observe that BLAS3 operation has high parallel scalability, while the speedup factor of "Orth" increases slowly as the number of cores increases, which directly leads to the higher overall scalability of PCAL than MOptQR. Moreover, as the width of the matrix variable increasing, the advantage of PCAL in parallel scalability becomes more obvious.
In the end, we test Problem 6 under n = 10000, p = 1000. Figure 5.11 illustrate the results of speedup factors associated with total running wall-clock time, "BLAS3", "Func" and "Orth" of PCAL and MOptQR, respectively. We can learn from this figure that the overall scalability of PCAL is again superior to that of MOptQR. 6. Conclusion. Optimization problems with orthogonality constraints have wide application in materials science, machine learning, image processing and so on. Particularly, when we apply Kohn-Sham density functional theory (KSDFT) to electronic structure calculation, the last step is to solve a Kohn-Sham total energy minimization with orthogonality constraints. There are plenty of existing algorithms based on manifold optimization, which work quite well when the number of columns of the matrix variable p is relatively few. With the increasing of p, a bottleneck of existent algorithms emerges, that is, lack of concurrency. The main reason leads to this bottleneck is that the orthonormalization process has low parallel scalability.
To solve this issue, we need to employ infeasible approaches. However, previous infeasible approaches including augmented Lagrangian method (ALM) is far less efficient than retraction based feasible methods. Even though the parallelization reduces the running time of ALM more significantly than that of manifold methods, ALM is still less efficient than manifold methods in parallel computing. The main purpose of this paper is to provide practical efficient infeasible algorithms for optimization problems with orthogonality constraints. Our main motivation is that the Lagrangian multipliers have closed-form expression at any stationary points. Hence, we use such expression to update multipliers instead of dual ascent step, at the same time, the subproblem for the prime variables only takes one gradient step instead of being solved to a given tolerance. The resultant algorithm, called PLAM, does not involve any orthonormalization. PLAM is comparable with the existent feasible algorithms under well chosen penalty parameter β. To avoid such restriction, we propose a modified version, PCAL, of PLAM. The motivation of PCAL is to use normalized gradient step instead of gradient step in updating prime variables. The numerical experiments show that PCAL works efficient, robust and insensitive with penalty parameter β. Remarkably, it outperforms the existent feasible algorithms in solving the KSDFT problems in MATLAB platform KSSOLV. We also run PCAL 28 and MOptQR, an excellent representative of retraction based optimization approach, in parallel with up to 96 cores.
Numerical experiments illustrate PCAL has higher scalability than MOptQR, and its superiority becomes more and more noticeable with the increasing of p.
The potential of PCAL has already emerged. In the future work, we will apply our PCAL to real KSDFT calculation.