W W Hager
email: hager@math.u.edu

G Bencteux

E Cancès

C Le Bris

Analysis of a Quadratic Programming Decomposition Algorithm *

Keywords: quadratic programming, orthogonality constraints, domain decomposition method, electronic structure calculations

We analyze a decomposition algorithm for minimizing a quadratic objective function, separable in x 1 and x 2 , subject to the constraint that x 1 and x 2 are orthogonal vectors on the unit sphere. Our algorithm consists of a local step where we minimize the objective function in either variable separately, while enforcing the constraints, followed by a global step where we minimize over a subspace generated by solutions to the local subproblems. We establish a local convergence result when the global minimizers nondegenerate. Our analysis employs necessary and sucient conditions and continuity properties for a global optimum of a quadratic objective function subject to a sphere constraint and a linear constraint. The analysis is connected with a new domain decomposition algorithm for electronic structure calculations.

Analyse d'un algorithme de décomposition pour un problème de programmation quadratique Résumé : On présente l'analyse numérique d'un algorithme de décomposition pour la minimisation d'une fonction coût quadratique, séparable en x 1 et x 2 , sous la contrainte que x 1 et x 2 sont orthogonaux sur la sphère unité. Notre algorithme consiste en une étape locale où la fonction coût est minimisée séparément en chacune de ses deux variables, en respectant les contraintes. Cette première étape est suivie d'une étape globale où on minimise la fonction coût sur un sous-espace généré par les solutions de l'étape locale. Un théorème de convergence locale est établi quand les minimiseurs globaux ne sont pas dégénérés. Notre analyse utilise les conditions nécessaires et susantes et les propriétés de continuité du minimum global d'une fonction coût quadratique minimisée sous une double contrainte sphérique et linéaire. Cette analyse est reliée à un nouvel algorithme de décomposition de domaine pour les calculs de structure électronique.

Mots-clés : programmation quadratique, contraintes d'orthogonalité, méthode de décomposition de domaine, calculs de structure électronique [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF] Introduction

In [START_REF] Barrault | Multilevel domain decomposition for electronic structure calculations[END_REF] and [START_REF] Bencteux | Domain decomposition and electronic structure computations: a promising approach[END_REF] we develop a multilevel domain decomposition algorithm for electronic structure calculations which has been extremely eective in computing electronic structure for large, linear polymer chains. Both the computational cost and memory requirement scale linearly with the number of atoms. Although this algorithm has been very eective in practice, a theory establishing convergence has not yet been developed. The algorithm in [START_REF] Barrault | Multilevel domain decomposition for electronic structure calculations[END_REF][START_REF] Bencteux | Domain decomposition and electronic structure computations: a promising approach[END_REF] was motivated by a related decomposition algorithm for a quadratic programming problem with an orthogonality constraint. In this paper, we develop a convergence theory for the decomposition algorithm.

Let H 1 and H 2 be symmetric n by n matrices. We consider the following quadratic optimization problem:

min F (x 1 , x 2 ) := x T 1 H 1 x 1 + x T 2 H 2 x 2 (1)
subject to x T 1 x 1 = 1 = x T 2 x 2 , x T 1 x 2 = 0.

In other words, nd orthogonal unit vectors x 1 and x 2 which minimize the separable quadratic objective function. Our algorithm for (1) consists of a local step where we minimize F over each variable separately, while enforcing the constraints, followed by a global step where we optimize over a subspace generated by the iterates of the local step.

There are two modes of the local step, a forward and a reverse mode. In consecutive iterations, we employ the forward mode followed by the reverse mode. If x k = (x k1 , x k2 ) is the iterate at step k, then the forward and the reverse modes of the local step are the following: forward y k1 ∈ arg min {F (z, x k2 ) : z = 1, z T x k2 = 0}, y k2 ∈ arg min {F (y k1 , z) : z = 1, z T y k1 = 0}, reverse y k2 ∈ arg min {F (x k1 , z) : z = 1, z T x k1 = 0}, y k1 ∈ arg min {F (z, y k2 ) : z = 1, z T y k2 = 0}.

(

) 2
Here and throughout the paper,

• denotes the Euclidean norm.

The problem which must be solved in electronic structure calculations is more general than [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF] and the multilevel algorithm developed in [START_REF] Barrault | Multilevel domain decomposition for electronic structure calculations[END_REF][START_REF] Bencteux | Domain decomposition and electronic structure computations: a promising approach[END_REF] is more complex than (2). For example, in electronic structure calculations, H 1 and H 2 could be of dierent dimensions and the orthogonality condition x T 1 x 2 = 0 in (1) would be replaced by the more general condition x T 1 Px 2 = 0 where P is rectangular. Nonetheless, the algorithm studied in this paper was the basis for the more general algorithm developed in [START_REF] Barrault | Multilevel domain decomposition for electronic structure calculations[END_REF][START_REF] Bencteux | Domain decomposition and electronic structure computations: a promising approach[END_REF], and our analysis is an initial step towards justifying and understanding the convergence properties of the more general algorithm.

One can think of either the forward or reverse modes as a block Gauss-Seidel iteration [5, p. 323]. In the forward mode, we rst hold the second block of variables x k2 xed and we optimize over the rst block of variables to obtain y k1 ; in the second step, we hold the rst block of variables xed at y k1 and we optimize over the second block to obtain y k2 .

RR n°6288

In general, the local steps converge to a limit which may not be a stationary point of (1). To achieve convergence to a stationary point for (1), each local step, either forward or reverse, is followed by a global step where we minimize F over the subspace spanned by the following 4 vectors, while enforcing the constraints of (1):

y k1 0 , y k2 0 , 0 y k1 , 0 y k2 (3)
After imposing the two normalization condition x T 1 x 1 = 1 and x T 2 x 2 = 1 and the orthogonality condition x T 1 x 2 = 0 on the subspace, we are left with a one dimensional curve of feasible points in the subspace spanned by the 4 vectors (3). This curve can be expressed in the following way:

z k (s) = 1 √ 1 + s 2 (y k + sd), (4)
where

d = ± y k2 -y k1
and y k = y k1 y k2 .

The vector z k (s) lies in the space spanned by the 4 vectors in (3) for each choice of s; the orthogonality condition holds since

(1 + s 2 )z T k1 z k2 = (y k1 ± sy k2 ) T (y k2 ∓ sy k1 ) T = ±s ∓ s = 0;
and the factor 1/ √ 1 + s 2 ensures that the two components of z k are unit vectors. Let F k (s) = F (z k (s)) be the objective function evaluated along the search direction. For convenience, the sign in the denition of d in the global step is chosen so that F k (0) ≤ 0. At iteration k in the global step, we set

x k+1 = z(s k ), (5)
where s k is the stepsize.

The motivation for optimizing over the subspace spanned by the 4 vectors (3) is the following: First, the subspace should include the original vectors (y k1 , 0) and (0, y k2 ) to ensure that the objective function value decreases. In order to further broaden the search space, we should consider vectors orthogonal to the original vectors. Since the vectors (y k2 , 0) and (0, y k1 ) are orthogonal to the original vectors, they are suitable for inclusion in the subspace. Finally, as we will see in Section 4, the optimality condition associated with the global step and with these 4 vectors provides a link between the subproblems which is exploited to obtain convergence.

