Novel reformulations and efficient algorithms for the generalized trust region subproblem
Rujun Jiang Duan Li
math.OC


Electron and Nucleon Localization Functions of Oganesson: Approaching the Thomas-Fermi  Limit
    Witold Nazarewicz
    
=============================================================================================




We present a new solution framework to solve the generalized trust region subproblem (GTRS)  of minimizing a quadratic objective  over a quadratic constraint. More specifically, we derive a convex quadratic reformulation (CQR)  via minimizing a linear objective over two convex quadratic constraints for the GTRS.
We show that an optimal solution of the GTRS can be recovered from an optimal solution of the CQR. We further prove that this CQR is equivalent to minimizing the maximum of the two convex quadratic functions derived from the CQR for the case under our investigation. Although the latter minimax problem is nonsmooth, it is well-structured and convex. We thus develop two  steepest descent algorithms corresponding to two different line search rules. We prove for both algorithms their global sublinear convergence rates. We also obtain a local linear convergence rate  of the first algorithm by estimating the  Kurdyka-Łojasiewicz exponent at any optimal solution under mild conditions.  We finally demonstrate the efficiency of our algorithms in our numerical experiments.







§ INTRODUCTION

We consider the following generalized trust region subproblem (GTRS),

    (P)    min   f_1(x):=1/2x^⊤ Q_1x+b_1^⊤ x
        s.t.   f_2(x):=1/2x^⊤ Q_2x+b_2^⊤ x+c≤ 0,

where Q_1 and Q_2 are n× n symmetric matrices (not necessary to be positive semidefinite),
b_1,b_2∈ℝ^n and c∈ℝ.

Problem (P) is known as the generalized trust region subproblem (GTRS)  <cit.>. When Q_2 is an identity matrix I and b_2=0, c=-1/2, problem (P) reduces to the classical trust region subproblem (TRS). The TRS first arose in the trust region method for nonlinear optimization <cit.>, and  has found many  applications including robust optimization <cit.> and the least square problems <cit.>.
As a generalization, the GTRS also admits its own applications such as time of arrival problems <cit.> and subproblems of consensus ADMM in signal processing <cit.>.  Over the past two decades,
numerous solution methods have been developed for TRS (see
<cit.> and references therein).

Various methods have been developed for  solving the GTRS under various assumptions
(see <cit.> and references therein).
Although it appears being nonconvex, the GTRS essentially enjoys its hidden convexity. The GTRS
can be solved via a semidefinite programming (SDP) reformulation, due to the celebrated S-lemma  <cit.>, which was first established in
 <cit.>. However, suffering from relatively large computational complexity, the SDP algorithm is not practical for large-scale  applications.
To overcome this difficulty, several recent papers <cit.> demonstrated that the TRS admits a second order cone programming (SOCP) reformulation.  Ben-Tal and den Hertog  <cit.> further showed an SOCP reformulation for the GTRS under a  simultaneously diagonalizing (SD) procedure of the quadratic forms.   Jiang et al. <cit.> derived an SOCP reformulation for the GTRS when the problem has a finite optimal value and further derived a closed form solution when the SD condition fails.
On the other hand, there is rich literature on iterative algorithms to solve the GTRS directly under mild conditions, for example, <cit.>.  Pong and Wolkowicz proposed an efficient algorithm based on minimum generalized eigenvalue of a parameterized
matrix pencil for the GTRS, which extended  the results in <cit.> and <cit.> for the TRS.
 Salahi and Taati <cit.>  also derived  a diagonalization-based algorithm under the SD condition of the quadratic forms. Recently, Adachi and Nakatsukasa <cit.> also developed a novel eigenvalue-based algorithm to solve the GTRS.


Our main contribution in this paper is to propose a novel convex quadratic reformulation  (CQR) for the GTRS that is simpler than <cit.>  and further a minimax problem reformulation and develop an efficient algorithm  to solve the minimax problem reformulation. Numerical results demonstrate that our method outperforms all the existing methods in the literature for sparse problem instances. We acknowledge that our CQR was inspired by  the following CQR in Flippo and Janson <cit.> for the TRS,

    min_x{1/2x^⊤(Q_1-λ_min(Q_1)I)x+b_1^⊤ x+1/2λ_min(Q_1):  x^⊤ x≤ 1},

where λ_min(Q_1) is the smallest eigenvalue of matrix Q_1. Unfortunately, this CQR was underappreciated in that time. Recently, people rediscovered this result; Wang and Xia <cit.>
 and Ho-Nguyen and Kilinc-Karzan <cit.> presented a linear time algorithm to solve the TRS by applying Nesterov's accelerated gradient descent algorithm to (<ref>).  We, instead, rewrite  the epigraph reformulation for (<ref>) as follows,

    min_x,t{t: 1/2x^⊤(Q_1-λ_min(Q_1)I)x+b_1^⊤ x+1/2λ_min(Q_1)≤ t,  x^⊤ x≤ 1}.

Motivated by the above reformulation, we demonstrate that the GTRS is equivalent to exact one of the following two convex quadratic  reformulations under two different conditions,

    (P_1)   min_x,t   {t: h_1(x)≤ t,  h_2(x)≤ t},
    (P_2)   min_x,t   {t: h_3(x)≤ t, f_2(x)≤ 0},

where h_1(x), h_2(x) and h_3(x), which will be defined later in Theorem <ref>,   and f_2(x) defined in problem (P),  are  convex but possibly not strongly convex, quadratic functions. To our best knowledge, our proposed CQRs are derived the first time for the GTRS.
The  reformulation (P_2) only occurs when the quadratic constraint is convex and thus can be solved by a slight modification of <cit.> in the accelerated gradient projection method by projecting, in each iteration, the current solution to the ellipsoid instead of the unit ball in the TRS case.

In this paper we focus on the problem reformulation (P_1). Although our CQR can be solved  as an SOCP problem <cit.>, it is not efficient when the problem size is large. Our main contribution is based on a recognition that problem (P_1) is equivalent to minimizing the maximum of the two convex quadratic functions in  (P_1),

    (M) min H(x):=max{h_1(x),h_2(x)}.

 We further derive efficient algorithms to solve the above minimax problem. To the best of our knowledge, the current literature lacks studies on such a problem formulation for a large scale setting except using a black box subgradient method with an O(1/ϵ^2) convergence rate <cit.>, which is really slow. Note that  Section 2.3 in Nesterov's book <cit.> presents a gradient based method with linear convergence rate for solving the minimization  problem (M) under the condition that both h_1(x) and h_2(x) are strongly convex. However, Nesterov's algorithms cannot be applied to solve our problem since in our problem setting at least one function of h_1(x) and h_2(x) is not strongly convex.
By using the special structure of problem (M),  we derive  a steepest descent method in Section 3. More specifically, we choose either the negative gradient when the current point is smooth,  or  a  vector in the subgradient set with the smallest norm (the steepest descent direction) when the current point is nonsmooth  as the descent direction, and derive two steepest descent algorithms with two different line search rules accordingly. In the first algorithm we choose a special step size, and in the second algorithm we propose
a modified Armijo line search  rule.
We also prove the global sublinear convergence rate for both algorithms.
 The
first algorithm even admits a global convergence rate of O(1/ϵ), in the same order as the gradient
descent algorithm, which is faster than the subgradient method. In addition, we demonstrate that the first algorithm also admits  a local linear convergence rate,  by a delicate analysis on the Kurdyka-Łojasiewicz (KL) <cit.> property for problem (M).
We illustrate in our numerical experiments the efficiency of the proposed algorithms  when compared with the state-of-the-art methods for GTRS in the literature.












 



The rest of this paper is organized as follows. In Section 2, we derive an explicit CQR for problem (P) under different conditions and show how to recover an optimal solution of problem (P) from that of the CQR.
In Section 3, we reformulate the CQR to a convex nonsmooth unconstrained minimax problem and derive two efficient solution algorithms. We provide convergence analysis for both algorithms. In Section 4, we demonstrate the efficiency of our algorithms from our numerical experiments. We conclude our paper in Section 5.

Notations  We use v(·) to denote the optimal value of problem (·).  The matrix transpose of matrix A is denoted by A^⊤ and  inverse of matrix A by A^-1, respectively.


§ CONVEX QUADRATIC REFORMULATION

In this section, we derive a novel convex quadratic reformulation for problem (P). To avoid some trivial cases, we assume, w.o.l.g., the Slater condition holds for problem (P), i.e., there exists at least one interior feasible point. When both f_1(x) and f_2(x) are convex,  problem (P) is already a convex quadratic problem.
Hence, w.l.o.g., let us assume that not both f_1(x) and f_2(x) are convex.

We need to introduce the following conditions to exclude some unbounded cases.

The set I_PSD:={λ:Q_1+λ Q_2≽0}∩ℝ_+ is not empty, where ℝ_+ is the nonnegative orthant.

The common null space of Q_1 and Q_2 is  trivial, i.e., (Q_1)∩(Q_2)={0}.




Before introducing our CQR, let us first recall the celebrated S-lemma by defining f̃_1(x)=f_1(x)+γ with an arbitrary constant γ∈ℝ.

The following two statements are equivalent:

1. The system of f̃_1(x)<0 and f_2(x)≤0 is not solvable;

2. There exists μ≥0 such that f̃_1(x)+μ f_2(x)≥0 for all x∈ℝ^n.


Using the S-lemma, the following lemma shows a necessary and sufficient condition under which problem (P) is bounded from below.
Problem (P) is bounded from below if and only if the following system has a solution for λ:

    Q_1+λ Q_2≽0, λ≥0, b_1+λ b_2∈(Q_1+λ Q_2).












We make Assumption <ref> without loss of generality, because otherwise we can prove an unboundedness from below of the problem  (see, e.g., <cit.> and <cit.>). Under Assumption <ref>, if Assumption <ref> fails, there exists no nonnegative λ such that Q_1+λ Q_2≽0 and problem (P) is unbounded from below due to Lemma <ref>.  So both Assumptions <ref> and  <ref> are made without loss of generality.

It has been shown in <cit.> that {λ:Q_1+λ Q_2≽0} is an interval and thus {λ:Q_1+λ Q_2≽0}∩ℝ_+ is also an interval (if not empty). Under Assumptions <ref> and <ref>,   we have the following three cases for I_PSD.

The set I_PSD=[λ_1,λ_2] with λ_1<λ_2.


The set  I_PSD=[λ_3,∞).


The set I_PSD={λ_4} is a singleton.


Note that Condition <ref> occurs only when Q_2 is positive semidefinite. Under Condition <ref>, Q_1 and  Q_2 may not be SD and may have 2×2 block pairs in a canonical form under congruence <cit.>. In this case, when λ is given,  the authors in <cit.> showed how to recover an optimal solution if the optimal solution is attainable, and how to obtain an ϵ optimal solution if the optimal solution is unattainable.
So in the following, we mainly focus on the cases where either Condition <ref> or <ref> is satisfied.Under Condition <ref> or <ref>, problem (P) is bounded from below.


Under Condition <ref> or <ref>, there exists λ_0 such that Q_1+λ_0 Q_2≻0 and λ_0≥0, which further implies b_1+λ_0 b_2∈(Q_1+λ_0 Q_2) as Q_1+λ Q_2 is nonsingular. With Lemma <ref>, we complete the proof.



 §.§ Convex quadratic reformulation for GTRS

It is obvious that problem (P) is equivalent to its epigraph reformulation as follows,

    ( P_0) min{t: f_1(x)≤ t, f_2(x)≤ 0}.

To this end, we are ready to present the main result of this section.
Under Assumption <ref>, by defining h_i(x)=f_1(x)+λ_if_2(x), i=1,2,3, we can reformulate problem (P)  to a convex quadratic problem under Conditions <ref> and <ref>, respectively:


  * Under Condition <ref>, problem (P) is equivalent to the following convex quadratic problem, 
    (P_1) min_x,t{t: h_1(x)≤ t,  h_2(x)≤ t};


  * Under Condition <ref>,  problem (P) is equivalent to the following convex quadratic problem,

    (P_2) min_x,t{t: h_3(x)≤ t,  f_2(x)≤ 0}=min_x{h_3(x):  f_2(x)≤ 0}.




Let us first consider the case where Condition <ref> holds. Due to Lemma <ref>, (P_1) is bounded from below. Together with the assumed Slater conditions,  problem (P_1)  admits the same optimal value as its  Lagrangian dual <cit.>. Due to the S-lemma, problem (P) also has the same optimal value as its Lagrangian dual <cit.>,

    (D) max_μ≥0min_x f_1(x)+μ f_2(x).

Under Condition <ref>, i.e., I_PSD=[λ_1,λ_2] with λ_1<λ_2, it is easy to show that (P_1) is a relaxation of (P_0) since they have the same objective function and the feasible region of (P_1) contains that of (P_0) (note that f_1≤ t and f_2≤0 imply that f_1(x)-t+uf_2(x)≤0 for all u≥0). Thus, 
    v( P_1)≤ v( P_0)=v( P).

The Lagrangian dual problem of (P_1) is

    (D_1) max_s_1,s_2≥0min_x,t t+(s_1+s_2)(f_1(x)-t)+(λ_1s_1+λ_2s_2) f_2(x).



For any primal and dual optimal solution pair (x^*,u^*) of (P) and (D), due to  u^*∈[λ_1,λ_2] as Q_1+μ^* Q_2≽0 from Lemma <ref>, we can always find a convex combination λ_1s̅_1+λ_2s̅_2=μ^* with s̅_1+s̅_2=1. Hence (x^*,s̅,t), with an arbitrary t, is a feasible solution  to  (D_1) and  the objective value of problem (D_1) at  (x^*,s̅,t) is the same with the optimal value of (D). This in turn implies

    v(D_1)≥ v(D).

Since (P_1) is convex and Slater condition is satisfied (because (P_1) is a relaxation of (P) and Slater condition is assumed for (P)), v(P_1)=v(D_1). Finally, by combining (<ref>) and (<ref>), we have v(P_1)=v(D_1)≥ v(D)=v(P)=v(P_0)≥ v(P_1). So all inequalities above become equalities and thus (P_1) is equivalent to (P).

Statement 2 can be proved in a similar way and is thus omitted.


Reformulation (P_2)  generalizes the approaches in <cit.> for the classical TRS with the unit ball constraint to the GTRS with a general convex quadratic constraint.




 To our best knowledge, there is no method in the literature to compute λ_1 and λ_2  in Condition <ref> for general Q_1 and Q_2. However, there exist  efficient methods in the literature to compute   λ_1 and λ_2  when a λ_0 is given such that Q_1+λ_0Q_2≻0 is satisfied. More specifically, the method mentioned in Section 2.4.1 in <cit.> gives a way to compute λ_1 and λ_2: first detect a λ_0 such that Q_0 瞿繙:=Q_1+λ_0Q_2≻0, and then compute λ_1 and λ_2 by some generalized eigenvalues  for a definite matrix pencil that are nearest to 0. Please refer to <cit.> for one of the state-of-the-art methods for detecting λ_0. We can also find another iterative method in Section 5 <cit.> to compute λ_0∈ int(I_PSD) by reducing the length of an interval [λ̅_1,λ̅_2]⊃ I_PSD.
We next report our new method to compute λ_1 and λ_2, which is motivated by <cit.>.
 Our first step is also to find a λ_0 such that  Q_0:=Q_1+λ_0Q_2≻0. Then we compute the maximum generalized eigenvalues for Q_2+μ Q_0 and -Q_2+μ Q_0, denoted by u_1 and u_2, respectively. Note that both  u_1>0 and u_2>0   due to  Q_0≻0 and Q_2 has a negative eigenvalue.  So we have

    Q_1+(1/u_1+λ_0) Q_2≽0  and   Q_1+(-1/u_2+λ_0) Q_2≽0.

Thus Q_1+η Q_2≽0 for all η∈[λ_0-1/u_2,λ_0+1/u_1], which implies λ_1=λ_0-1/u_2 and λ_2=λ_0+1/u_1.
In particular, when one of Q_1 and Q_2 is positive definite, we can skip the step of detecting the definiteness, which would save significant time in implementation.


In fact, when λ_0 is given, we only need to compute one extreme eigenvalues, either λ_1 or λ_2, to obtain our convex quadratic reformulation.  Define x(λ)=-(Q_1+λ Q_2)^-1(b_1+λ b_2) for all λ∈ int(I_PSD) and define γ(λ)=f_2(x(λ)).
 After we have computed λ_0 such that λ_0∈ int(I_PSD), under Assumption <ref>,  we further have Q_1+λ_0Q_2≻0, which makes (Q_1+λ Q_2)^-1 well defined.
In fact, there are Newton type methods in the literature (e.g., <cit.>) for solving the GTRS by finding the optimal λ through  γ(λ)=0.  However, each step in <cit.> involves solving a linear system -(Q_1+λ Q_2)^-1(b_1+b_2), which is time consuming for high-dimension settings. Moreover, the Newton's method does not converge in the so called  hard case[The definition here follows <cit.>. In fact, the definitions of hard case and easy case of the GTRS are similar to those of the TRS. More specifically, if the null space of the Hessian matrix, Q_1+λ^* Q_2, with  λ^* being the optimal Lagrangian multiplier of problem (P), is orthogonal to b_1+λ^*  b_2, we are in the hard case; otherwise we are in the easy case.]. On the other hand, for easy case, an initial λ in  I_PSD is also needed as a safeguard to guarantee the positive definiteness of Q_1+λ Q_2 <cit.>.
It is shown in <cit.> that γ(λ) is either a strictly decreasing function or a constant in int(I_PSD). Following  <cit.>, we have the following three cases: if γ(λ_0)>0, the optimal λ^* locates in [λ_0,λ_2]; if γ(λ_0)=0, x(λ_0) is an optimal solution; and if γ(λ_0)<0, the optimal λ^* locates in [λ_1,λ_0].
Hence we have the following corollary, whose proof is similar to that in Theorem <ref> and thus omitted.
Assume that Assumption <ref> holds and define h_i(x)=f_1(x)+λ_if_2(x), i = 0, 1, 2.
 Under Condition <ref>, the following results hold true.


  * If γ(λ_0)>0, problem (P) is equivalent to the following convex quadratic problem,

    (P_1) min_x,t{t: h_0(x)≤ t,  h_2(x)≤ t}.


  * If γ(λ_0)=0, x(λ_0)=-(Q_1+λ_0 Q_2)^-1(b_1+λ_0 b_2)  is the optimal solution.

  * If γ(λ_0)<0, problem (P) is equivalent to the following convex quadratic problem,

    (P_1) min_x,t{t: h_1(x)≤ t,  h_0(x)≤ t}.


Since both (P_1)  and (P_1) have a similar form to (P_1) and can be solved in a way similar to the solution approach for (P_1), we only discuss how to solve (P_1) in the following.






 §.§ Recovery of optimal solutions

In this subsection, we will discuss the recovery of an optimal solution to problem (P) from an optimal solution to reformulation (P_1).
Before
that, we first introduce the following lemma. Let us assume from now on h_i(x)=1/2x^⊤ A_ix+a_i^⊤ x+r_i, i=1,2. 
If Condition <ref> holds, A_1 and A_2 are simultaneously diagonalizable. Moreover, we have   d^⊤ A_1d>0  for all nonzero vector
d∈(A_2).

Note that Condition <ref> and Assumption <ref> imply that Q_1+λ_1+λ_2/2Q_2≻0, i.e.,  A_1+A_2/2≻0. Let A_0=A_1+A_2/2 and A_0=L^⊤ L be its Cholesky decomposition, where L is a nonsingular symmetric matrix. Also let (L^-1)^⊤ A_1L^-1=P^⊤ DP be the spectral decomposition, where P is an orthogonal matrix and D is a diagonal matrix. Then we have  (L^-1P^-1)^⊤ A_1L^-1P^-1=D and 
    (L^-1P^-1)^⊤ A_2L^-1P^-1=(L^-1P^-1)^⊤ A_0L^-1P^-1-(L^-1P^-1)^⊤ A_1L^-1P^-1=I-D.
Hence A_1 and A_2 are simultaneously diagonalizable by the congruent matrix L^-1P^-1.

Now let us assume S=L^-1P^-1 and thus S^⊤ A_1S=(p_1,…,p_n) and S^⊤ A_2S=(q_1,…,q_n) are both diagonal matrices. Define  K={i:q_i=0,i=1,…,n}.  Since A_1+A_2≻0, p_i>0 for all i∈ K. Let e_i be  the n-dimensional
vector with ith entry being 1 and all others being 0s. We have (Se_i)^⊤ A_1Se_i=p_i>0 for all i∈ K. On the other hand, A_2Se_i=0 for all i∈ K. Hence    d^⊤ A_1d>0 holds for all nonzero vector
d∈(A_2).


From Lemma <ref>, Condition <ref> implies the boundedness of problem (P) and thus the optimal solution is always attainable <cit.>.
 In the following theorem, we show how to recover the optimal solution of problem (P) from an optimal solution of problem (P_1).

 Assume that Condition <ref> holds and x^* is an optimal solution of problem (P_1). Then an optimal solution of problem (P) can be obtained in the following ways:


  * If h_1(x^*)= t and h_2(x^*)≤ t, then x^* is an optimal solution to (P);

  * Otherwise h_1(x^*)<t and h_2(x^*)= t.
For any vector v_l∈(A_2), let θ̃ be a solution of the following equation,
 
    h_1(x^*+θ v_l)=1/2v_l^⊤ A_1v_lθ^2 +(v_l^⊤ A_1x^*  +a_1^⊤ v_l)θ+h_1(x^*)=t.