Notice that when H 1 and H 2 are 2 by 2 matrices, the 4 vectors in (3) span R 4 . Hence, in the 2 by 2 case, the global step yields a global optimum for [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF]. More generally, we nd that the local steps steer the iterates into a subspace associated with the eigenvectors of H 1

INRIA

for k = 0, 1, 2, . . .

If k is even, perform a forward step:

y k1 ∈ arg min {F (z, x k2 ) : z = 1, z T x k2 = 0}, y k2 ∈ arg min {F (y k1 , z) : z = 1, z T y k1 = 0},
Else perform a reverse step:

y k2 ∈ arg min {F (x k1 , z) : z = 1, z T x k1 = 0}, y k1 ∈ arg min {F (z, y k2 ) : z = 1, z T y k2 = 0}. Global step: Set x k+1 = z(s k ) where z k (s) = 1 √ 1+s 2 (y k + sd), d = ± y k2 -y k1
, and

s k ∈ arg min {F k (s) : s ∈ [0, -ρF k (0)]}, F k (s) = F (z k (s)).
The sign of d is chosen so that F k (0) ≤ 0. Ideally, the stepsize s k is the global minimum of F k (s) over all s. However, the convergence analysis for this optimal step is not easy since y k -x k+1 could be on the order of 1 for all k. For example, if H 1 = H 2 , then F k (s) is constant, independent of s; consequently, any choice of s is optimal, and there is no control over the iteration change. To ensure global convergence of the algorithm, we restrict the stepsize to an interval [0, -ρF k (0)], where ρ is a xed positive scalar. In other words, we take

s k ∈ arg min {F k (s) : s ∈ [0, -ρF k (0)]}. (6)
Notice that s k = 0 when F k (0) = 0 and the global step is skipped. In practice, we observe convergence when s k is a global minimizer of F k . The constraint on the stepsize is needed to rigorously prove convergence of the iteration. For reference, the complete algorithm is recapped in Figure 1.

As we show later in (41), F k (0) tends to zero. Hence, the constraint s ∈ [0, -ρF k (0)] on the stepsize in the line search [START_REF] Nocedal | Numerical Optimization[END_REF] implies that the iteration dierence x k+1 -y k tends to zero. Another approach for controlling the stepsize is to employ a trust region scheme [4, p. 129] where we minimize F in the subpace (3) and a ball of radius ρ k centered at y k . If ρ k tends to zero, then the change x k+1 -y k again tends to zero. The update (6) amounts to a trust region step with a special choice for the trust region radius.

Since F is a pure quadratic, the objective function satises

F (x 1 , x 2 ) = F (-x 1 , x 2 ) = F (x 1 , -x 2 ) = F (-x 1 , -x 2 ).

RR n°6288

Hence, if y kj is a minimum in a subproblem at iteration k, then so is -y kj . In order to carry out the analysis, it is convenient to choose the signs so that following inequalities hold:

x T k2 H 1 y k1 ≥ 0 and y T k1 H 2 y k2 ≥ 0 (forward mode) x T k1 H 2 y k2 ≥ 0 and y T k2 H 1 y k1 ≥ 0 (reverse mode) [START_REF] Sorensen | Newton's method with a model trust region modication[END_REF] With this sign convention, the multipliers associated with the orthogonality constraints in the local step are always nonnegative as shown in Section 4.

Our analysis establishes local, and in some cases global, convergence of the decomposition algorithm of Figure 1 to a stationary point. In Corollary 1 we show that if y = (y 1 , y 2 ) is a local minimizer for (1), then there exist scalars λ 1 , λ 2 , and µ satisfying the rst-order condition

H 1 0 0 H 2 y 1 y 2 = λ 1 I µI µI λ 2 I y 1 y 2 . ( 8
)
where λ 1 and λ 2 lie between the smallest and second smallest eigenvalues of H 1 and H 2 respectively. The condition [START_REF] Strang | Linear Algebra and Its Applications[END_REF] together with the requirement that y 1 and y 2 are feasible in (1) form the KKT (Karush-Kuhn-Tucker) conditions.

Solutions of the subproblems (2) satisfy the following conditions: There exist scalars λ k1 , µ k1 , λ k2 , and µ k2 such that forward

H 1 0 0 H 2 y k1 y k2 = λ k1 y k1 + µ k1 x k2 µ k2 y k1 + λ k2 y k2 , reverse H 1 0 0 H 2 y k1 y k2 = λ k1 y k1 + µ k1 y k2 µ k2 x k1 + λ k2 y k2 . (9)
A fundamental dierence between the rst-order optimality conditions for the original optimization problem (1) and the subproblems (2) is that µ k1 may not equal µ k2 in the subproblems. Hence, a key objective in the analysis is to show that the multipliers in the subproblems approach a common limit. As will be seen in the analysis that follows, we are able to bound the dierence µ k1 -µ k2 in terms of F k (0), which tends to zero. Hence, as F k (0) tends to zero, both the iteration change x k+1 -y k tends to zero, according to [START_REF] Nocedal | Numerical Optimization[END_REF], and the multiplier dierence µ k1 -µ k2 tends to zero.

The local step in the decomposition algorithm requires the solution of a quadratic program of the following form:

min x T Hx subject to x = 1, a T x = 0, (10)
where a ∈ R n with a = 1 and H is symmetric. In the decomposition algorithm, H is H i and a is either x ki or y ki , i = 1 or 2. Since H is symmetric, we can perform an orthogonal change of variables to diagonalize H. Hence, without loss of generality, we can assume that H is diagonal with the ordered eigenvalues

1 ≤ 2 ≤ . . . ≤ n (11)

INRIA

The analysis of the decomposition algorithm is based on an analysis of how the multiplier for the constraint a T x = 0 in (10) depends on a. If H is a multiple of the identity, then this multiplier vanishes, and the dependence of the multiplier on a is trivial. Except in this special case, the dependence of the multiplier on a is nontrivial. For almost every choice of a, the multiplier is unique and depends continuously on a. Suppose that H is not a multiple of the identity and let s denote the smallest eigenvalue of H which is strictly larger than 1 .

The degenerate choices of a, where uniqueness and continuity are lost correspond to those a = 0 which satisfy the equations i= 1

a 2 i s -1 = i> s a 2 i i -s , a i = 0 when i = s . (12)
We Our main result is the following:

Theorem 1. If the global minimizers of (1) are all nondegenerate, then for any starting guess suciently close to the solution set, there exists a subsequence of the iterates of the decomposition algorithm of Figure 1 that approaches a stationary point for [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF].

The proof of Theorem 1 will be given in Section 4.

Remark 1. In the special case where

H 1 = H 2 and 3 -2 ≥ 2 -1 , it is shown in [1]
that the decomposition algorithm is globally convergent for any starting point. On the other hand, we observe in Section 5 that when 3 -2 < 2 -1 , then for specially chosen starting points, the algorithm could converge to a stationary point which is not a global minimum.