Then {x̃: x̃=x^*+θ̃v_l,v_l∈(A_2), θ̃ is a solution of (<ref>)} forms the set of optimal solutions of  (P).



Note that at least one of h_1(x^*)≤ t and h_2(x^*)≤ t takes equality. Then we prove the theorem for the following two cases:



  * If h_1(x^*)= t and h_2(x^*)≤ t,  then f_1(x^*)+λ_2f_2(x^*)≤ f_1(x^*)+λ_1f_2(x^*). Hence f_2(x^*)≤0 due to λ_2-λ_1>0.


  * Otherwise, h_1(x^*)<t and h_2(x^*)=t. In this case, for all d∈(A_2) we have d^⊤ A_1d>0 due to Lemma <ref>. We also claim that a_2^⊤ d=0. Otherwise, setting d such that a_2^⊤ d<0  (This can be done since we have a_2^⊤(-d)<0 if a_2^⊤ d>0.) yields

    h_2(x^*+d)=h_2(x^*)+1/2d^⊤ A_2d+(x^*)^⊤ A_2 d+a_2^⊤ d=h_2(x^*)+a_2^⊤ d<t,

where the second equality is due to d∈(A_2) and h_1(x^*+d)<t for any sufficiently small d. This implies that (x^*,t) is not optimal, which is a contradiction.
Equation (<ref>) has two solutions due to the positive parameter before the quadratic term, i.e., v_l^⊤ A_1v_l>0 and the negative constant, i.e., h_1(x^*)-t<0.
With the definition of θ̃,
we know  h_1(x̃)=t and h_2(x̃)=t. This further implies  f_1(x̃)=t and f_2(x̃)=0, i.e., x̃ is an optimal solution to (P).




In Item  2 of the above proof, a_2^⊤ d=0 indicates that problem (P) is in the hard case.




We next illustrate our recovery approach for the following simple     example,

    min{3x_1^2-1/2x_2^2-x_2:-x_1^2+1/2x_2^2+x_2+1≤0}.

Note that, for this example, Condition <ref> holds, λ_1=1 and λ_2=3. Then we have the following CQR,

    min{t:2x_1^2+1≤ t, x_2^2+2x_2+3≤ t}.

An optimal solution of the CQR  is x=(0,-1)^⊤ , t=2. However, this x is  not feasible to (P).
Using the approach in Theorem <ref>, we obtain an optimal solution x̃=(√(2)/2,-1)^⊤   to problem (P). In fact, this instance is in the hard case since the optimal Lagrangian multiplier, λ^*=3, is at the end of the interval {λ:Q_1+λ Q_2≽0, λ≥0} and a-λ^*b∈(Q_1+λ^* Q_2).




We finally point out that our method can be extended to the following variants of GTRS with equality constraint and interval constraint,

    [                                   (EP)                                    min            f_1(x):=1/2x^⊤ Q_1x+b_1^⊤ x                                   (IP)                                    min            f_1(x):=1/2x^⊤ Q_1x+b_1^⊤ x;                                                                          s.t.       f_2(x):=1/2x^⊤ Q_2x+b_2^⊤ x+c=0,                                   s.t. c_1≤ f_2(x):=1/2x^⊤ Q_2x+b_2^⊤ x≤ c_2. ]
   It is shown in <cit.> that (IP) can be reduced to  (EP) with minor computation. It is obvious that all our previous results for inequality constrained GTRS hold for (EP) if we   remove the non-negativity requirement for λ in I_PSD, i.e., I_PSD={λ:Q_1+λ Q_2≽0}.
We thus omit detailed discussion for (EP)  to save space.

In the last part of this section, we compare the   CQR in this paper with CQR for general QCQP in <cit.>.  The authors in <cit.> considered the following general QCQP,

    (QP) minb̃_0^⊤ x     s.t. 1/2x^⊤Q̃ _ix+b̃_ix+c̃_i≤0, i=1,…,m, x∈ X,

where X is a polyhedron.
They further showed that the SDP relaxation of (QP) is equivalent to  the following CQR for (QP):

    (CQP)  minb̃_0^⊤ x     s.t. x∈ G,

where G={x:F_s(x)≤0 for every s∈ T},  F_s(x)=∑_i=1^ms_i(1/2x^⊤Q̃_ix+b̃_ix+c̃_i) and

    T:={s∈ℝ^m: s≥0,τ∈ℝ,  [        ∑_i=1^ms_ic̃_i   1/2(∑_i=1^ms_ib̃_i); 1/2(∑_i=1^ms_ib̃_i^⊤)      ∑_i=1^ms_i/2Q̃_i ]≽0}.
 For the  quadratic problem (P_1), because the variable t is linear in the objective and the constraints, we can reduce T  to

    T:={s:s_1+s_2=1,s≥0,∑_i=1^2s_i/2Q_i≽0,∑_i=1^2s_ib_i∈(∑_i=1^2s_iQ_i)},

where the restriction s_1+s_2=1 does not affect the feasible region G since F_s(x)≤0 is equivalent to  k F_s(x)≤0 with any positive scaling k for s.  Note that   h_1(x)=F_s^1(x) and h_2(x)=F_s^2(x) with s^1=(1,0)^⊤ and s^2=(0,1)^⊤. For any s∈ T,   h_1(x)≤0 and h_2(x)≤0 imply F_s(x)≤0 because
F_s(x) is a convex combination of f_1(x)+λ_1f_2(x) and f_1(x)+λ_2f_2(x). Hence, by the strong duality and with analogous proof to that in Theorem <ref>, the two  feasible regions  of problems (P_1) and (CQP) are equivalent and we further have v(P_1)= v(CQP).













§ EFFICIENT ALGORITHMS IN SOLVING THE MINIMAX PROBLEM REFORMULATION OF THE CQR 


In this section, we propose efficient algorithms to solve the GTRS under Condition <ref>. As shown in Theorem <ref> and Corollary <ref>,  the GTRS is equivalent to (P_1) or either   (P_1) or (P_1). The three problems have similar forms and can be solved by the following proposed method in this section. Hence, to save space, we only consider  solution algorithms for (P_1) in this section.

The convex quadratic problem (P_1)  can be cast as an SOCP problem and solved by many existing solvers, e.g.,  CVX <cit.>, CPLEX <cit.> and MOSEK <cit.>. However, the SOCP reformulation is not very efficient when the dimension is large (e.g., the SOCP solver will take about 1,000 seconds to solve a problem of dimension 10,000).
Fortunately, due to its simple  structure, (P_1) is equivalent to the following  minimax problem of two convex quadratic functions
    (M)   min{H(x):=max{h_1(x),h_2(x)}}.

Hence we aim to derive an efficient method to solve the above minimax problem, thus solving the original GTRS.  Our method is a steepest descent method to find a critical point with 0∈∂ H(x). It is well known that such a critical point is an optimal solution of problem (M).


The following theorem tells us how to find the steepest descent direction.
Let g_1=∇ h_1(x) and g_2=∇ h_2(x). If g_1 and g_2 have opposite directions, i.e., g_1=-tg_2 for some constant t>0 or if g_i=0 and h_i(x)≥ h_j(x) for i≠ j, i,j∈{1,2}, then x is a global optimal solution. Otherwise we can always find the steepest descent direction d in the following way:


  * when  h_1(x)≠ h_2(x),  d=-g_1 if h_1(x)>h_2(x) and otherwise d=-g_2;

  * when h_1(x)=h_2(x), d=-(α g_1+(1-α)g_2), where α is defined in the following three cases:


  * α=0, if  g_1^⊤ g_1≥ g_1^⊤ g_2≥ g_2^⊤ g_2,

  * α=1, if g_1^⊤ g_1≤ g_1^⊤ g_2≤ g_2^⊤ g_2,
  * α=g_2^⊤ g_2-g_1^⊤ g_2/g_1^⊤ g_1+g_2^⊤ g_2-2g_1^⊤ g_2, if g_1^⊤ g_2≤ g_2^⊤ g_2 and  g_1^⊤ g_2≤ g_1^⊤ g_1.




If
 h_1(x)=h_2(x) and g_1=-tg_2, then  0∈∂ H(x). Hence, by the definition of subgradient, we have
 
    H(y)≥ H(x)+0^⊤ (y-x)=H(x), ∀ y,
which further implies that x is the optimal solution.

If g_i=0 and h_i(x)≥ h_j(x) for i≠ j, i,j∈{1,2}, then for all y≠ x, we have H(y)≥ h_i(y)≥ h_i(x)=H(x), i.e.,  x is a global optimal solution.

Otherwise we have the following three cases:


  * When  h_1(x)≠ h_2(x), (suppose, w.l.o.g., h_1(x)>h_2(x)),  for all y∈ℬ(x,δ) with ℬ(x,δ)⊂ {x:h_2(x)<h_1(x)}), H(x)=h_1(x) and thus H(x) is differentiable at x and smooth in its
neighbourhood. Hence, d=-g_1 if h_1(x)>h_2(x). Symmetrically, the case with h_2(x)>h_1(x) can be proved in the same way.

  * When h_1(x)=h_2(x), the steepest descent direction can be found by solving the following problem:

    min_y=1max_g∈∂ H(x)g^Ty.

The above  problem is equivalent to min_g∈∂ H(x)g^2 <cit.>, which is
exactly the following problem in minimizing a quadratic function of α,

    min_0≤α≤1 (α g_1+(1-α)g_2)^⊤(α g_1+(1-α)g_2).

The first order derivative  of the above objective function is g_2^⊤ g_2-g_1^⊤ g_2/g_1^⊤ g_1+g_2^⊤ g_2-2g_1^⊤ g_2. Then if the derivative is in the interval [0,1], the optimal α is given by g_2^⊤ g_2-g_1^⊤ g_2/g_1^⊤ g_1+g_2^⊤ g_2-2g_1^⊤ g_2. Otherwise,
(<ref>) takes its optimal solution on its boundary. In particular, 

  * when
  g_2^⊤ g_2-g_1^⊤ g_2/g_1^⊤ g_1+g_2^⊤ g_2-2g_1^⊤ g_2>1, i.e., g_1^⊤ g_1<g_1^⊤ g_2 and g_2^⊤ g_2>g_1^⊤ g_2,  we have α=1,

  * when
 g_2^⊤ g_2-g_1^⊤ g_2/g_1^⊤ g_1+g_2^⊤ g_2-2g_1^⊤ g_2<0, i.e., g_1^⊤ g_2>g_2^⊤ g_2,   we have α=0.




The above theorem shows that the descent direction at each point with h_1(x)=h_2(x) is either the one with the smaller norm between ∇ h_1(x) and ∇ h_2(x)
 or the negative convex combination d of  ∇ h_1(x) and ∇ h_2(x)
 such that ∇ h_1(x)^⊤ d=∇ h_1(x)^⊤ d. 
 
We next present an example in Figure <ref> to illustrate the necessity of involving the subgradient (in some cases, both gradients are not descent directions). Consider h_1(x)=x_1^2+x_2^2 and h_2(x)=(x_1-1)^2+x_2^2. The optimal solution of this problem is (0.5,0)^⊤. The gradient method can only converge to some point in the intersection curve of h_1(x)=h_2(x), i.e., x_1=0.5, but not the global optimal solution. For example, when we are at x̅=(0.5,0.1)^⊤, the gradients for h_1(x̅) and h_2(x̅) are g_1=(1,0.2)^⊤ and g_2=(-1,0.2)^⊤, respectively. Neither -g_1 nor -g_2 is a descent direction at H(x̅); H(x̅+ϵ g_i)>H(x̅) for any small ϵ>0, i=1,2, due to g_1^⊤ g_2=-0.96<0 and h_1(x̅)=h_2(x̅). (The direction -g_1 is a descent direction, at x̅, for h_1(x)  but ascent for h_2(x) and thus ascent for H(x); the same analysis holds for -g_2.)
The way we use to conquer this difficulty is to choose the steepest descent direction in the subgradient set at  points in the intersection curve. If we use the subgradient direction d=-1/2(g_1+g_2)=-(0,0.2)^⊤, then d is a descent direction since
h_1(x̅+ϵ d)=H(x̅)+2ϵ g_1^⊤ d+ϵ^2 d^⊤ d<H(x̅) and h_2(x̅+ϵ d)=H(x̅)+2ϵ g_1^⊤ d+ϵ^2 d^⊤ d<H(x̅)  for any ϵ with 0<ϵ<2.


Using the descent direction presented in Theorem <ref>, we propose two  algorithms to solve the minimax problem (M), respectively, in Algorithms <ref> and <ref>:
 we first compute a descent direction by Theorem <ref>, apply then two different line search rules for choosing the step size, and finally terminate the algorithm if some termination criterion is met.   The advantage of our algorithms is that each iteration is  very cheap, thus yielding, with an acceptable iterate number, a low cost in CPU time. The most expensive operation  in each iteration is to compute several matrix vector products, which could become cheap  when the matrices are sparse.




 §.§ Line search with a special step size













In the following, we first derive a local linear convergence rate for Algorithm <ref> and then demonstrate a  global sublinear convergence rate for Algorithm <ref>. We analyze the local convergence rate by studying the growth in the neighbourhood of any optimal solution to H(x) in (M). In fact, H(x) belongs to a more general class of piecewise quadratic functions. Error bound and KL property, which are two widely used techniques for convergence analysis, have been studied in the literature,  for several kinds of piecewise quadratic functions, see <cit.>. However, these results are based on piecewise quadratic functions separated by polyhedral sets, which is not the case of H(x).
Li et al. <cit.> demonstrated  that KL property holds for the maximum of finite polynomials, but their KL exponent depends on the problem dimension and is close to one, which  leads to a very weak sublinear convergence rate.
Gao et al. <cit.> studied the KL exponent for the TRS with the constraint replaced by an equality constraint x^⊤ x=1.
However, their technique depends on the convexity of the function x^⊤ x and cannot be applied to analyze our problem. Up to now,  the KL property or error     bound for H(x) has not been yet investigated in the literature related to the linear convergence of optimization algorithms. A significant result of this paper is to estimate the KL exponent of 1/2 for function H(x) when min_x H(x)> max_i{min_x h_1(x),min_x h_2(x)}. With this KL  exponent, we are able to illustrate the linear convergence of our first algorithm with the proposed special step size.

For completeness, we give a definition of KL property in the following. By letting ℬ(x,δ)={y:y-x≤δ}  where · denotes the Euclidean norm of a vector, we have the following definition of KL inequality.
<cit.>
Let f:ℝ^n→ℝ∪{+∞} be a proper lower semicontinuous function satisfying that the restriction of f
to its domain is a continuous function.  The function f is said to have the KL property if for any  ∀ x^*∈{x:0∈∂ f(x)},
there exist C,ϵ>0 and θ∈[0,1) such that 
    Cy≥|f(x)-f^*(x)|^θ,   ∀ x∈ B(x^*,ϵ), ∀ y∈∂ f(x),

where θ is known as the KL exponent.




Under Condition <ref>, we know that there exists λ_0≥0 such that Q_1+λ_0Q_2≻0 and thus  b_1+λ_0b_2∈( Q_1+λ_0Q_2) due to the non-singularity of Q_1+λ_0Q_2. Hence from Lemma <ref>, problem (P) (and thus problem (P_1)) is bounded from below. It is shown in <cit.> that when the two matrices are SD and problem (P) is bounded from below,  the optimal solution of problem (P) is attainable. This further implies that problem (P_1) is  bounded from below with its optimal solution attainable. Assuming that x^* is an optimal solution, the following theorem shows that the KL inequality holds with an exponent of 1/2  at x^* under some mild conditions.

Assume that min h_1(x)<min H(x) and  min h_2(x)<min H(x). Then the KL property in Definition <ref> holds with exponent θ=1/2. 

Note that min h_1(x)<min H(x) and  min h_2(x)<min H(x) imply that, for any x^*∈{x:∂ H(x)=0}, ∇ h_1(x^*)≠0 and ∇ h_2(x^*)≠0, respectively.
Assume L= max{λ_max(A_1),λ_max(A_2)}. We carry out our proof by considering the following two cases.


  * For any point with h_1(x)≠ h_2(x), w.l.o.g., assuming h_1(x)>h_2(x) gives rise to  ∂ H(x)={∇ h_1(x)}. Hence

    |H(x)-H(x^*)|   =   1/2(x-x^*)^⊤ A_1(x-x^*)+(x^*)^⊤ A_1(x-x^*)+a_1^⊤(x-x^*)
       ≤   1/2Lx-x^*^2+∇ h_1(x^*)x-x^*.



On the other hand, ∇ h_1(x)=A_1x+a_1 and

    ∇ h_1(x)^2   =   ∇ h_1(x)-∇ h_1(x^*)+∇ h_1(x^*)^2
       =   (x-x^*)^⊤ A_1A_1(x-x^*)+∇ h_1(x^*)^2+2(∇ h_1(x^*))^⊤ A_1(x-x^*) 
       ≥   ∇ h_1(x^*)^2-2L∇ h_1(x^*)x-x^*.

Define  ϵ_0=min{1,∇ h_1(x^*)/4L}. As ∇ h_1(x^*)≠0, for all x∈ℬ(x^*,ϵ_0),  we then have

    |H(x)-H(x^*)|≤1/2Lϵ_0^2+ ∇ h_1(x^*)ϵ_0≤9/32L∇ h_1(x^*)^2

and

    ∇ h_1(x)^2≥∇ h_1(x^*)^2-2L∇ h_1(x^*)ϵ_0≥1/2∇ h_1(x^*)^2.

Hence |H(x)-H(x^*)|^1/2≤√(9/32L)∇ h_1(x^*)≤3/4√(L)∇ h_1(x). So we have the following
 inequality,

    |H(x)-H(x^*)|^θ≤ C_0y,
 for all y∈∂ H(x) (here {∇ h_1(x)}=∂ H(x)) with θ=1/2, C_0=3/4√(L).


  * Consider next a point x with h_1(x)=h_2(x). Define h_α(x)=α h_1(x)+(1-α)h_2(x), for some parameter α∈[0,1]. Let I={i|(∇ h_1(x^*))_i≠0}.  The optimality condition 0∈∂ H(x^*) implies that there exists some α_0∈[0,1] such that α_0∇ h_1(x^*)+(1-α_0)∇ h_2(x^*)=0.
 Note that  ∇ h_1(x^*)≠0 and ∇ h_2(x^*)≠0 as assumed and thus α_0∈(0,1). Define j= argmax_i{|(∇ h_1(x^*))_i|,i∈ I}, M_1=(∇ h_1(x^*))_j and M_2=(∇ h_2(x^*))_j.
Note that ∂ H_α_0(x^*)=0 implies that α_0 M_1=(1-α_0) M_2.  W.o.l.g, assume  M_1≥ M_2 and thus α_0≤1/2.


Since A_1x (A_2x, respectively) is a continuous function of x, there exists an ϵ_1>0 (ϵ_2>0, respectively) such that for any x∈ℬ(x^*,ϵ_1) (x∈ℬ(x^*,ϵ_2), respectively), 3/2M_1≥|(∇ h_1(x))_j| >1/2M_1 (3/2M_2≥|(∇ h_2(x))_j| >1/2M_2, respectively).
Let ϵ_3=min{ϵ_1,ϵ_2}. Then we have the following two subcases.


  * For all x∈ℬ(x^*,ϵ_3) and α∈[0,1/4α_0], we  have

    ∇ h_α(x)   ≥   -α |(∇ h_1(x))_j|+(1-α) |(∇ h_2(x))_j|
       ≥   -3/2α M_1+1/2(1-α)M_2
       ≥   -3/8α_0 M_1+3/8(1-α_0)M_2+(1/8+1/8α_0)M_2
       =   (1/8+1/8α_0)M_2.

The third inequality is due to the fact that -3/2α M_1+1/2(1-α)M is a decreasing function of α and the last equality is due to α_0 M_1=(1-α_0) M_2.
Symmetrically,  for α∈[1-1-α_0/4,1], we have |(∇ h_α(x))|≥(3/8-1/4α_0)M_1.
Combining these two cases and  α_0≤1/2
yields ∇ h_α(x)≥1/8M_2.

On the other hand

    |H(x)-H(x^*)|   =   1/2(x-x^*)^⊤ A_1(x-x^*)+(x^*)^⊤ A_1(x-x^*)+a_1^⊤(x-x^*)
       ≤   1/2Lx-x^*^2+∇ h_1(x^*)x-x^*
       ≤   (1/2Lϵ_3^2+∇ h_1(x^*))x-x^*.

Letting ϵ_4=min{ϵ_3,1} leads to M_2^2/32Lϵ_3^2+64∇ h_1(x^*)|H(x)-H(x^*)|≤∇ h_α(x)^2. So 
    |H(x)-H(x^*)|^θ≤ C_1∇ h_α(x), ∀α∈[0,1/4α_0]∪[1-1-α_0/4,1], ∀ x∈ℬ(x^*,ϵ_4)
where θ=1/2 and C_1=√(32Lϵ_3^2+64∇ h_1(x^*))/M_2.


  * Next let us consider the case with α∈[α_0/4,1-1-α_0/4].
In this case,
 defining A_α=α A_1+(1-α)A_2
and a_α=α a_1+(1-α)a_2 gives rise to

    ∇ h_α(x)^2   =   ∇ h_α(x)-∇ h_α(x^*)+∇ h_α(x^*)^2
       =   (x-x^*)^⊤ A_α A_α(x-x^*)+∇ h_α(x^*)^2+2(∇ h_α(x^*))^⊤ A_α(x-x^*)