In our local convergence result Theorem 1, the requirement for the starting point ensures that the iterates avoid degenerate points for either H The paper is organized as follows. In Section 2 we develop necessary and sucient optimality conditions for a quadratic optimization problem with both a sphere and an ane constraint, and we develop necessary optimality conditions for [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF]. In Section 3 we apply the optimality theory to obtain an optimal solution for the local subproblem, and we show that the multipliers in the subproblems possess a continuity property. The optimality theory also yields the solution to the original problem (1) when H 1 and H 2 commute. In Section 4 we prove our local convergence result Theorem 1. In Section 5 we investigate the global convergence of the decomposition algorithm using a series of numerical examples.

RR n°6288 2 Optimality Conditions

Each step of the domain decomposition algorithm requires the solution of a sphere constrained, quadratic programming problem with a linear constraint. This leads us to consider a problem with the structure

min f (x) := 1 2 x T Hx -h T x subject to x T x = 1, Ax = b, ( 13
)
where A is m by n, h ∈ R n , and b ∈ R m . The local steps of our decomposition algorithm correspond to the case h = 0, m = 1, and b = 0. Our analysis in this section, however, applies to the more general quadratic cost function and linear constraints appearing in (13).

The following result gives necessary and sucient conditions for a point to be a global minimum. Without the linear constraint, this result is known (see [START_REF] Sorensen | Newton's method with a model trust region modication[END_REF]). We give a slightly dierent analysis which also takes into account linear constraints. Recall that at a local minimizer where a constraint qualication holds, the Hessian of the Lagrangian is typically positive semidenite over the tangent space associated with all the constraints, both the linear constraint Ax = b and the sphere constraint x T x = 1. If y is a global minimizer for (13) and λ is the Lagrange multiplier for the sphere constraint, then the second-order necessary optimality condition is that the rst-order condition (8) holds and d T (H -λI)d ≥ 0 whenever Ad = 0 and y T d = 0.

In (15), we claim that the condition y T d = 0 can be dropped and the Hessian of the Lagrangian is positive semidenite over a larger space, the null space of A.

Proposition 1. Suppose that y is feasible in (13). A necessary and sucient condition for y to be a global minimizer is that there exist λ ∈ R and µ ∈ R m such that

Hy = h + yλ + A T µ (14)
and

d T (H -λI)d ≥ 0 whenever Ad = 0. (15) Proof. Let L : R × R m × R n → R be the Lagrangian dened by L(λ, µ, x) = f (x) + λ 2 (1 -x T x) + µ T (b -Ax).
First, suppose that there exist λ ∈ R and µ ∈ R m such that ( 14) and (15) hold. For any feasible x for (13), a Taylor expansion of L around y yields

f (x) = L(λ, µ, x) = L(λ, µ, y) + ∇ x L(λ, µ, y)(x -y) + 1 2 (x -y) T (H -λI)(x -y) = f (y) + 1 2 (x -y) T (H -λI)(x -y). (16)

INRIA

The rst-derivative term in (16) vanishes due to (14). Since x is feasible, A(xy) = 0.

Hence, (15) and ( 16) imply that f (x) ≥ f (y), which shows that y is a global minimizer for (13).

Conversely, suppose that y is a global minimizer for (13). Condition ( 14) is the usual rst-order optimality condition at y. This condition holds if the following constraint qualication is satised (e. g. see [START_REF] Nocedal | Numerical Optimization[END_REF]): For each vector d in the tangent space T at y, there exists a feasible curve approaching y along the direction d, where

T = {d ∈ R n : y T d = 0, Ad = 0}.
Given d ∈ T , such a feasible curve is given by the formula

x(t) = x 0 + y + td -x 0 y + td -x 0 x 0 -y , ( 17
)
where t is a scalar and x 0 is the point satisfying the linear equation Ax = b which is closest to the origin. Since the expression in parentheses in (17) lies in the null space of A and since Ax 0 = b, it follows that Ax(t) = b for each choice of t. Since x 0 is orthogonal to the null space of A, it follows from the Pythagorean theorem that x(t) is a unit vector for each choice of t. Dierentiating x(t), we obtain

x (0) = d -(y -x 0 ) (y -x 0 ) T d y -x 0 2 .
Since d ∈ T , y T d = 0. Since x 0 is orthogonal to the null space of A and Ad = 0, we have x T 0 d = 0. Hence, x (0) = d and there exists a feasible curve approaching y in the direction d. This veries the constraint qualication for (13); consequently, the rst-order condition (14) is satised for some λ ∈ R and µ ∈ R m . By (16), the rst-order optimality condition (14), and the global optimality of y, we have

(x -y) T (H -λI)(x -y) = 2(f (x) -f (y)) ≥ 0 (18)
whenever x is feasible in (13). Suppose that Ad = 0. If in addition, d T y = 0, then d ∈ T . Earlier we observed that when d ∈ T , x(t) is feasible in (13) for all choices of t. Since x(t) is feasible, we can substitute x = x(t) in (18). Since x(t) -y = td + O(t 2 ), it follows from (18), after dividing by t 2 and letting t approach 0 + , that (15) holds. If d T y = 0, then d ∈ T , and y + td < 1 for a suitable choice of t near 0. Increase the magnitude of t until y + td = 1. Substituting x = y + td in (18) gives (15).

We now obtain bounds on the location of the multiplier λ associated with (13).

Proposition 2. If the eigenvalues of H are arranged in increasing order as in (11) and if A has rank k ≥ 1, then λ ≤ k+1 when (15) holds. Moreover, if h = 0 and b = 0 and y is a global minimizer in (13), then λ ≥ 1 .

RR n°6288

Proof. If W is the k + 1 dimensional space spanned by the eigenvectors associated with the k + 1 smallest eigenvalues of H, then we have

k+1 = max{v T Hv : v ∈ W, ||v|| = 1}. ( 19
)
Since A has rank k ≥ 1, the dimension of the null space of A is n -k, and there exists a unit vector v which lies both in the null space of A and in W. Since Av = 0, the second-order condition ( 15) and ( 19) yield

k+1 ≥ v T Hv ≥ λ.
If h = 0 and b = 0, then the rst-order condition (14) implies that

λ = y T Hy ≥ min{v T Hv : v = 1} = 1 .
Next, we focus on the original 2-variable problem (1).

Corollary 1. If (y 1 , y 2 ) is a local minimizer for (1), then there exist λ 1 , λ 2 , and µ such that [START_REF] Strang | Linear Algebra and Its Applications[END_REF] holds. If y is a global minimizer for (1), then for i = 1, 2, we have

λ i ∈ [ i1 , i2 ],
where ij is the j-smallest eigenvalue of

H i , i1 ≤ i2 ≤ . . . ≤ in . (20)
Proof. The gradients of the constraints for (1) at y are multiples of the 3 vectors

y 1 0 , 0 y 2 , y 2 y 1 .
Since these vectors are orthogonal, they are linearly independent. Since the linear independence constraint qualication is satised, the rst-order optimality condition (8) holds for suitable choices of λ 1 , λ 2 , and µ. If y is a global minimizer of (1), then y i is a global minimizer of the problem

min x T H i x subject to x T x = 1, z T x = 0,
where z = y 2 when i = 1 and z = y 1 when i = 1. We apply Proposition 2 with k = 1 to obtain

λ i ∈ [ i1 , i2 ], i = 1, 2.

INRIA 3 The local step and continuity

In each step of the domain decomposition algorithm, we must solve a quadratic programming problem of the form (10). After an orthogonal change of variables, we can assume, without loss of generality, that H is diagonal with the ordered eigenvalues (11) on the diagonal and a = 1. Using Propositions 1 and 2, we now determine the optimal solutions to (10). In the special case h = 0 and A = a T , the rst-order optimality conditions (14) reduce to Hy = yλ + aµ

E i = {j : j = i }. If µ = 0, then by (21) we must have a i = 0 for all i ∈ E 1 . If i ∈ E 1 , then y i = µa i i -1 and a T y = µ i ∈E1 a 2 i i -1 = 0, (22)
which violates the orthogonality condition a T y = 0. Hence, µ = 0 and all y satisfying the following conditions are solutions to (10):

y i = 0 if i ∈ E 1 , a T y = 0, y = 1. (23)
Observe that there is an innite set of solutions y while the multipliers λ and µ are unique.

Case 2: 1 < 2 and a 1 = 0. If λ > 1 , then the second-order condition (15) is violated by the vector d 1 = 1 and d i = 0 for i > 1. Hence, λ = 1 . As in Case 1, the orthogonality condition a T y = 0 is violated unless µ = 0. The solution is again given by (23) and the multipliers are λ = λ 1 and µ = 0.

Case 3: 1 < 2 and a 1 = 0. We rst show that λ > 1 . Suppose, to the contrary, that λ = 1 . The rst component of (21) implies that µ = 0. Hence, (21) reduces to Hy = 1 y. Since H is diagonal and i > 1 for i > 1, we conclude that y i = 0 for i > 1. Hence, y 1 = ±1 since y is a unit vector. However, a vector of this form violates the orthogonality condition a T y = 0 when a 1 = 0. This gives a contradiction, so we have λ > 1 .

(a) a i = 0 for some i ∈ E 2 . We show that λ < 2 . Suppose, to the contrary, that λ = 2 . Since a 1 = 0 and a i = 0 for some i ∈ E 2 , the second-order condition (15) is violated by taking d to be completely zero except for components 1 and i. Since 1 < λ < 2 , (21) can be solved for y:

y = µ(H -λI) -1 a. (24)

RR n°6288

If µ = 0, then y = 0, which violates the constraint y T y = 1. We combine the expression (24), with the orthogonality condition a T y = 0, and the fact that µ = 0 to obtain the equation

g(λ) := n i=1 a 2 i i -λ = 0. ( 25
)
Observe that g is strictly monotone increasing on the interval ( 1 , 2 ) and g( + 1 ) = -∞ since a 1 = 0 while g( - 2 ) = +∞ since a 2 = 0. There exists a unique zero λ of g in ( 1 , 2 ). The solution to (10) is

y i = µa i i -λ where µ 2 = n i=1 a 2 i ( i -λ) 2 -1 = g (λ) -1 . ( 26
)
The equation for µ 2 is obtained from the requirement that y T y = 1. Notice that both the solution y and the multiplier µ are unique to within sign. (b) a i = 0 for all i ∈ E 2 and g( 2 ) < 0. We show that λ = 2 and µ = 0. By (21), we have

y i = µa i i -λ when i ∈ E 2 .
If µ = 0, then the orthogonality condition a T y = 0 reduces to (25), which has no solution on ( 1 , 2 ) since g is monotone on this interval, g( + 1 ) = -∞, and g( 2 ) < 0. Hence, µ = 0. If λ < 2 , then (24) implies that y = 0, which violates the constraint y T y = 1. Hence, λ = 2 and µ = 0. The solution consists of all vectors y satisfying

y i = 0 if i ∈ E 2 , y = 1. (27)
Notice that λ and µ are again unique. (c) a i = 0 for all i ∈ E 2 and g( 2 ) > 0. First, suppose that µ = 0. Since g( + 1 ) = -∞ while g( 2 ) > 0, g in (25) has a unique zero on ( 1 , 2 ). Hence, one solution to (21) is given by (26). We now consider the possibility that µ = 0 at a global minimum. We will show that this leads to a contradiction. Consequently, there is a unique (to within sign) global minimizer for (10) given by (26). If µ = 0, then by ( 21), ( i -λ)y i = 0 for all i, which implies that y i = 0 for i ∈ E 2 since 1 < λ ≤ 2 . Since y = 1, it follows that λ = 2 (or else y = 0, violating the condition y = 1). We now show that the second-order condition (15) is violated for the choice

d i = a i i -γ ,

INRIA

where γ is the unique zero of g on the interval ( 1 , 2 ). This choice for d satises the condition a T d = 0 since g(γ) = 0. Since γ < 2 , we have

d T (H -λI)d = d T (H -2 I)d = n i=1 d 2 i ( i -2 ) = i ∈E2 a 2 i ( i -2 ) ( i -γ) 2 < i ∈E2 a 2 i ( i -γ) ( i -γ) 2 = g(γ) = 0.
This violates the second-order condition (15).

(d) a i = 0 for all i ∈ E 2 and g( 2 ) = 0. This is the degenerate case introduced in Section 1. We rst observe that λ = 2 . Suppose, to the contrary, that λ < 2 . By (21), y is given by (24). If µ = 0, then the orthogonality condition gives (25), which has no solution on ( 1 , 2 ) since g is monotone and g( 2 ) = 0. Consequently, µ = 0 and (24) implies that y = 0, violating the constraint y T y = 1. Thus λ = 2 .

By (21),

y i = µa i i -2 for i ∈ E 2 . (28)
For i ∈ E 2 , the rst-order condition (14) provides no information concerning y i since both sides of the equation vanish identically:

( i -2 )y i = µa i = 0.
The general solution is the following. First, choose any value for y i , i ∈ E 2 , such that i∈E2

y 2 i ≤ 1.
Then choose µ in (28) such that y = 1. In other words, we choose µ so that

µ 2 = 1 -i∈E2 y 2 i g ( 2 )
.

Notice that λ is unique in the degenerate case, while both µ and y are not unique.

Lemma 1. For the optimization problem (10) and a global minimizer y, the multiplier λ associated with the constraint y T y = 1 is a Lipschitz continuous function of a on the unit sphere. With appropriate sign, the corresponding multiplier µ associated with the orthogonality constraint a T y = 0 is continuous at any nondegenerate a.

RR n°6288

Proof. In Case 1, Lipschitz continuity is trivially satised, so we focus on the situation where 1 < 2 . Since the intersection of the hyperplanes a 1 = 0 or a 2 = 0 with the unit sphere are sets of measure zero on the surface of the sphere, Lipschitz continuity over the complement implies Lipschitz continuity over the entire sphere (by continuity). Hence, we restrict our attention to 1 < 2 , a 1 = 0, and a 2 = 0. In this case, λ is the unique solution to (25) on the interval ( 1 , 2 ). Dierentiating (25) gives

∂λ ∂a i = 2a i (λ -i )g (λ)
.

(

) 29
If for some i, we have j∈Ei

a 2 j ≥ 1/2, (30)
then

| i -λ|g (λ) = 1 | i -λ|   j∈Ei a 2 j   + | i -λ| j ∈Ei a 2 j ( j -λ) 2 (31) ≥ 1 | i -λ| j∈Ei a 2 j ≥ 1 2| i -λ| ≥ 1 2( n -1 )
.

It follows from (29) that when (30) holds, ∂λ ∂a i ≤ 4( n -1 ).

If a i = 0, then ∂λ/∂a i = 0 by (29). Now, suppose that a i = 0 and (30) is violated. By (29) and (31), we have

∂λ ∂a i = 2|a i | 1 | i -λ|   j∈Ei a 2 j   + | i -λ|   j ∈Ei a 2 j ( j -λ) 2   ≤ 2 |a i | | i -λ| + | i -λ| |a i |   j ∈Ei a 2 j ( n -1 ) 2   ≤ 2 |a i | | i -λ| + | i -λ| |a i | 1 2( n -1 ) 2 . (32)

INRIA

The last inequality is due to the assumption that (30) is violated, which implies that j ∈Ei a 2 j ≥ 1/2.

|a i | | i -λ| ≥ 1. (33)
By (32), the partial derivative ∂λ/∂a i is bounded by 2 in magnitude since both terms in the denominator of (32) are positive and one of the terms is greater than or equal to 1.

Conversely, if (33) is violated, then we drop the rst term in the denominator of (32) to obtain

∂λ ∂a i ≤ 4( n -1 ) 2 .
At this point, we have shown that there exists a constant β with the property that if a lies on the unit sphere with a 1 = 0 and a 2 = 0, then ∂λ ∂a i ≤ β

for each i. Observe that the solution λ to (25) does not change if a is multiplied by a nonzero scalar. Hence, for any nonzero a (not necessarily on the unit sphere) with a 1 = 0 and a 2 = 0, we have

∂λ ∂a i ≤ β/ a .
Consequently, there exists constants r 1 < 1 < r 2 with the property that ∂λ ∂a i ≤ 2β

whenever r 1 ≤ a ≤ r 2 , a 1 = 0, and a 2 = 0. Given two arbitrary points on the unit sphere, we can construct a piecewise linear path between them with the property that the line segments all lie within the shell formed by the spheres of radius r 1 and r 2 . Due to the bound on the partial derivatives of λ, the change in λ across each line segment is bounded by 2β times the length of the line segment. Since the number of line segments is bounded, independent of the location of the points, we deduce that λ is a Lipschitz continuous function of a on the unit sphere. Now consider the multiplier µ. If 1 = 2 , then µ = 0 (Case 1) and there is nothing to prove. Next, we focus on Case 3a where µ is given by (26). For λ ∈ ( 1 , 2 ), g is bounded away from zero. For example,

g (λ) = n i=1 a 2 i ( i -λ) 2 ≥ 1 ( n -1 ) 2 .

RR n°6288

Let λ(a) denote the unique multiplier associated with any given a. Since λ is a Lipschitz continuous function of a, it follows that µ is continuous at any point a where λ(a) = i for all i. Since λ ∈ [ 1 , 2 ], the only potential points of discontinuity are those points b for which λ(b) = i , i = 1 or i ∈ E 2 . If b 1 = 0 or b i = 0 for any i ∈ E 2 , then µ(a) approaches 0 as a approaches b due to the pole in the denominator of g . Hence, µ is continuous at b.

If b 1 = 0 and λ(b) = 1 , then we use (25) to solve for a 2 1 /( 1 -λ):

a 2 1 λ -1 = i>1 a 2 i i -λ (34)
By assumption, λ(a) approaches 1 as a approaches b. Since the right side of ( 34) is continuous when λ is near 1 , the limit of the left side as a approaches b is

i>1 b 2 i i -1 > 0 since b 1 = 0 and b = 1.
Consequently, by (26), µ tends to 0 as a approaches b, the same limit given in Case 2.

Finally, suppose that b i = 0 for all i ∈ E 2 and λ(b) = 2 . Again, by (25), we have

1 λ -2 i∈E2 a 2 i = i ∈E2 a 2 i i -λ (35)
According to the statement of the lemma, we only need to prove continuity at nondegenerate b, in which case the right side does not vanish at λ = 2 . Hence, as a approaches b, the right side approaches the limit

i ∈E2 b 2 i i -2 = 0.
Consequently, by (26), µ tends to 0 as a approaches b, the same limit given in Case 3b. Note that Case 3c is not a point of discontinuity of µ since λ < 2 . This completes the proof. Now let us consider the original problem [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF]. If H 1 and H 2 commute, then they are simultaneously diagonalizable by the same eigenvector matrix [8, p. 249]. In this case, we can perform an orthogonal change of variables to reduce H 1 and H 2 to diagonal matrices.

The solution is as follows:

Corollary 2. Suppose H i , i = 1 and 2, are diagonal with diagonal element ij , 1 ≤ j ≤ n, arranged in increasing order. The minimum cost in (1) is whichever is smaller. In the rst case, an associated solution to (1) is y 11 = 1, y 22 = 1, and y ij = 0 otherwise.

In the latter case, an associated solution to (1) is y 12 = 1, y 21 = 1, and y ij = 0 otherwise.

Proof. First, let us consider the case where the diagonal elements of H i are strictly separated:

i1 < i2 < . . . < in for i = 1, 2.
By the rst-order optimality conditions [START_REF] Strang | Linear Algebra and Its Applications[END_REF] and by the diagonal structure of the H i , we have ( 1j -λ 1 )y 1j = µy 2j and ( 2j -λ 2 )y 2j = µy 1j .

We combine these equations to obtain

( 1j -λ 1 )( 2j -λ 2 ) -µ 2 y 1j = 0 = ( 1j -λ 1 )( 2j -λ 2 ) -µ 2 y 2j . (36)
By Corollary 1, the multipliers λ 1 and λ 2 satisfy λ i ∈ [ i1 , i2 ]. Hence, the coecients of y 1j and y 2j in (36) are strictly increasing functions of j ∈ [2, n]. It follows that these coecients can vanish for at most one j ∈ [2, n] and possibly for j = 1. When the coecients of y 1j and y 2j do not vanish in (36), we must have y 1j = y 2j = 0. In summary, at the global optimum, all the components of y ij vanish except possibly y 11 , y 21 , y 1j , and y 2j for some j ∈ [2, n]. We focus on the case j = 2 since j > 2 leads to a larger cost. Dene x 2 1j = v j and x 2 2j = w j for j = 1, 2. The optimization problem (1) with H i diagonal and x ij = 0 for j > 2 reduces to

min v 1 11 + v 2 12 + w 1 21 + w 2 22 subject to v 1 + v 2 = 1 = w 1 + w 2 , v 1 w 1 = v 2 w 2 , v 1 , v 2 , w 1 , w 2 ≥ 0.
The equation v 1 w 1 = v 2 w 2 is the orthogonality condition x 11 x 21 = -x 12 x 22 squared. We substitute v 1 = 1 -v 2 and w 1 = 1 -w 2 to reduce the optimization problem to

min 11 + 21 + v 2 ( 12 -11 ) + w 2 ( 22 -21 ) (37) subject to v 2 + w 2 = 1, v 2 ≥ 0, w 2 ≥ 0.
We substitute v 2 = 1 -w 2 in the objective function to further reduce the optimization problem to

min 21 + 12 + w 2 ( 11 + 22 -21 -12 ) subject to 0 ≤ w 2 ≤ 1.

RR n°6288

Since the cost function is linear in w 2 , the minimum is achieved at either w 2 = 0 (w 1 = 1, v 1 = 0, v 2 = 1) with objective function value 21 + 12 or w 2 = 1 (w 1 = 0, v 1 = 1, v 2 = 0) with objective function value 11 + 22 .

When the diagonal elements are not strictly separated, the solution given in the statement of the corollary remains valid. This can be proved as follows: First, perturb the diagonal elements to make them strictly separated. By the previous analysis, we know that the solution given in the statement of the corollary is valid. Next, let the perturbation tend to zero. The limit of these perturbed solutions is a solution of the original unperturbed problem.

Remark 2. Assuming the eigenvalues are all distinct, then the degenerate choices for a in (10) correspond to those vectors a for which

a 2 1 2 -1 = n i=3 a 2 i i -2 , a 2 = 0, i =2 a 2 i = 1. (38)
The solution to (1) given by Corollary 2 has the property that the nonzeros lie in the rst two components of the vectors, while a degenerate a must have nonzero in components greater than or equal to 3. In fact, it follows from (38) that

n i=3 a 2 i ≥ 3 - 2 3 -1 .
Remark 3. Let us consider the special case H 1 = H 2 = H. Since H 1 and H 2 commute, we can apply Corollary 2. Let j denote the j-th smallest eigenvalue of H. As shown in (37), the optimization problem (1) reduces to

min 2 1 + (v 2 + w 2 )( 2 -1 ) subject to v 2 + w 2 = 1, v 2 ≥ 0, w 2 ≥ 0.
Since v 2 + w 2 = 1, the objective function value is 1 + 2 , independent of the choice of v 2 and w 2 satisfying the constraints. Hence, when H 1 = H 2 there are an innite number of solutions to (1). y 1 is any unit vector in the span of the eigenvectors associated with 1 and 2 , and y 2 is any orthogonal unit vector in the same eigenspace. Note that if H is 2 by 2, then all feasible points are optimal and F (x 1 , x 2 ) is the trace of H whenever x 1 and x 2 are feasible in (1).

Convergence of the decomposition algorithm

The proof of Theorem 1 is organized into four steps. In Step 1, we analyze the global step and show that the iteration dierence x k+1 -y k tends to 0. In Step 2, we show that the multipliers µ kj in the local step are almost monotone decreasing since the violation in INRIA monotonicity decays to zero as the iteration number k tends to innity.

Step 3 and 4 focus on the limit of the multipliers µ kj as k tends to innity.

Step 3 considers the limit 0, while

Step 4 considers a positive limit.

Step 1. Analysis of the global step Suppose that iteration k corresponds to the forward mode. Let y k denote the result of the local step, and let x k+1 be the result of the global step based on the starting point y k . Since x k1 is feasible in the rst subproblem of the forward mode (2), we have

F (y k1 , x k2 ) ≤ F (x k1 , x k2 ).
Since x k2 is feasible in the second subproblem, we have

F (y k1 , y k2 ) ≤ F (y k1 , x k2 ).
Combining these relations gives

F (y k ≤ F (x k ). (39)
A similar analysis for the reverse mode also gives F (y k ) ≤ F (x k ).

The components of z k (s) lie on the unit sphere for all choices of s. Consequently, F k (s) is bounded by a nite constant M , uniformly in k and s. Dene the constants δ = min{ρ, 1/M } and sk = -δF k (0).

Since sk lies on the interval [0, -ρF k (0)] appearing in (6), we have

F (x k+1 ) = F k (s k ) ≤ F k (s k ). (40)
Expanding in a Taylor series around s = 0, there exists ξ k ∈ [0, sk ] such that

F k (s k ) = F k (0) + F k (0)s k + 1 2 s2 k F k (ξ k ) ≤ F k (0) + F k (0)s k + 1 2 s2 k M = F k (0) + δF k (0) 2 ( 1 2 δM -1) ≤ F k (0) -δ 2 F k (0) 2 = F (y k ) -δ 2 F k (0) 2 ≤ F (x k ) -δ 2 F k (0) 2 ,
where the last inequality is (39). Combining this with (40) gives

F (x k+1 ) ≤ F (x k ) - δ 2 F k (0) 2 .
Summing this inequality over k yields

F (x k ) ≤ F (x 0 ) - δ 2 k-1 i=0 F i (0) 2 .

RR n°6288

Since the feasible points for (1) lie on the unit sphere, the objective function value is bounded from below. Hence, we have

lim k→∞ F k (0) = 0. (41)
By the denition of z k and the fact that s k ∈ [0, -ρF k (0)] where F k (0) approaches 0, we also conclude that

x k+1 -y k ≤ c|F k (0)| ( 42
)
where c is a constant which is independent of k.

Step 2. The change in the multiplier µ.

By the orthogonality between y k1 and x k2 (forward), between y k1 and y k2 (forward and reverse), and between x k1 and y k2 (reverse), the rst-order optimality conditions (9) for the subproblems (2) yield:

forward            λ k1 = y T k1 H 1 y k1 µ k1 = x T k2 H 1 y k1 λ k2 = y T k2 H 2 y k2 µ k2 = y T k1 H 2 y k2 reverse            λ k1 = y T k1 H 1 y k1 µ k1 = y T k2 H 1 y k1 λ k2 = y T k2 H 2 y k2 µ k2 = x T k1 H 2 y k2 (43)
By our sign convention [START_REF] Sorensen | Newton's method with a model trust region modication[END_REF], the multipliers µ kj are nonnegative. By the denition of F k (s), we have

F k (0) = ±2(y T k1 H 1 y k2 -y T k2 H 2 y k1 ). (44)
We multiply the rst equation in (9) by y T k2 to obtain y T k1 H 1 y k2 = µ k1 y T k2 x k2 . We multiply the second equation by y T k1 to obtain y T k1 H 2 y k2 = µ k2 . Hence, in the forward mode, it follows from (44) that

µ k2 = y T k1 H 2 y k2 = y T k1 H 1 y k2 ∓ F k (0)/2 = µ k1 y T k2 x k2 ∓ F k (0)/2, (45)
which implies that

µ k2 ≤ |µ k1 y T k2 x k2 | + |F k (0)|/2 ≤ µ k1 + |F k (0)|/2 (
µ k1 = y T k-1,2 Hy k1 = y T k1 (µ k-1,2 y k-1,1 + λ k-1,2 y k-1,2 ) = µ k-1,2 y T k1 y k-1,1 .
Since y k1 and y k-1,1 are unit vectors, we deduce that µ k1 ≤ µ k-1,2 . Hence, (51) holds with replaced by ≥. Recall that in our theoretical analysis, the stepsize was restricted to an interval [0, -ρF k (0)], for some xed ρ > 0, to ensure that the iterates approach each other in the limit. However, in all our numerical experiments, we found that there was no need to restrict the stepsize to obtain convergence. Hence, it appears that the restriction on the stepsize in ( 6) is an artifact of the analysis presented in this paper.

In (51) we show that the multipliers associated with the local steps in the decomposition algorithm almost decay monotonically. In Remark 4, we show that the decay is monotone when H 1 = H 2 . Numerically, we found that when H 1 = H 2 , the convergence of the multipliers may not be monotone. An illustration is given in Figure 2 where we randomly generate two 100 by 100 diagonal matrices H 1 and H 2 with entries between -1 and +1, and we plot the multipliers as a function of the iteration number. Due to the initial growth in the multipliers during the rst 300 iterations, the convergence is not monotone and the inequalities in (51) can not be replaced by ≥ in general. In our experiments, the convergence speed when H 1 = H 2 was closely related to the distribution of the smallest 3 eigenvalues of the matrix, independent of the matrix dimension.

To illustrate the typical convergence, we consider the 3 by 3 diagonal matrix

H 1 = H 2 =   1 0 0 0 2 0 0 0 2 + α   , (56)
where α > 0 is a parameter which we vary to explore the convergence. The starting guess is

x T 02 = 1 √ 3 [1 1 1] T . ( 57
Since H 1 = H 2 , the reverse step can be skipped, and decomposition algorithm operates in forward mode without a global step. The components of the iterates are always nonzero, and the iterates given by (26) can be expressed in the form

x k1i = µ k1 x k-1,2i i -λ k1 , x k2i = µ k2 x k,1i i -λ k2 . ( 58
)
The iterates converge to a pair (x 1 , x 2 ) of the form

x T 1 = [a b 0] T and x T 2 = [-b a 0] T with a 2 + b 2 = 1
, which is a valid solution to (1) according to Corollary 2 and Remark 3.

Since the third component of the solution vanishes, we will study how quickly the third component approaches 0. By (58), the ratio between the second and third components can be expressed as Since 2 > 3 for the matrix (56), the rational function ( 2 -λ)/( 3 -λ), λ ∈ [1, 2], attains its maximum value 1/(1 + α) at λ = 1. Hence, by induction, we have

x k13 x k12 = x k-1,23 x k-1,22 1e
x k13 x k12 ≤ 1 1 + α 2k . ( 59
)
Thus as α approaches 0, the bound on the rate at which the third component approaches 0, relative to the second component, grows, and the convergence could be much slower. And as α becomes large, the bound decreases and the convergence rate increases. In summary, the convergence speed seems to depend on the ratio of the gap between 2 and 3 relative to the gap between 1 and 2 . As the ratio approaches 0, the convergence could be slower, as seen in (59).

In Figure 3 we show the convergence of x k13 as a function of the iteration number k for various choices of α. Notice that as α approaches 0, the optimization problem (1) becomes more poorly conditioned since the points

x T 1 = [1 0 0] and x T 2 = [0 0 1]
are feasible with objective function value 3 + α ≈ 3, when α ≈ 0. Thus there are feasible points which are separated from the optimal solution, but with nearly the same cost as the optimal solution. In the extreme case where α = 0, the algorithm nds the global minimum RR n°6288

of the objective function at the rst iteration. The two vectors x 1 and x 2 have their three components dierent from zero. We now give an example where the decomposition algorithm does not converge to the global minimum when the starting guess is suciently poor. In Remark 1, we point out that the decomposition algorithm is convergent for any starting guess when H 1 = H 2 and 3 -2 ≥ 2 -1 . Suppose that 3 -2 < 2 -1 , and the starting guess is

x T 02 = 1 √ 2 [1 0 1 0 0 . . . ] T .
According to case 3c of Section 3,

x T 11 = 1 √ 2
[-1 0 1 0 0 . . . ] T and x 12 = x 02 .

Hence, the algorithm converges after one iteration to the stationary point

x T 1 = 1 √ 2 [-1 0 1 0 0 . . . ] T and x 2 = 1 √ 2 [1 0 1 0 0 . . . ] T .
This starting guess, however, is exceptional. If the second component of x 02 is changed to any nonzero value α, then the iterates quickly converge to the global minimum. For example, if α = 10 -14 and

H 1 = H 2 = diag [-0.9, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, +0.1, +0.2, +0.3],
then the error is reduced to 10 -10 within 44 iterations.

For the case H 1 = H 2 , the decomposition algorithm may converge to a stationary point which is not the global minimum when the starting guess is degenerate and the degenerate iterates are chosen in a very special way. As an example, suppose that H 1 and H 2 are diagonal matrices, with diagonal elements arranged in increasing order, and that all the components of the starting point x 02 are nonzero except for component 2 which is 0. The nonzero components of x 02 are chosen to make it degenerate for H 1 . The initial iterate y 11 is described by case 3d of Section 3. There are an innite number of solutions to the local subproblem. We choose the solution for which the second component is zero (the remaining components are nonzero). Take 11 close enough to 12 to ensure that y 111 is near 1 in magnitude. In this case, g( 22 ) < 0 in the second local step. By case 3b of Section 3, all the components of the iterate y 12 are zero except for the second component which is 1. Since F 1 (0) = 0, the global step has no eect; we have x 21 = y 11 and x 22 = y 12 . Thereafter, x k1 is the rst column of the identity and x k2 is the second column of the identity. If H 1 and H 2 are chosen so that 11 + 22 > 21 + 12 , then the iteration has reached a stationary point which is not the global minimum. In contrast, with any perturbation in the second component of x 02 , we obtain convergence to the global minimum.

INRIA

If the eigenvectors corresponding to the smallest eigenvalues of H 1 and H 2 are orthogonal, then these orthogonal eigenvectors are the solution of [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF]. By randomly choosing the remaining orthogonal eigenvectors of H 1 and H 2 , we obtain noncommuting matrices for which the solution of (1) is known. As a specic example, we took H i = Q i D i Q T i where D i is a diagonal matrix with diagonal elements chosen randomly on [-1, 1], and Q i , i = 1 or 2, is an orthogonal matrix of the form

Q 1 = e 1 U 1 and Q 2 = e 2 U 2 .
Here e i denotes the i-th column of the identity matrix, and U i is an n by n -1 matrix with randomly chosen entries such that Q i is orthogonal, i = This diers from the rst-order optimality conditions (8) associated with the original optimization problem (1) because µ 1 may not equal µ 2 . If the local step is followed by the global step, then according to the analysis of Section 4, F k (0) tends to zero (see (41)), which implies that (see (44)) y T 1 H 1 y 2 = y T 2 H 2 y 1 . Since µ 1 = y T 1 H 1 y 2 and µ 2 = y T 2 H 2 y 1 , the global step ensures that µ 1 = µ 2 . Consequently, the rst-order optimality condition for [START_REF] Barrault | Développement de méthodes rapides pour le calcul de structures électroniques[END_REF] is satised.

The complexity of the analysis is connected with the proof of convergence. To show that the iterate x k2 converges to the same limit as y k2 , we studied the properties of multipliers for the subproblem (10). We showed that the multiplier λ for the sphere constraint lies between 1 and 2 , the two smallest eigenvalues of H; moreover, λ depends Lipschitz continuously on a. In contrast, the multiplier µ associated with the orthogonality constraint is continuous at nondegenerate choices for a.

A less technical explanation for the performance of the decomposition algorithm is that the local steps steer the iterates into a low dimensional subspace associated with the eigenvectors of the smallest eigenvalues of H 1 or H 2 , while the global step nds the best point in this subspace.

RR n°6288

endFigure 1 :

 1 Figure 1: The decomposition algorithm.

1 = 2 .

 12 By Proposition 2, the multiplier λ of Proposition 1 is λ = 1 = 2 . Dene the set

  numerical experiments were performed to investigate the convergence rate of the decomposition algorithm and to explore the connections between the theoretical analysis and the practical convergence. The experiments we describe were performed using Scilab (www.scilab.org). The solution of each local step (10) was obtained by computing an eigenvector associated with the smallest non-zero eigenvalue of the matrix PHP, where P = Iaa T is the projection into the subspace perpendicular to a. The global step was implemented using the Scilab routine optim with default parameter values.

Figure 2 :

 2 Figure 2: Convergence of the multipliers, random 100 by 100 diagonal matrices H 1 and H 2 , ρ = 1.

Figure 3 :

 3 Figure 3: Convergence of x k13 for matrices (56) and starting point (57).

  say that a is degenerate for H if (12) holds, and conversely, a is nondegenerate if (12) is violated or H is a multiple of the identity. The degenerate choices of a compose a set of measure 0. We say that (y 1 , y 2 ) is nondegenerate for (1) if y 1 is nondegenerate for H 2

and y 2 is nondegenerate for H 1 . If H 1 and H 2 commute, then the solution to (1), given in Section 3, is nondegenerate.

  The denominator contains both |a i |/| i -λ| and its reciprocal | i -λ|/|a i |. Suppose that

  46)since y k2 and x k2 are unit vectors. In a similar fashion, for the reverse mode at iteration k + 1, we haveµ k+1,1 ≤ µ k+1,2 + |F k+1 (0)|/2.

	Since y k2 and y k-1,2 are unit vectors, this shows that µ k2 ≤ µ k1 . Multiplying the rst
	equation in (55) by y k-1,2 gives
	(47)
	INRIA

  1 or 2. For all starting points, we observed convergence to the global minimum. Convergence to local, non global, minima has also been observed in the case where the matrices H 1 and H 2 do not commute, but have the same eigenvector associated with the smallest eigenvalue. 6 Conclusions A decomposition algorithm is developed for a quadratic programming problem with sphere and orthogonality constraints. The algorithm consists of local steps, both forward and reverse, and a global step where we minimize over a subspace. Without the global step, any limit (y 1 , y 2 ) of the local step satises

	H 1 0 0 H 2	y 1 y 2	=	λ 1 I µ 1 I µ 2 I λ 2 I	y 1 y 2	.

+ 22 or

+ 21 , INRIA

-λ 1k

-λ 1k = x k-1,13 x k-1,12 2 -λ 2,k-1 3 -λ 2,k-1 2 -λ 1k 3 -λ 1k . INRIA

http://www.inria.fr ISSN 0249-6399

Acknowledgments

Constructive comments by the reviewers are gratefully acknowledged. Their comments signicantly improved the paper.

work supported by the National Science Foundation under Grants 0619080 and 0620286.

If iteration k corresponds to the forward mode, then the multiplier µ k1 corresponds to a = x k2 and H = H 1 in (10). The multiplier µ k-1,1 corresponds to a = y k-1,2 and H = H 1 in (10). By (42) x k2 -y k-1,2 ≤ c|F k-1 (0)|. We apply Lemma 1 and (41). For any (small) η > 0, we have

when k is suciently large, which implies that

The analogous result for the reverse mode is

for k suciently large.

Combining (46)(50), it follows that when k is large enough that |F j (0)| ≤ η for all j ≥ k, we have

where the notation µ k1 µ k2 means that µ k2 ≤ µ k1 + η. Hence, in each iteration, the µ multiplier either decreases or makes an increase which is bounded by η. By (43), the multipliers are bounded by the largest absolute eigenvalues of H 1 and H 2 .

Step 3. The case lim inf µ k1 = 0.

When lim inf µ k1 = 0, there exists a subsequence of the iterates with the property that µ k1 tends to 0. By (51) and the fact that η can be taken arbitrarily small, we conclude that the corresponding subsequence of the multipliers µ k2 also approaches 0. Since y k lies in a compact set, we can extract subsequences converging to a limit that we denote by y. By (43), the corresponding subsequence of multipliers λ k1 and λ k2 also approach limits denoted λ 1 and λ 2 . By (9), we have

Since y 1 and y 2 are orthogonal unit vectors, we conclude that y is a stationary point for (1) corresponding to the multiplier µ = 0.

Step 4. The case µ = lim inf µ k1 > 0.

We extract a subsequence, denoted ν j1 , of the multiplier sequence µ k1 which converges to µ:

Given any η > 0, choose K large enough that (51) holds for all k ≥ K. Also, choose K larger, if necessary, so that |µ -ν j1 | ≤ η for all j ≥ K.

(52)

RR n°6288

Hence, for any j ≥ K, we have µ -η ≤ ν j1 ≤ µ + η.

(53) By ( 51), (53), and the fact that the ν j1 form a subsequence of the µ k1 , it follows that

Let k denote an index in the original sequence with the property that µ k2 = ν j2 . By (51), we have

By ( 48) and ( 52), it follows that

Since η is arbitrary, it follows that ν j1 and ν j2 approach the same limit µ.

Again, by extracting subsequences, there exist limits y 1 , y 2 , λ 1 , and λ 2 such that

By (45), we have

where µ > 0. Since y 2 and x 2 are unit vectors, we deduce that x 2 = y 2 . When x 2 is replaced by y 2 in (54), we see that y is a stationary point. Remark 4. In the special case H 1 = H 2 = H, both the analysis and the algorithm simplify. As noted earlier, F k (0) = 0 in this case so the global step is skipped. Moreover, the monotonicity property (51) holds without the reverse iteration. Hence, the decomposition algorithm can simply employ forward steps, for which the associated rst-order optimality condition is H 0 0 H y k1 y k2 = λ k1 y k1 + µ k1 y k-1,2 µ k2 y k1 + λ k2 y k2 .

(

Multiplying the rst equation by y T k2 gives µ k1 y T k2 y k-1,2 = y T k2 Hy k1 = y T k1 (µ k2 y k1 + λ k2 y k2 ) = µ k2 .

INRIA