and since h_1(x)=h_2(x) and h_1(x^*)=h_2(x^*),

    |H(x)-H(x^*)|   =   1/2(x-x^*)^⊤ A_α(x-x^*)+(x^*)^⊤ A_α(x-x^*)+a_α^⊤(x-x^*)
       =   1/2(x-x^*)^⊤ A_α(x-x^*)+(∇ h_α(x^*))^⊤(x-x^*).

Define μ_0=λ_min(A_α). Then

    ∇ h_α(x)^2-2μ_0|H(x)-H(x^*)|
       =   (x-x^*)^⊤ A_α (A_α-μ_0 I)(x-x^*)+∇ h_α(x^*)^2
          +2(∇ h_α(x^*))^⊤(A_α-μ_0 I)(x-x^*) 
       =   (A_α-μ_0 I)(x-x^*)+∇ h_α(x^*)^2+μ_0(x-x^*)^⊤(A_α-μ_0 I)(x-x^*)≥0.


We next show that μ_0 is  bounded from below. Define μ_1 (μ_2, respectively) as the smallest nonzero eigenvalue of A_1 (A_2, respectively).
Note that α A_1+(1-α)A_2 is positive definite for all α∈[α_0/4,1-1-α_0/4] as assumed in Condition <ref>. Then A_1 and A_2 are simultaneously diagonalizable as shown in Lemma <ref>. Together with the facts that A_1≽0 and A_2≽0, there exists a nonsingular matrix P such that P^⊤ A_1P=D_1≽μ_1(δ ) and P^⊤ A_2P=D_2≽μ_2(δ ), where δ_i=1 if D_ii>0 and δ_i=0 otherwise.
Since α∈[α_0/4,1-1-α_0/4], λ_min(A_α)≥min{αμ_1, ㄗ(1-α)μ_2}≥min{α_0/4μ_1,1-α_0/4μ_2}>0.
From  ∇ h_α^2-2μ_0|H(x)-H(x^*)|≥0, we know ∇ h_α^2-μ_0|H(x)-H(x^*)|≥0.

Let θ=1/2, C_2=√(1/(2μ_0)). We have

    C_2∇ h_α(x)≥|H(x)-H(x^*)|^θ,
     ∀α∈[α_0/4,1-1-α_0/4], x∈ℬ(x^*,ϵ_4).


Combining cases (a) and (b) gives rise to 
    |H(x)-H(x^*)|^θ≤ C_3∇ h_α(x)
 with θ=1/2, C_3=max{C_1,C_2}, for all x∈ℬ(x^*,ϵ_4).

Combining cases 1 and 2 yields that the KL inequality holds  with θ=1/2 and C=max{C_0,C_3}   for all x∈ℬ(x^*,ϵ) with ϵ=min{ϵ_0,ϵ_4}.


Note that the assumption  min h_1(x)<min H(x) and  min h_2(x)<min H(x) means that we are in the easy case of GTRS as in this case λ^* is an interior point of I_PSD and  Q_1+λ^* Q_2 is nonsingular, where λ^*  is the optimal Lagrangian multiplier of the GTRS <cit.>. However, there are two situations for the hard case. Let us consider the KL property for H(x) at the optimal solution x^*. When h_i(x^*)>h_j(x^*), for i=1 or 2 and j={1,2}/{i}, in the neighbourhood x^*,   H(x) is just h_i(x), and the KL is also 1/2 <cit.>. In such a case, our algorithm performs asymptotically like the gradient descent method for unconstrained quadratic minimization.
However, when  h_i(x^*)=h_j(x^*) (note that min h_j(x)<H(x^*) can still hold in this situation),
 the KL exponent is not always 1/2 for H(x).  Consider the following counterexample with
h_1(x)=x_1^2 and h_2(x)=(x_1+1)^2+x_2^2-1. The optimal solution is (0,0) and is attained by both h_1 and h_2.
Let x_2=-ϵ, where ϵ is a small positive number. Consider the curve where h_1(x)=h_2(x), which further implies x_1=-ϵ^2/2.
Then we have

    (1-β ) ∇ h_1+β∇ h_2=2[ -(1-β) ϵ^2/2+β(-ϵ^2/2+1);                       βϵ ]=2[ - ϵ^2/2+β;        βϵ ],

and thus

    min_y∈∂ H(x)y^2   =   min_β4(β^2ϵ^2+β^2-ϵ^2β+ϵ^4/4)
       =   min_β4((1+ϵ^2)(β-ϵ^2/2(1+ϵ^2))^2-ϵ^4/4(1+ϵ^2)+ϵ^4/4)
       =   ϵ^6/2(1+ϵ^2)=ϵ^6/2+O(ϵ^8).

Thus, min_y∈∂ H(x)y= O(ϵ^3). On the other hand,

    H(x)-H(x^*)=x_1^2=ϵ^4/4.

The KL inequality cannot hold with θ=1/2, but it holds with  θ=3/4 since min_y∈∂ H(x)y= O(ϵ^3) and
H(x)-H(x^*)=O(ϵ^4).



It  is interesting to compare our result with a recent result on KL exponent of the quadratic sphere constrained optimization problem <cit.>. In <cit.>, the authors showed that the KL exponent is 3/4 in general and 1/2 in some special cases, for the following problem, 
    (T)   min1/2x^⊤ Ax+b^⊤ x     s.t.  x^⊤ x=1.

The above problem is equivalent to the TRS when the constraint of the TRS is active, which is the case of interest in the literature. For the TRS, the case that the constraint is inactive is trivial: Assuming x^* being the optimal solution, (x^*)^⊤ x^*<1 if and only if the objective function is convex and the optimal solution of the unconstrained quadratic function 1/2x^⊤ Ax+b^⊤ x locates in the interior of the unit ball. The authors in   <cit.> proved that the KL exponent is 3/4 in general and particularly the KL exponent is 1/2 if A-λ^*I is nonsingular, where λ^* is the optimal Lagrangian multiplier. The later case is a subcase of  the easy case  for the TRS and the case that KL exponent equals 3/4 only occurs in some special situations of the hard case.  On the other hand, our result shows the KL exponent is  1/2  for the minimax problem when the associated GTRS is in the easy case.
So our result can be seen as an extension of the resents on KL exponent for  problem (T) in <cit.>. One of our future research is to verify if the KL exponent is 3/4 for H(x) when the  associated GTRS is in the hard case.


For convergence analysis with error bound or KL property, we still need a sufficient descent property to achieve the convergence rate. We next propose an algorithm with such a property. We further show that our algorithm converges locally linearly with the descent direction  chosen  in Theorem <ref> and the step size specified  in the following theorem.


Assume that the conditions in Theorem <ref> hold and that the initial point x^0∈ℬ(x^*,ϵ).  Assume that h_i is the active function when h_1(x_k)≠ h_2(x_k) and h_j, j={1,2}/{i}, is thus inactive.  Let the descent direction   be chosen in Theorem <ref> and the associated step size be chosen as follows.


  * When  h_1(x_k)= h_2(x_k),


  * if  there exists   g_α=α∇ h_1(x_k)+(1-α)∇ h_2(x_k) with α∈[0,1] such that ∇ h_1(x_k)^⊤ g_α=∇ h_2(x_k)^⊤ g_α, then set d_k=-g_α and β_k=1/L, where L=max{λ_max(A_1),λ_max(A_2)};

  * otherwise  set d_k=-∇ h_i(x_k) for i such that  ∇ h_1(x_k)^⊤∇ h_2(x_k)≥∇ h_i(x_k)^⊤∇ h_i(x_k), i=1,2, and β_k=1/L.



  * When h_1(x_k)≠ h_2(x_k) and the following quadratic equation for γ,

    𝐚x^2+𝐛x+𝐜=0,

where 𝐚=1/2γ^2∇ h_i(x_k)^⊤ (A_i-A_j)∇ h_i(x_k), 𝐛=(∇ h_i(x_k)^⊤ -∇ h_j(x_k)^⊤ )∇ h_i(x_k) and 𝐜= h_i(x_k)-h_j(x_k),
 has no positive solution or any positive solution γ≥1/L, set d_k=-∇ h_i(x_k) with  and β_k=1/L;

  * When h_1(x_k)≠ h_2(x_k) and the quadratic equation (<ref>)  has a positive solution γ<1/L,  set β_k=γ  and d_k=-∇ h_i(x_k).

 Then  the sequence {x_k} generated by Algorithm <ref> satisfies, for any k≥1,

    H(x_k)-H(x^*)≤(√(%s/%s)2C^2L-12C^2L)^k-1(H(x^0)-H(x^*)),

and  
    dist(x_k,X)^2≤2/L(H(x_k)-H(x^*)≤2/L(√(%s/%s)2C^2L-12C^2L)^k-1(H(x^0)-H(x^*)) .


For simplicity, let us denote g_i=∇ h_i(x_k)  for i=1,2. We claim the following sufficient descent property for  steps 1,  2 and 3:

    H(x_k)-H(x_k+1)≥L/2x_k-x_k+1^2.

  Hence,
if the step size is 1/L (i.e., steps 1 and 2), we have
    H(x_l)-H(x^*)≤ Cd_l^2= C^2L^2x_l-x_l+1^2≤ 2C^2L(H(x_l)-H(x_l+1)),

where the first inequality is due to the KL inequality in Theorem <ref>, the second equality is due to x_l+1=x_l-1/Ld_l and the last inequality is due to the sufficient descent property. Rearranging the above inequality yields

    H(x_l+1)-H(x^*)≤2C^2L-1/2C^2L(H(x_l)-H(x^*)).

And since our method is a descent method, we have
H(x_l+1)-H(x^*)≤ H(x_l)-H(x^*) for all iterations. Suppose that there are p iterates of step size 1, q iterates of step size 2, and r iterates of step size 3.
From the definitions of the steps, every step 3 is followed by  a step 1 and thus r≤ p+1 if we terminate our algorithm at step 1 or 2. So for all k≥1, after k=p+q+r steps, we have

    H(x_k)-H(x^*)≤(2C^2L-1/2C^2L)^p+q(H(x^0)-H(x^*))≤(2C^2L-1/2C^2L)^k-1/2(H(x^0)-H(x^*)).


The sufficient descent property further implies that 
    L/2∑_k^∞x_k-x_k+1^2≤ H(x_k)-H(x^*).
 Hence, with ∑_k^∞x_k-x_k+1^2≥ dist(x_k,X)^2, we have L/2 dist(x_k,X)^2≤ H(x_k)-H(x^*).
Thus  
    dist(x_k,X)^2≤2/L(H(x_k)-H(x^*)).


By noting g_i=A_ix_k+a_i, we have

    h_i(x_k+1)-h_i(x_k)   =    1/2(x_k+d_k)^⊤ A_i(x_k+d_k)+a_i^⊤(x_k+d_k)-[1/2(x_k)^⊤ A_ix_k+a_i^⊤ x_k]
       =   1/2d_k^⊤ A_id_k+(A_ix_k+a_i)^⊤ d_k
       =   1/2d_k^⊤ A_id_k+g_i^⊤ d_k.

We next prove our claim (<ref>) according to the three cases in our updating rule:

  * When h_1(x_k)= h_2(x_k), noting that h_i  is active at x_k+1 as assumed,  we have
    H(x_k)-H(x_k+1)= h_i(x_k)-h_i(x_k+1).



  * If there exists an α such that g_α^⊤ g_1=g_α^⊤ g_2, we have g_α^⊤ g_i=g_α^⊤ g_α. And by noting that d_i=-g_α, we further  have

    h_i(x_k+1)-h_i(x_k)   =   1/2L^2d_k^⊤ A_id_k+1/Lg_i^⊤ d_k
       ≤   1/2Lg_α^⊤ g_α-1/Lg_α^⊤ g_α
       =   -1/2Lg_α^⊤ g_α.

Substituting g_α=L(x_k-x_k+1) to the above expression, we have the following sufficient descent property,

    H(x_k)-H(x_k+1)=h_i(x_k)-h_i(x_k+1)≥L/2x_k-x_k+1^2.


  * If there does not exist an α such that g_α^⊤ g_1=g_α^⊤ g_2, then we must have g_1^⊤ g_2>0. And thus we must have g_1^⊤ g_1≥ g_1^⊤ g_2≥ g_2^⊤ g_2 or g_2^⊤ g_2≥ g_1^⊤ g_2≥ g_1^⊤ g_1. If g_i^⊤ g_i≥ g_i^⊤ g_j≥ g_j^⊤ g_j, we set d_k=-g_j.
Then

    H(x_k+1)-H(x_k)   ≤   max{h_i(x_k+1)-h_i(x_k), h_j(x_k+1)-h_j(x_k)}
       ≤   max{1/2L^2g_j^⊤ A_ig_j-1/Lg_i^⊤ g_j, 1/2L^2g_j^⊤ A_ig_j-1/Lg_j^⊤ g_j}
       ≤   max{1/2L^2g_j^⊤ A_ig_j-1/Lg_j^⊤ g_j, 1/2L^2g_j^⊤ A_ig_j-1/Lg_j^⊤ g_j}
       ≤   max{1/2Lg_j^⊤ g_j-1/Lg_j^⊤ g_j, 1/2Lg_j^⊤ g_j-1/Lg_j^⊤ g_j}
       =   -1/2Lg_j^⊤ g_j=-L/2x_k-x_k+1^2.


Symmetrically, if g_i^⊤ g_j>0 and g_j^⊤ g_j≥ g_i^⊤ g_j≥ g_i^⊤ g_i, setting  d_k=-g_i yields the same sufficient descent property.


  * When h_1(x_k)≠ h_2(x_k) and the quadratic equation (<ref>) for γ has no positive solution or has a positive solution γ≥1/L, we have h_i(x_k+1)>h_j(x_k+1)  for x_k+1=x_k+β_kd_k, where d_k=-∇ h_i(x_k)  and β_k=1/L. Moreover,

    H(x_k+1)-H(x_k)   =    h_i(x_k+1)-h_i(x_k)
       =   1/2L^2g_i^⊤ A_ig_i-1/Lg_i^⊤ g_i
       ≤   -1/2Lg_i^⊤ g_i.

Hence H(x_k)-H(x_k+1)≥1/2Lg_i^⊤ g_i≥L/2x_k-x_k+1^2.

  * When h_1(x_k)≠ h_2(x_k) and the quadratic equation (<ref>)  has a positive solution γ<1/L. With β_k=γ  and d_k=-∇ h_i(x_k),
 it is easy to see that the step size γ makes h_1(x_k+1)=h_2(x_k+1). Then we have

    H(x_k+1)-H(x_k)   =   h_i(x_k+1)-h_i(x_k)
       =   1/2γ^2d_k^⊤ A_id_k+γ g_i^⊤ d_k
       ≤   1/2Lγ^2g_i^⊤ g_i-γ g_i^⊤ g_i
       =    (L/2-1/γ)x_k-x_k+1^2,

 which further implies H(x_k)-H(x_k+1)≥L/2x_k-x_k+1^2 due to γ≤1/L.



It is worth to note that Step 3 in our algorithm is somehow similar to the retraction step in manifold optimization <cit.>. In manifold optimization, in every iteration, each point is retracted to the manifold. In Step 3, every  point  is drawn to the curve that h_1(x)=h_2(x).


We will next show that in general a global sublinear convergence rate, in the same order with the gradient descent algorithm,  can also be theoretically guaranteed for Algorithm 1.


Assume that x^* is an optimal solution. Then we have 
    H(x_N)-H(x^*)≤L/Nx_0-x^*^2.

That is, the required iterate number for H(x_N)-H(x^*)≤ϵ is at most O(1/ϵ). 

From the proof in Theorem <ref>, for any step size γ≤1/L, we have

    H(x_k+1)-H(x_k)≤-γ g^⊤ g+1/2Lγ^2g^⊤ g≤-γ/2g^⊤ g.

From the convexity of H(x) and g∈∂ H(x_k), we have

    H(x_k+1)   ≤    H(x_k)-γ/2g^⊤ g
       ≤   H(x^*)+g^⊤(x_k-x^*)-γ/2g⊤ g
       =   H(x^*)+1/2γ(x_k-x^*^2-x_k-x^*-γ g^2)
       =   H(x^*)+1/2γ(x_k-x^*^2-x_k+1-x^*^2).

Since H(x_k+1)≥ H(x^*), we have x_k-x^*^2-x_k+1-x^*^2≥0. Let us use  indices i_k, k=0,…,K to denote the indices in Steps 1 and 2. By noting that γ=1/L, we have

    H(x_i_k+1)≤ H(x^*)+L/2(x_i_k-x^*^2-x_i_k+1-x^*^2).
 Note that every Step 3 is followed by S tep 1. Hence N≤ 2K+1.  Adding the above inequalities from i_0 to i_K, we have

    ∑_k=0^KH(x_i_k)-H(x^*)
       ≤   L/2∑_k=0^K(x_i_k-x^*^2-x_i_k+1-x^*^2)
       ≤   L/2(x_i_0-x^*^2-x_i_K+1-x^*^2+∑_k=1^K(-x_i_k-1+1-x^*^2+x_i_k-x^*^2))
       ≤   L/2(x_i_0-x^*^2-x_i_K+1-x^*^2)
       ≤   L/2x_i_0-x^*^2
       ≤   L/2x_0-x^*^2,

where in the second inequality we use the fact,

    -x_i_k-1+1-x^*^2+x_i_k-x^*^2≤-x_i_k-1+1-x^*^2+x_i_k-1-x^*^2≤⋯≤0 .
 Since H(x_k) is non-increasing, by noting that N≤ 2K+1, we have

    H(x_N)-H(x^*)   ≤   1/K+1∑_k=0^KH(x_i_k)-H(x^*)
       ≤   L/Nx_0-x^*^2.






 §.§ Line search with the modified Armijo rule

An alternative way to choose the step size in the classical gradient descent type methods  is the line search with the Armijo rule. A natural thought is then to extend the Armijo rule in our minimax problem (M) as  in the proposed Algorithm <ref>.
In particular, we set the following modified Armijo rule to choose the smallest nonnegative integer k such that the following inequality holds for the step size β_k=ξ s^k with 0<ξ≤1 and 0<s<1,

    f(x_k+β_k p_k)≤ f(x_k)+σβ_k p_k^⊤ g,

where 0≤σ≤0.5, g=-d and  d is the steepest descent direction defined in Theorem <ref>. Particularly, we set the search direction p_k=d at iterate k. Our numerical result in the next section shows that Algorithm <ref> has a comparable performance when compared with  (or even better than) Algorithm <ref>.
For the sake of completeness, we present the convergence result  for Algorithm <ref> in the following.
 Before that, we generalize the definition of a critical point to a (ρ,δ) critical point.

A point x is called a (ρ,δ) critical point of H(x)=max{h_1(x),h_2(x)} if  ∃g<δ, for some g∈∂ H_ρ(x), where ∂ H_ρ(x) is defined as follows:

  * ∂ H_ρ(x)={α∇ h_1(x)+(1-α)∇ h_2(x):α∈[0,1]}, if |h_1(x)-h_2(x)|≤ρ;

  * ∂ H_ρ(x)={∇ h_1(x)}, if h_1(x)-h_2(x)> ρ;

  * ∂ H_ρ(x)={∇ h_2(x)}, if h_2(x)-h_1(x)> ρ.  

The following proposition shows the relationship of a critical point  and  a (ρ,δ) critical point. As this result is pretty obvious, we omit its proof.

Assume that  {x_k} is a sequence in ℝ^n and that (ρ^t,δ^t)→(0,0), for t→∞  and that there exists a positive integer  K(t), such that  x_k is a (ρ^t,δ^t) critical point of H(x) for all k≥ K(t) and t≥1. Then, every accumulation point of the  sequence  {x_k} is  a critical point of H(x).


Slightly different from Algorithm <ref>, our goal in Algorithm <ref> is to find a (ρ,δ) critical point. With Proposition <ref>, we conclude that Algorithm <ref> outputs a solution that is sufficiently close to a critical point of H(x).


Assume that  i) d= argmin_y∈∂  H_ρ(x_k)y with ρ>0,ii) the termination criterion is d<δ for some δ>0 and iii) x^* is an optimal solution. Then for any given  positive numbers ρ and δ,   Algorithm <ref> generates a (ρ,δ) critical point in at most

    H(x_0)-H(x^*)/σ smin{1/L,ξ,ρ/2G^2}δ^2

iterations, where G is some positive constant only depending on the initial point and problem setting.


 Consider the following different cases with d≥δ.


  * If |h_1(x_k)-h_2(x_k)|<ρ, then as assumed  d>δ and  from Theorem <ref>, we know that d= argmin_α∈[0,1]α∇ h_1(x_k)+(1-α)∇ h_2(x_k) is just the parameter α which we choose in Algorithm <ref>.
It suffices to show that the step size β_k is bounded from below such that

    H(x_k+1)-H(x_k)≤ -σβ_kd^⊤ d.
 This further suffices
to show that  β_k is bounded from below such that for i=1 or 2,

    h_i(x_k+1)-h_i(x_k)= -β_k∇ h_i(x_k)^⊤ d+1/2β_k^2d^⊤ A_id≤ -σβ_kd^⊤ d.

By noting that ∇ h_i^⊤ d≥ d^⊤ d from Remark <ref>, the  second inequality in (<ref>) holds true for all β_k≤2(1-σ)/L.  Then the step size chosen by the modified Armijo rule satisfies β_k≥ smin{2(1-σ)/L,ξ}, which further implies that

    H(x_k)-H(x_k+1)≥σβ_kg^⊤ g=σβ_k g^2≥σ smin{2(1-σ)/L,ξ}δ^2.


  * If h_1(x_k)-h_2(x_k)>ρand ∇ h_1(x_k)>δ, we have g=∇ h_1(x_k).
Because H(x_k) is decreasing, under Condition <ref>, h_1(x_k)+ h_2(x_k)=1/2x_k^⊤ (A_1+A_2)x_k+(a_1+ a_2)^⊤ x_k≤ 2h_1(x_k)=2H(x_k)≤2H(x_0) and thus x_k is bounded due to A_1+A_2=2(Q_1+λ_1+λ_2/2Q_2)≻0. This further implies that ∇ h_i(x_k)=A_ix_k+b_i  is bounded for all k. So there exists some positive constant only depending on the initial point and problem parameters such that ∇ h_i(x_k)≤ G, i=1,.  Hence d≤ G because d is a convex combination of ∇ h_1(x_k) and ∇ h_2(x_k).
Then we have

    h_1(x_k+1)-h_1(x_k)≤ -β_k∇ h_1(x_k)^⊤ d+1/2β_k^2d^⊤ A_1d

and for any β_k≤1/L,
    h_2(x_k+1)-h_2(x_k)   ≤    -β_k∇ h_2(x_k)^⊤ g+1/2β_k^2g^⊤ A_ig
       ≤   β_kGg+1/2β_k^2Lg^2
       ≤   β_kG^2(1+1/2β_kL)
       ≤   3/2β_kG^2.

On the other hand,  when β_k≤1/L,

    h_1(x_k+1)-h_1(x_k)≤ -β_k∇ h_1(x_k)^⊤ g+1/2β_k^2g^⊤ A_1g≤ -β_kg^⊤ g+1/2β_k^2Lg^⊤ g≤-1/2β_kg^⊤ g.

Note that for all β_k≤ρ/2G^2,    3/2β_kG^2+1/2β_kg^⊤  g≤ρ.
Thus for   β_k≤min{1/L, ρ/2G^2},
we have

    h_1(x_k+1)≤ h_1(x_k)-1/2β_kg^⊤ g,
          h_2(x_k+1)≤ h_2(x_k)+3/2β_kG^2≤ h_1(x_k)-ρ+3/2β_kG≤ h_1(x_k)-1/2β_kg^⊤ g.

Hence we have

    H(x_k+1)-H(x_k)   =   max{h_1(x_k+1),h_2(x_k+1)}-h_1(x_k)
       =   max{h_1(x_k+1)-h_1(x_k), h_2(x_k+1)-h_1(x_k)}
       ≤   max{h_1(x_k+1)-h_1(x_k), h_2(x_k+1)-h_2(x_k)}
       ≤   -1/2β_kg^⊤ g.

So the Armujo rule implies  β_k≥ smin{1/L,ξ,ρ/2G^2}, i.e., β_k is lower bounded. Then according to the modified Armijo rule, we have

    H(x_k)-H(x_k+1)≥σβ_kg^⊤ g≥σ smin{1/L,ξ,ρ/2G^2}δ^2.


  * Symmetrically, the case with h_2(x_k)-h_1(x_k)>ρ  yields the same result as in (<ref>).

The above three cases show that H(x_k)-H(x_k+1)≥σ smin{1/L,ξ, ρ/2G^2}δ^2 (as 1-σ≥1/2, the decrease in case 1 also admits this bound). Since the decrease in each iterate is larger than σ smin{1/L,ξ,ρ/2G^2}δ^2, the total iterate number is bounded by

    H(x_0)-H(x^*)/σ smin{1/L,ξ,ρ/2G^2}δ^2.


At the current stage, we cannot demonstrate a theoretical convergence rate for Algorithm <ref>  as good as the sublinear rate O(1/ρ) for Algorithm 1 in Theorem <ref>. But our numerical tests show that Algorithm <ref>  converges as fast as Algorithm <ref>. Proposition <ref> and Theorem <ref> offer our main convergence result for Algorithm <ref> as follows.

Assume that (ϕ_k,ψ_k)→0 and  that {x^(k)} is a sequence of  solutions generated by Algorithm <ref> with ρ=ϕ_k and δ=ψ_k. Then
any accumulation point of {x^(k)} is an optimal solution of problem (M).




§ NUMERICAL TESTS

In this section, we illustrate the efficiency of our algorithm with numerical experiments. All the numerical tests were implemented in
Matlab 2016a, 64bit and were run on a Linux machine with 48GB
RAM, 2600MHz cpu and 64-bit CentOS release 7.1.1503.  We compare both Algorithms <ref> and <ref>  with the ERW algorithm in <cit.>.
We  disable the parallel setting in the Matlab for fair comparison. If the parallel setting is allowed, our algorithm has a significant improvement, while the ERW algorithm does not.

We use the following same  test problem as <cit.> to show the efficiency of our algorithms,
    (IP)   min   x^⊤ Ax-2a^⊤ x
        s.t   c_1≤ x^⊤ Bx≤ c_2,

where A is an n× n  positive definite matrix and B is an  n× n (nonsingular) symmetric indefinite matrix.
We first reformulate problem (IP) to  a formulation of problem (P) in the following procedure, which is motivated from <cit.> (the proof in <cit.> is also based on the monotonicity of γ(λ), which is defined in Section 2.1), in order to  apply the CQR for problem (P) and  then invoke Algorithms <ref> and <ref> to solve the CQR.

  Let x_0=-A^-1a. Then
the followings hold.  

  * If x_0^⊤ Bx_0<c_1, problem (IP) is equivalent to

    (IP_1)   min{ x^⊤ Ax-2a^⊤ x: s.t. c_1≤ x^⊤ Bx};


  * Else if c_1≤ x_0^⊤ Bx_0≤ c_2, problem (IP)  admits an interior solution x_0;

  * Otherwise c_2<x_0^⊤ Bx_0,  problem (IP) is equivalent to
    (IP_2)   min{ x^⊤ Ax-2a^⊤ x: s.t.  x^⊤ Bx≤ c_2}.



Item 2 is obvious. Item 1 and Item 3 are symmetric. So in the following, we only prove Item 1.

In our problem set, matrix A is positive definite and B is indefinite. Hence, in the definition I_PSD={λ:Q_1+λ Q_2≽0}, we have λ_1<0, λ_2>0. Thus from Case 1 in Section 2.2.2 in <cit.> we know, when x_0^⊤ Bx_0<c_1,  problem (IP) is equivalent to

    ( EP_1)   min{ x^⊤ Ax-2a^⊤ x: s.t. c_1=x^⊤ Bx}.

Since  x_0^⊤ Bx_0<c_1, the optimal solution of ( IP_1) must be at its boundary <cit.>. This further yields  that problem (IP) is equivalent to ( IP_1).


Theorem 4.1 helps us  solve problem (IP) as an inequality constrained GTRS instead of solving two GTRS with equality constraints.
Before showing the numerical results, let us illustrate some functions used in our initialization. To obtain the CQR, the generalized eigenvalue problem is solved  by   in Matlab, which was developed in <cit.> for computing the maximum generalized eigenvalues for sparse definite matrix pencils.
In our numerical setting  is usually faster than the Matlab function , though  will outperform   when the condition number is large or the density is low. We use the Matlab command  and  to generate Q_1 and Q_2. We set the density of matrices  at 0.01 and use three levels of condition number for matrix Q_1, i.e., 10, 100 and 1000 and, in such settings,    always dominates   (this may be because  is developed for computing extreme generalized eigenvalues for arbitrary matrices and does not utilize the definiteness and symmetry properties of the matrix pencils in our problem setting).
In general, the main cost in estimating L is to compute the maximum eigenvalues of matrices A_1 and A_2, which may be time consuming for large-scale matrices. To conquer this difficulty, we can  estimate a good upper bound with very cheap cost instead. Specially, we can run the function  with precision 0.1, which is much more efficient than  computing the true maximum eigenvalue with  , and, assuming M is the output,  M+0.1 is then a good upper bound for L. In our numerical tests, we just use   to estimate L since our main goal is to illustrate the efficiency of Algorithm 2. In Algorithm <ref>, to avoid some numerical accuracy problem, we approximate  h_1(x_k)=h_2(x_k) by |h_1(x_k)-h_2(x_k)|/(|h_1(x_k)|+|h_2(x_k)|)≤ϵ_1. Also we use |h_1(x_k)-h_2(x_k)|/(|h_1(x_k)|+|h_2(x_k)|)≤ϵ_1 instead of  |h_1(x_k)- h_1(x_k)|≤ρ in Algorithm 2 for stableness consideration. In our numerical tests for both Algorithms <ref> and <ref>, we use the following termination criteria  (if any one of the following three conditions is met, we terminate our algorithm), which are slightly different from the presented algorithms for robust consideration: 

  * H(x_k-1)-H(x_k)<ϵ_2,

  * |h_1(x_k)-h_2(x_k)|/(|h_1(x_k)|+|h_2(x_k)|)≤ϵ_1, α∇ h_1(x_k)+(1-α)∇ h_2(x_k)≤ϵ_3,

  * ∇ h_i(x_k)≤ϵ_3 and |h_1(x_k)-h_2(x_k)|/(|h_1(x_k)|+|h_2(x_k)|)>ϵ_1, where i≠ j and i,j∈{1,2},ㄛ

where ϵ_1, ϵ_2 and ϵ_3>0 are some small positive numbers for termination of the algorithm.   Particularly, we set   ϵ_1=10^-8, ϵ_2=10^-11 and ϵ_3=10^-8 in Algorithm 1, and ϵ_1=10^-8, ϵ_2=10^-11, ϵ_3=10^-8,  σ=10^-4 and ξ=1 (for the modified Armijo rule) in Algorithm 2.


To improve the accuracy of the solution, we apply the Newton refinement process in Section 4.1.2 in <cit.>. More specifically, assuming x^* is the solution returned by our algorithm, we update x^* by

    δ=(x^*)^⊤ Bx^*/2Bx^*^2Bx^*,   x^*=x^*-δ.
In general, the ERW algorithm can achieve a higher precision than our method (after the Newton refinement process); the precision in their method is about 10^-14, while ours is slightly less precise than theirs. Letting v_1 denote the optimal value of ERW algorithm and v_2 denote the optimal value of our algorithm, we  have at least |v_2-v_1|/|v_1|≈10^-10 for most cases.
The iteration number reduces to 1/5 if we reduce the precision of from ϵ_1=10^-8, ϵ_2=10^-11, ϵ_3=10^-8 to ϵ_1=10^-5, ϵ_2=10^-8, ϵ_3=10^-5. This observation seems reasonable as our method is just a first order method.
 

We report our numerical results in Table 1.  We use “Alg1” and “Alg2” to denote Algorithms <ref> and <ref>, respectively.   For each n and each condition number, we generate 10 Easy Case and 10 Hard Case 1 examples.
Please refer to Table 1 in <cit.> for the detailed definitions of Easy Case and Hard Cases 1 and 2. There is a little difference about the definitions of easy and hard cases between <cit.> and  <cit.>. Our analysis in the above sections uses the definitions in  <cit.>. In fact, the Easy Case and Hard Case 1 are the easy case  and Hard Case 2 is the hard case mentioned in the above sections and <cit.>. We use the notation “time" to denote the average CPU time (in unit of second) and “iter” to denote the average iteration numbers for all the three algorithms. For “Alg1” and “Alg2”, “time" is  just the time for Algorithms <ref> and <ref>, respectively. The notation “time_eig" denotes the average CPU time for computing the generalized eigenvalue for our algorithm. So the total time for solving problem (P) should be the summation of the time of reformulate (P) into (M) and the time of Algorithm 1 or 2, whose main cost is just “time"+“time_eig". And “fail"  denotes the failure times in the 10 examples in each case for the ERW algorithm. One reason of the failures may be that the ERW algorithm terminates in 10 iterations even when it does not find a good approximated solution. We point out that for randomly generated test examples, our method always succeeds in finding an approximated solution to prescribed precision while  the ERW algorithm fails frequently in Hard Case 1.   Another  disadvantage of the ERW algorithm is the requirement of an efficient prior estimation of the initialization, which is unknown in general. In our numerical test, we assume that such an initialization is given as the same as <cit.> does.

We also need to point out that in the Hard Case 2,  our algorithms do not outperform the ERW algorithm which uses the  shift and deflation technique.
 The main time cost of shift and deflate operation
is the computation of the extreme generalized eigenvalue of the matrix pencil (A,B) and its corresponding generalized eigenvectors.
In the test instances, as the dimension of the eigenspace of the extreme generalized eigenvalue is  one, the   shift and deflation technique directly finds the optimal solution by calling   once. Our algorithm reduces to an unconstrained quadratic minimization in Hard Case 2. However, the condition number of this unconstrained  quadratic minimization is so large that our algorithm performs badly as the classical gradient method. To remedy this disadvantage, we can add a step with almost free-time cost that claims that   either  we are in Hard Case 2 and output an optimal solution or we are in Easy Case or Hard Case 1. Recall that the hard case (or equivalently, Hard Case 2)  states that b_1+λ^*b_2 is orthogonal to the null space of Q_1+λ^*Q_2 which means that  λ^*  must be a boundary point of I_PSD. Suppose λ_i=λ^*. Then we must have  that  x^*=min H(x) and H(x^*)=h_i(x^*) for some i=1 or 2.   In fact, if ∇ h_i(x)=0 and h_i(x)≥ h_j(x), j∈{1,2}/{i} for some x, then x is optimal and we are in the hard case. So ∇ h_i(x)=0 and h_i(x)≥ h_j(x) is sufficient and necessary for x to be optimal to (M)  and be in the hard case. Hence we can construct an optimal solution for problem (M) as  x̅=(Q_1+λ_i Q_2)^†(b_1+λ _ib_2)+∑_i^k α_j v_j  (where A^† denotes the Moore–Penrose pseudoinverse of A) if v_j,j=1,…,k  are the generalized eigenvectors of matrix pencil (Q_1,Q_2) with respect to the generalized eigenvalue λ_i such that h_i(x̅)≥ h_j(x̅) and α≥0. This equals to identifying if a  small dimensional convex quadratic programming problem (with variable α) has an optimal value less than h_i((Q_1+λ_i Q_2)^†(b_1+λ _ib_2)). And if such α does not exist, we are in the easy case (or equivalently, Easy Case or Hard Case 1). This technique is very similar to the shift and deflation technique in <cit.>. Hence we can  solve Hard Case 2 within almost the same CPU time as the ERW algorithm. So we do not make further comparison for Hard Case 2.


Our numerical tests show that  both Algorithms 1 and 2 are much more efficient than the ERW algorithm in
Easy Case and for most cases in Hard Case 1.  The efficiency of our algorithms is mainly due to that  we only call the generalized eigenvalue solver once and every iteration only involves several matrix vector products (which are very cheap for sparse matrices).  We also note that, in Easy Case, Algorithm 1 is faster than Algorithm 2 when the condition number is small and slower than Algorithm 2 when the condition number is large. This may be because that  Algorithm 2  is equipped with the modified Armijo rule, which makes it more aggressive in choosing the step size and thus yields a fast convergence. In Hard Case 1,
Algorithm 2 is still much more efficient than the ERW algorithm while Algorithm 1 is slower than the ERW algorithm in about half the cases. This is because Algorithm 2 has a moderate iterate number due to the aggressiveness in choosing the step size and Algorithm 1 has a much large iterate number for these cases.  Moreover, our algorithms always succeed, while the ERW algorithm fails frequently in Hard Case 1.
A more detailed analysis with condition number for Algorithm 1 will be given in the following.

We note that several examples (of the 10 examples) in Easy Cases admit a much larger iteration number than average. This motivates us to analyze the main factor that affects the convergence rate (reflected by the iteration number) of Algorithm 1 (the analysis for Algorithm 2 seems hard  due to the non-smoothness of the problem). We then find that the main factor  is √(λ_maxα/2λ_min nnzα^2), as evidenced by the fact that examples in Easy Case and Hard Case 1 with more iterates all have a larger √(λ_maxα/2λ_min nnzα^2), where  λ_maxα denotes the maximum eigenvalue of matrix α A_1+(1-α)A_2 and λ_min nnzα denotes the smallest nonzero eigenvalue of matrix α A_1+(1-α)A_2  with α being defined in Theorem <ref> in the last iteration. In fact, when x^k→ x^*∈{x:∂ H(x)=0} (in our examples, the optimal solution is unique), let the value of α at iterate k be α^k, then
α^k→α^*, where α^*  is the solution of α∇ h_1(x^*)+(1-α)∇ h_2(x^*)=0.
From the definition of KL exponent,  we have

    C×min_αα∇ h_1(x^k)+(1-α)∇ h_2(x^k)≥|H(x^k)-H(x^*)|^1/2.

Intuitively, the smallest value of  C should be at least 
    |H(x^k)-H(x^*)|^1/2/min_αα∇ h_1(x^k)+(1-α)∇ h_2(x^k)→|α (h_1(x^k)-h_2(x^*))+(1-α) (h_1(x^k)-h_2(x^*))|^1/2/min_αα∇ h_1(x^k)+(1-α)∇ h_2(x^k)
 which is upper bounded by √(λ_maxα/2λ_min nnzα^2). Thus, the asymptotic value of C can be  roughly seen as √(λ_maxα/2λ_min nnzα^2). Hence both Easy Case and Hard Case 1 admit local linear convergence and the convergence rate is

    (√(1-1/2C^2L))^k= (√(1-λ_min nnzα^2/Lλ_maxα))^k

from Theorem <ref>.
We also observe  from our numerical tests that in most cases the values of   λ_maxα are similar and that λ_min nnzα in Easy Case is much larger than λ_min nnzα in Hard Case 1 and λ_maxα in Easy Case is very close to λ_maxα  in Hard Case 1. Hence,   √(1-λ_min nnzα^2/(Lλ_maxα)) in Easy Case  is usually much smaller than that in Hard Case 1.   (As Q_2 is random in our setting, the larger the condition number of Q_1 is, the larger expectation of  √(1-λ_min nnzα^2/(Lλ_maxα)) is.) This  explains why the condition number of matrix Q_1 measures, to a large degree, the hardness of our algorithms in solving problem (M).
  Since Easy Case has a smaller √(1-(λ_min nnzα^2/Lλ_maxα)) than Hard Case 1 for the same condition number and problem dimension, Easy Case can be solved faster than Hard Case 1.
   This coincides
with our numerical results, i.e., Easy Case admits a smaller iterate number than  Hard Cases 1.


We also tried to apply MOSEK <cit.> to solve the CQR. But our numerical results showed that MOSEK is much slower than both our algorithms and the ERW algorithm, which took about 833 seconds for Easy Case and 960 second for Hard Case 1 with   n=10000 and cond=10. So we do not run further numerical experiments with MOSEK. We also tested the SOCP reformulation <cit.> under the simultaneous digonalization condition of the quadratic forms of the GTRS and the DB algorithm in <cit.> based on the simultaneous digonalization condition of the quadratic forms. The simultaneous digonalization condition naturally holds for problem (IP) when A is positive definite. Our preliminary result shows that our method is much more efficient than the two methods based on  simultaneous digonalization when n≥10000 and density=0.01 and thus we also do not report the numerical comparison in this paper. We believe this is mainly because the  simultaneously digonalization procedure of the matrices involves  matrix inverse, matrix matrix product, a full Cholesky decomposition and a spectral decomposition (of a dense matrix),  which is more time consuming than the operations of  matrix vector products in our algorithm.  Hence we do not report the numerical results based on the simultaneous digonalization technique.




§ CONCLUDING REMARKS

In this paper, we have derived a simple convex quadratic reformulation for the GTRS, which only involves a linear objective function and two convex quadratic constraints under mild assumption.
 We further reformulate the   CQR to an unconstrained minimax problem under Condition <ref>, which is the case of interest. The minimax reformulation is a well structured convex, albeit non-smooth, problem.  By investigating its inherent structure, we have proposed two efficient matrix-free algorithms to solve this minimax reformulation. Moreover, we have offered a theoretical guarantee of global sublinear convergence rate for both algorithms and  demonstrate a local linear convergence rate for Algorithm 1 by proving  the KL property for the minimax problem with an exponent of 1/2 under some mild conditions.
Our numerical results have demonstrated clearly out-performance of our algorithms over the state-of-the-art algorithm for the GTRS.

As for our future research, we would  like to show whether the CQR and the minimax reformulation and the algorithms for the minimax problem can be extended to solve GTRS  with additional linear constraints.
As  the analysis in numerical section indicates that our algorithms have similar performance with unconstrained quadratic minimization, i.e., both algorithms admit a locally linear convergence rate with the steepest descent method, we would like to generalize existing algorithms that are efficient in solving unconstrained quadratic minimization to solve our minimax reformulation, e.g., the conjugate gradient method or Nesterov's accelerated gradient descent algorithm.
Another line of future research is to investigate whether our algorithm can be extended to general minimax problems with more (finite number of) functions. It is also interesting to verify whether the KL property still holds and whether the KL exponent is still 1/2 when more functions are involved.



§ ACKNOWLEDGEMENTS

This research was partially supported by Hong Kong Research Grants Council under
Grants 14213716 and 14202017. The second author is also grateful to the support from Patrick Huen
Wing Ming Chair Professorship of Systems Engineering and Engineering Management. The authors would also like to thank Zirui Zhou and Huikang Liu  for their insightful discussions.

abbrv

