SIAM J. OPTIM.
Vol. 27, No. 1, pp. 269–291

c 2017 Society for Industrial and Applied Mathematics

SOLVING THE TRUST-REGION SUBPROBLEM
BY A GENERALIZED EIGENVALUE PROBLEM∗SATORU ADACHI† , SATORU IWATA† , YUJI NAKATSUKASA‡ , AND AKIKO TAKEDA§
Abstract. The state-of-the-art algorithms for solving the trust-region subproblem (TRS) arebased on an iterative process, involving solutions of many linear systems, eigenvalue problems, subspace optimization, or line search steps. A relatively underappreciated fact, due to Gander, Golub,
and von Matt [Linear Algebra Appl., 114 (1989), pp. 815–839], is that TRSs can be solved by
one generalized eigenvalue problem, with no outer iterations. In this paper we rediscover this factand discover its great practicality, which exhibits good performance both in accuracy and efficiency.
Moreover, we generalize the approach in various directions, namely by allowing for an ellipsoidal constraint, dealing with the so-called hard case, and obtaining approximate solutions efficiently whenhigh accuracy is unnecessary. We demonstrate that the resulting algorithm is a general-purpose
TRS solver, effective both for dense and large-sparse problems, including the so-called hard case.
Our algorithm is easy to implement: its essence is a few lines of MATLAB code.
Key words. trust-region subproblem, generalized eigenvalue problem, elliptic inner product,
hard case
AMS subject classifications. 49M37, 65K05, 90C25, 90C30
DOI. 10.1137/16M1058200

1. Introduction. The trust-region subproblem (TRS) [4], [26, Chap. 4] is to(1)1minimize g > p + p> Ap,
2
kpkB ≤∆

where A, B ∈ Rn×n are symmetric, g ∈ Rn , ∆ > 0, and B is symmetric positive
definite. Trust-region methods are a popular approach to dealing with general nonlinear optimization problems to minimize f (x), in which each iteration requires an(approximate) solution for TRS (1). In a trust-region method, the objective functionof TRS (1) is a quadratic model of f near the current approximate solution x̃, inwhich A is the Hessian and g is the gradient of f at x̃. Note that we allow
pfor the constraint kpkB ≤ ∆ in an ellipsoidal norm, defined by kpkB = kB 1/2 pk = p> Bp for apositive definite matrix B, not necessarily equal to the identity I. An appropriate andnontrivial choice B 6= I can be important, for example, when working in a properly
scaled trust region to solve the nonlinear problem efficiently [23]. For example, in[4] it is argued that a good choice would be to take B = |A|, the Hermitian polar
factor [17] of A. For more on trust-region methods and TRS, see the book [4].
∗ Received by the editors January 25, 2016; accepted for publication (in revised form) November
21, 2016; published electronically February 22, 2017.
http://www.siam.org/journals/siopt/27-1/M105820.html
Funding: The work of the second, third, and fourth authors was supported by Grants-in-Aid for
Scientific Research (Challenging Exploratory Research, 26540007) from Japan Society for Promotionof Science. The work of the third author was supported by Japan Society for Promotion of Science
as Postdoctoral Fellow for Research Abroad.
† Department of Mathematical Informatics, The University of Tokyo, Tokyo 113-8656, Japan(satoru adachi@mist.i.u-tokyo.ac.jp, iwata@mist.i.u-tokyo.ac.jp).
‡ Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK (nakatsukasa@maths.ox.ac.
uk).
§ Department of Mathematical Analysis and Statistical Inference, The Institute of Statistical
Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan (atakeda@ism.ac.jp).

269

270

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDAThe necessary and sufficient optimality condition for TRS (1) is the following
theorem.
Theorem 1.1. A vector p∗ is an optimal solution to the TRS (1) if and only ifthere exists λ∗ ≥ 0 such that(2)kp∗ kB ≤ ∆,(3)

(A + λ∗ B)p∗ = −g,(4)λ∗ (∆ − kp∗ kB ) = 0,

(5)

A + λ∗ B  0.

Equations (2)–(5) are mentioned in [4, Thm. 7.4.1] and [15] as a necessary condition. Theorem 1.1 is well known when B = I, and it can be obtained by using achange of variables in the discussion in Gay [11] and Moré and Sorensen [25] (see also[26, Chap. 4]) to allow for a general B  0.
Most existing algorithms for the TRS focus on B = I and attempt to find a TRS
solution (λ∗ , p∗ ) in Theorem 1.1 by an iterative process, during which a parametersuch as an estimate for λ∗ is updated. Algorithms for TRS can be classified into threebroad categories as we summarize below, together with the representative studies.
1. Accurate methods for dense problems. The classical algorithm by Moré andSorensen [25] iteratively solves symmetric positive-definite linear systems via
the Cholesky factorization. During the iteration the estimate for λ∗ is adjusted. Safeguard techniques are sometimes necessary to ensure convergenceto the solution. This is the standard approach for moderate-sized problems,
say n ≤ 1000.
2. Accurate methods for large-sparse problems. Sorensen’s algorithm [35], refinedand implemented by Rojas, Santos, and Sorensen [32, 33], iteratively com >
putes the smallest eigenvalue of a parameterized matrix of the form αg gA ,where α is a parameter adjusted during the iteration to find the solution.
Again, a safeguard technique is needed to guarantee convergence. Anotherapproach due to Rendl and Wolkowicz [30] solves TRS via semidefinite programming (SDP), for which the standard solver based on the interior-point
method also involves an iteration of linear systems [2]. A modification of theMoré–Sorensen algorithm using Taylor series is presented by Gould, Robinson, and Thorne [15], and other accurate algorithms include subspace projection methods, such as Hager [16] and Erway and Gill [6].
3. Approximation methods. The Steihaug–Toint algorithm [36, 39] is a wellknown method that employs a truncated conjugate gradient step (see [4,sect. 7.5]). Also known is the approach by Gould et al. [13] using truncatedLanczos steps. These algorithms allow one to pursue the practical goal of
reducing the overall cost of solving the nonlinear optimization problem via
TRS, in which a sufficient reduction in the objective value (1) is usuallyadequate, and finding the exact minimizer of TRS is unnecessary. Other
approaches that iteratively solve the TRS approximately include the dogleg
method [26, sect. 4.1] and DC-based algorithms [38].
All these algorithms require iterations of a computational routine, and the numberof iterations is often unpredictable and potentially large.
One exception to the iteration-based methods was proposed in 1989 by Gander,
Golub, and von Matt [10], which reduces TRS to a single quadratic eigenvalue problem, which they linearize to obtain a standard eigenvalue problem of size 2n. However,TRS VIA GEP

271

in that paper they report that their eigenvalue-based approach is slower and less accurate than Moré–Sorensen’s algorithm. This is perhaps why this approach appears tohave received less attention by TRS researchers than those mentioned above; anotherreason may be that the paper [10] revolves around (constrained) eigenvalue problems
rather than TRS.
The algorithm we advocate here, however, results in an extension of [10], whichturns out to be both efficient and accurate. Apparently, the slow speed and loss of
accuracy reported in [10] was largely due to the relatively unrefined eigenvalue solver
available those days. We demonstrate that the accuracy and speed are both dramatically improved by the highly developed eigensolvers available today, such as theshift-and-invert Arnoldi method [20]. Our MATLAB experiments demonstrate that
our algorithm has excellent accuracy and competitive speed compared with existingalgorithms in the standard case B = I with large-sparse A.
Our approach, based on a different derivation from [10], further extends the algorithm by allowing B 6= I and preserving symmetry. Moreover, the paper [10] does
not discuss how to deal with the “hard case,” which are TRS whose λ∗ is equal to thelargest eigenvalue of the pencil A + λB (which implies by (3) that g ⊥ N (A + λ∗ B),
where N (A+λ∗ B) denotes the null space of A+λ∗ B). In this paper we treat the hardcase in detail, and illustrate that our algorithm solves the hard case efficiently withhigh accuracy. We also demonstrate its effectiveness for computing an approximate
solution efficiently.
Note that the lack of iterations makes our algorithm easy to implement (its essence
is a few lines of MATLAB code), and the runtime predictable. Moreover, difficulties
associated with iterations—such as the need for safeguard techniques—simply disappear. Overall, this paper introduces a practical, general-purpose TRS algorithm,
useful for a wide range of problems, from dense-medium sized to large-sparse; we haveexperimented successfully with problems of size up to 107 .
The algorithms in [25, 30] are designed for B = I, and so are most publically available codes [30, 33]. Algorithms applicable to B 6= I are discussed in [4,sect. 7.5.6], [13, 15], but these algorithms are approximate and/or iterative. The traditional approach to deal with B 6= I is to reduce the problem to an equivalent one
1with B = I by a change of variables: defining p̃ := B 2 p, one can reduce (1) to an111
equivalent TRS with B = I, i.e., the minimization of (B − 2 g)> p̃ + 12 p̃> (B − 2 AB − 2 )p̃subject to kp̃k ≤ ∆. However, this reduction involves computing the matrix squareroot or the Cholesky factor of B and their inverse, which can be expensive and nu1
1
merically unstable. Furthermore, the conversion to B − 2 AB − 2 generally destroys
the problem structure: e.g., the inverse of an irreducible tridiagonal matrix is dense11(matrix-free algorithms such as [32] avoids forming B − 2 AB − 2 by working with the(sparse) Cholesky factor of B). Our algorithm works directly with A, B and thus
takes full advantage of the sparsity.
Besides choosing B to reflect the geometry of the problem such as B ≈ |A|,
another situation where an ellipsoidal norm arises is when a standard TRS with
B = I is solved via the Steihaug–Toint conjugate gradient-based algorithm [36, 39]
with preconditioning. We note that in [5] it is suggested that precautions are needed
when such change of the inner product is employed, as it can alter the quality of thesolution. This also indicates the importance of the choice of B. This paper makes noclaims on how to choose B, but rather takes it as given, and we focus on solving theTRS (1).
Our algorithm can be regarded as an extension of a related eigenvalue-based algo-

272

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

rithm in [19] for the point-ellipsoid distance problem, which minimizes kxk2 subject
to the constraint that x lies on the boundary of an ellipsoid. This problem has aconvex objective but a nonconvex constraint, and a generalized eigenproblem is derived in [19] that yields a global solution. We note that Rimon and Boyd [31] alsosuggests using the algorithm by Gander, Golub, and von Matt [10] mentioned above
for computing point-ellipsoid distances.
Along with [19], this work was initially inspired by the fact that some quadratic
optimization problems lead to polynomial equations, as follows. Consider TRS in
the simple case B = I and let A = V DV > be the eigenvalue decomposition witheigenvalues di , i = 1, . . . , n. We focus on the solution with kp∗ k = ∆, which is
generally the more difficult case than interior solutions kp∗ k < ∆ (see section 2.1).
Then by (3) we have p∗ = −V (D+λI)−1 V > g, and writing out the condition kp∗ k = ∆and defining ĝ = V > g, one sees that the solution can be obtained via a rational
equation in λ of the form
(6)

∆2 =

n
X

ĝj2
.
(dj + λ)2
j=1This equation has been presented in the literature [11, 25, 26], but is usually treatedas a “difficult” nonlinear equation that needs to be solved through an iterative process. Nonetheless, (6) is nothing else than a rational rootfinding problem, which canmathematicallybe reduced to a polynomial rootfinding problem by multiplying out
Qn
2
(d
+
λ)
,
which
can then be solved by a single eigenvalue problem via linearizaj
j=1
tion such as the companion matrix, without iterations (except those used within the
eigensolver). Although this “polynomialization” is usually not recommended due tonumerical instability (and we will not pursue it), this observation does suggest that amethod based on iterations is perhaps unnecessary.
This paper is organized as follows. In section 2 we detail the optimality conditions
and the KKT conditions. We then describe our main algorithm in section 3. Section 4
discusses how to deal with the hard case. We summarize the algorithm in pseudocode
in section 5. In section 6 we compare our algorithm with existing methods. Section 7
shows numerical experiments to illustrate and compare the performance.
Notation. Throughout this paper, R(X) denotes the range of a matrix X, andN (X) the null space. In addition, X  0 indicates X is a positive semidefinite matrix,
In is the n × n identity, and On , Om×n are the n × n and m × n zero matrices; we
simply write I, O if the dimensions are clear from the context. We always denote a
TRS solution by p∗ with associated Lagrange multiplier λ∗ . Computed approximants
wear a hat, so, for example, p̂∗ is a computed approximant to p∗ . We denote by u theunit roundoff, u ≈ 1.1 × 10−16 in double precision arithmetic.
2. Interior and boundary solutions. Equation (4) shows that a TRS solution(λ∗ , p∗ ) belongs to either (or both) of the following two cases, known as complementaryslackness: λ∗ = 0, or ∆ − kp∗ kB = 0. In view of this, we separate the TRS solution
into the following two types:(i) an interior solution with kp∗ kB < ∆,(ii) a boundary solution kp∗ kB = ∆.
We must have λ∗ = 0 in case (i). The “hard case” arises in the boundary case (ii), if
we further have det(A + λ∗ B) = 0; see section 4.
Below we first discuss how to deal with an interior solution, and then treat theboundary solutions (ii). Our algorithm by default proceeds by executing both steps

TRS VIA GEP

273

to find one or two candidates, and obtains a TRS solution by comparing the objective
value.
2.1. Interior solutions. An interior solution clearly satisfies λ∗ = 0 by (4).
Obtaining an interior solution is done by solving the linear system Ap0 = −g for p0 .
We then check if the constraint kp0 kB < ∆ is satisfied. If it is, then p0 is a candidate
for an interior TRS solution.
If further A  0, then p0 is indeed the TRS solution, If A has negative eigenvalues,by (5), p0 is not a TRS solution. Similarly, if A  0 has a zero eigenvalue, then solving
the linear system Ap0 = −g becomes nontrivial, but a TRS solution must exist on
the boundary (case (ii)). In practice, checking the positive definiteness of A may be
a costly operation (requiring O(n3 ) if A is dense), so instead, we simply attempt tosolve Ap0 = −g using a direct solver or preconditioned CG, and keep p0 as a candidate
if kp0 kB < ∆ and discard it if not.
2.2. Boundary solutions. The first three, (2)–(4), of the TRS optimality conditions in Theorem 1.1 represent the KKT conditions. The last, (5), specifies whichof the many (up to 2n, as we discuss in section 3.2) KKT multipliers corresponds tothe solution; we will show it is the largest real one.
By (3), for any KKT multiplier λ, unless the matrix A + λB is singular we canwrite p as a function of λ as
(7)p(λ) = −(A + λB)−1 g.

Since we focus on the case kpkB = ∆, plugging (7) into this equation we obtain(8)g > (A + λB)−1 B(A + λB)−1 g = ∆2 .

Now we regard A + λB as a matrix pencil with eigenvalues µ1 ≤ µ2 ≤ · · · ≤ µn , andlet W be the matrix that achieves the simultaneous diagonalization by congruence [12,sect. 8.7]
(9)

W > (A, B)W = (Λ, I),

where Λ = −diag(µ1 , . . . , µn ); note the minus sign as here we define µi as the eigenvalues of the pencil A + λB, not A − λB. Then (8) can be written as(10)(W > g)> (Λ + λI)−2 (W > g) − ∆2 = 0.

We write W = [w1 , . . . , wn ] and define(11)h(λ) :=

nX(wj> g)2− ∆2 .
2
(λ
−
µ
)
j
j=1

Equation (10) can be written as a zero-finding problem h(λ) = 0, which is a rational
equation with respect to λ, analogous to (6). This equation has 2n solutions in C
(counting multiplicities), and it is possible that they are all real.
We do not work directly with h(λ) = 0 as in (11), because this would require thecomputation of the eigenvector matrix W , which is quite expensive, and the resulting
method is not as accurate as the one we propose. Our aim is to construct a generalized
eigenvalue problem whose eigenvalues contain the values of λ satisfying (11), as we
describe next.

274

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

3. Eigenvalue-based algorithm for boundary TRS solutions. We now
come to the heart of the paper where we discuss finding a boundary TRS solution via
one generalized eigenvalue problem.
3.1. Key matrix pencils. The starting point is to introduce two matrix pencils
whose eigenvalues include the desired KKT multiplier λ∗ . Define the (2n+1)×(2n+1)
matrix pencil
 2

∆
0
g>
−B
A + λB (12)
M (λ) =  0
g A + λBOn
and the 2n × 2n matrix pencil−B
A + λB(13)M̃ (λ) =
A + λB>.
− gg∆2

The crucial facts that we prove are (i) the eigenvalues of these pencils include
the values of λ satisfying the KKT conditions, and (ii) we can obtain the solution(λ∗ , p∗ ) via the largest real eigenpair. We first examine the connection between the
eigenvalues and λ∗ , and then discuss how to obtain p∗ .
3.2. Finding λ∗ . We show that all the Lagrange multipliers at the KKT points(λ, p) on the boundary kpkB = ∆, including the desired λ∗ , are contained in the
eigenvalues of M (λ) and M̃ (λ).
Lemma 3.1. For any (λ, p) satisfying (A + λB)p = −g and kpkB = ∆, we havedet M (λ) = det M̃ (λ) = 0.
Proof. If det(A + λB) = 0 at λ with (A + λB)x = 0, then it follows from g ∈0n = 0, hence det M (λ) = det M̃ (λ) = 0.
R(A+λB) that M (λ) 0n+1
=
0
and
M̃
(λ)
x
x
We now deal with λ such that det(A + λB) 6= 0. Define p(λ) = −(A + λB)−1 g asin (7), and1.
X(λ) = p(λ) I
I
Then X(λ) is unimodular, i.e., det X(λ) ≡ 1, and we have
(14)

(15)

det M (λ) = det X(λ)T M (λ)X(λ)
 2∆ − p(λ)> Bp(λ) −p(λ)> B−Bp(λ)−B= det 0
A + λB
0
A + λB O= (−1)n det(A + λB)2 {∆2 − p(λ)> Bp(λ)}.

If (λ, p(λ)) satisfies kp(λ)kB = ∆ and det(A + λB) 6= 0, then recalling (8) we have
∆2 − p(λ)> Bp(λ) = 0, and so det M (λ) = 0. This proves det M (λ) = 0.
To prove det M̃ (λ) = 0, we define

1
T =

− ∆12 g >
In




In

275

TRS VIA GEP

and see that
(16)

 2
∆
T > M (λ)T = −B
A + λB 2∆
A + λB =− ∆12 gg >M̃ (λ).

It follows that det M (λ) = ∆2 det M̃ (λ), and, hence together with the above result,det M̃ (λ) = det M (λ) = 0.
Fortunately, both matrix pencils M (λ) and M̃ (λ) are regular matrix pencils, thatis, their determinants are nonzero for some λ and the number of eigenvalues is equal
to their size. To see this, observe that M̃ (∞) := OBn OBn is nonsingular, and thatM̃ (λ) is obtained from M (λ) by taking its Schur complement. Therefore, the numberof eigenvalues is finite, more precisely, 2n + 1 and 2n, respectively, and 2n of themmatch the 2n roots of the rational equation (11). Note that 2n, the size of the matrix
pencil M̃ (λ), is the smallest possible, since (11) is a rational equation that can be
reduced to a degree-2n polynomial equation with 2n solutions. Each real eigenvalue
λ of M (λ) corresponds to a KKT point p = (A + λB)−1 g with KKT multiplier λ.
Lemma 3.1 shows that the TRS solution on the boundary satisfies det M (λ) = 0
and det M̃ (λ) = 0, both representing a generalized eigenvalue problem. The eigenvalues λ contain the Lagrange multipliers at the KKT points, so the multiplier for
the TRS solution must be one of the 2n finite eigenvalues (note that M (λ) has one
eigenvalue at infinity, which is not the one of interest).
We next show that among these 2n finite eigenvalues, λ∗ corresponds to the largest
real one lying in [µn , ∞).
Theorem 3.2. For a TRS solution (λ∗ , p∗ ) on the boundary kp∗ kB = ∆ satisfying (2)–(5), the multiplier λ∗ is equal to the largest real eigenvalue of M (λ) (excluding
λ = ∞) and M̃ (λ), and λ∗ ∈ [µn , ∞), where µn is the largest eigenvalue of the pencil
A + λB.
Proof. The fact λ∗ = λmax ∈ [µn , ∞) has been shown in [7, 19] for special cases:[7] for B = I and [19] for A = I; see also [24]. By a change of variables we can extend
these results to the TRS (1).
Alternatively, we can directly obtain the result as follows. First, from A+λ∗ B  0
in (5) we must have λ∗ ≥ µn .
Recall from (11) that kp(λ∗ )kB = ∆ is equivalent to h(λ∗ ) = 0, unless λ∗ = µn .
To prove that λ∗ is the largest real eigenvalue of M (λ), we consider the two cases
λ∗ > µn and λ∗ = µn separately. First, when λ∗ > µn , we have h(λ∗ ) = 0 anddet(A + λ∗ B) > 0 in (15). Note from (11) that h(λ) is strictly decreasing on (µn , ∞).
Hence, M̃ (λ) has exactly one real eigenvalue larger than µn , which must be λ∗ .
Next, when λ∗ = µn (this is the “hard case”), from (3) we see that g ⊥ N (A +λ∗ B), and so h(λ) has no pole at λ = µn , and strictly decreasing on (µn−` , ∞), where`(≥ 1) is the integer such that µn−` < µn−`+1 = · · · = µn . Hence, h(λ∗ ) is formallyPn−` (wj> g)2Pn−` (wj> g)2
2
2= j=1 (µn −µdefined as h(λ∗ ) = j=1 (µn −µ2 − ∆ . Now h(λ∗ ) + ∆2 is equalj)j)
to the smallest kpk2B such that (A + λB)p = −g (we prove this in (24)). Therefore,
if h(λ∗ ) > 0, then no solution exists with (3) and kp∗ kB = ∆. Hence, h(λ∗ ) ≤ 0,
and since h(λ) is strictly decreasing on (µn−` , ∞), there is no solution for h(λ) = 0
with λ > µn , that is, λ∗ is the largest real eigenvalue of M (λ) and M̃ (λ) also in thiscase.
3.3. Obtaining p∗ . We next show that, generically, the TRS solution p∗ can beobtained from the eigenvector of M (λ) and M̃ (λ) corresponding to λ∗ .

276

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

Theorem 3.3. The eigenvectors of M (λ) and M̃ (λ) for the largest real eigenvalue
λ = λ∗ correspond one-to-one as 1 >  
− ∆2 g y2y1
 = 0.
y1(17)M̃ (λ∗ )= 0 ⇔ M (λ∗ ) y2y2Moreover, if g > y2 6= 0, then a TRS solution p∗ can be obtained by∆2p∗ = − > y1 .
g y2

(18)

Proof. We prove that (17) holds for every
finite eigenvalue λ of M (λ). By thecongruence relation (16), we see that if yy12 is an eigenvector of M̃ (λ), then
 1 >  
− ∆2 g y2y = 0.
y1(19)M̃ (λ) 1 = 0 ⇒ M (λ) 
y2
y2
h t i
For the converse, suppose that M (λ) yy1 = 0. Then since T in (16) is nonsingular
2
with


1 >
1
∆2 g
,
In
T −1 = In
we have  2t∆ (t + ∆12 g >y2 ) = 0,yT −1 y1  = T −> M̃ (λ) 1
M̃ (λ)
y2
y2
h
i
 
 
t
thus M̃ (λ) yy12 = 0 and t = − ∆12 g > y2 . Now noting that M (λ) yy1 = 0 with yy12 6= 0
 
 2
t
−> ∆


M (λ) y1 = T
y2



2

can hold only for λ = ∞, we obtain the first statement (17).
We next prove (18). From M̃ (λ∗ ) yy12 = 0, that is,

 

 
−B
A
y1
O B y1
>
(20)
= −λ∗
,
y2
B O y2A − gg∆2>

the lower block gives (A+λ∗ B)y1 = g(g∆2y2 ) , which is a scalar multiple of g. Therefore,in view of (A + λ∗ B)p∗ = −g from (3), we obtain (18) as required, provided thatg > y2 6= 0.
Practical extraction of p∗ . In practice, computing the solution p∗ from thenormalization (18) may introduce unnecessary numerical errors, and we choose thesimpler normalization(21)p∗ = −sign(g > y2 )∆y1,
ky1 kB

which is on the boundary to working precision: kp̂∗ kB = ∆+O(u). The normalization(21) is obtained by scaling the vector y1 (which is parallel to the solution p∗ ) to haveB-norm ∆.
When g > y2 = 0, finding the solution is not straightforward; this is described in
section 4.

277

TRS VIA GEP

3.4. Observations for a practical implementation.
3.4.1. The rightmost eigenvalue is real. We have seen that λ∗ at the TRS
solution is equal to the largest real eigenvalue of M (λ) and M̃ (λ). For the eigensolver
it helps to further know that λ∗ is indeed the rightmost eigenvalue, that is, there isno nonreal eigenvalue that lies to the right of λ∗ .
Proposition 3.4. The rightmost eigenvalues of M̃ (λ) and M (λ) (excluding ∞)are both real and equal to λ∗ .
Proof. It suffices to show that the rightmost eigenvalue of M̃ (λ) is real. Supposethat λ̃ = α + βi where α, β ∈ R, α > µn , and β > 0 is a nonreal eigenvalue. Thenα − βi must also be an eigenvalue. Now if α > λ∗ ≥ µn , then recalling (11) we havePn (wj> g)2
1
2
is strictly negativeh(λ) = j=1 (λ−µ2 −∆ , and since the imaginary part ofj)(λ̃−µ )2jfor all j, we conclude that the imaginary part of h(λ̃) must also be strictly negative.
Hence, h(λ̃) cannot be 0, so λ̃ is not an eigenvalue of M̃ (λ).
The above discussion proves the statement except when α = µn = λ∗ , that is, toexclude λ̃ = λ∗ + βi in the hard case. In this case, we note that the terms wj> g mustbe zero for all indices j such that λj = λ∗ . For the remaining indices, wj> g cannot be
1) < 0.
all zeros (since W is nonsingular), and again we have Im( (λ̃−µ)2jWe have shown that the TRS solution can be obtained from the rightmost (whichis the largest real) eigenpair of M (λ) or M̃ (λ), which is generally much cheaper
to compute than the whole eigenvalues via the standard QR or QZ algorithms; forexample, the Arnoldi algorithm provides an effective means for computing extremal
eigenvalues of large-sparse matrices.
3.4.2. Which matrix pencil to use, M (λ) or M̃ (λ)?. The generalized eigenvalue problem M (λ) has one additional eigenvalue at ∞, but the matrices involved are
explicitly sparse. On the other hand, M̃ is one size smaller with no eigenvalue at infinity, but contains gg > , which is rank-one but dense. Thus the choice between M andM̃ should be made based on the properties of the eigensolver available. Some eigensolvers, such as MATLAB’s eigs, which is based on (shift-and-invert) Arnoldi [20],allow the user to provide just a routine that multiplies the matrices to a vector; inthis case the rank-one structure can be exploited. For this reason we use M̃ in our
MATLAB experiments.
3.4.3. Comparison with the algorithm by Gander, Golub, and vonMatt. The algorithm by Gander, Golub, and von Matt [10] considers the case B = I
and finds the largest eigenvalue of the 2n × 2n standard but nonsymmetric eigenvalue
problem

 
 
A
−I y2
y
(22)
= −λ∗ 2 .
gg >
y
y1
− ∆2
A
1
This is a nonsymmetric matrix, with a dense bottom-left block. We can obtain (22)
when B = I by the equivalence transformation of right-multiplying IOn IOn to M̃ .
Therefore, we arrived at a different derivation of (22) and generalized it to B 6=
I. In [10] solving TRS via (22) is not recommended over Moré–Sorensen’s secularequation approach [22], observing the inefficiency and inaccuracy with (22) in theirexperiments.
As we shall see in our experiments, the poor performance seems to have been
largely due to the relatively undeveloped eigenvalue solvers available at the time:

278

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

with state-of-the-art
  eigensolvers the algorithm is both fast and accurate. Between
solving M̃ (λ∗ ) yy12 = 0 and (22), the speed comparison depends on the architecture
etc: in our experiments (22) was often slightly faster, so a reasonable option is to turnto (22) when B = I. In any case, our experiments illustrate that by solving (20),
our algorithm has efficiency comparable to the state-of-the-art algorithms. We do notclaim to have a faster algorithm than [10] when B = I (when employing the same
eigensolver), but instead we observe its practicality and extend it to B 6= I, and treat
the hard case.
4. Dealing with the “hard case”. We saw in Theorem 3.3 that a TRS solution(λ∗ , p∗ ) can be obtained via the rightmost eigenpair of M (λ) provided that g > y2 6= 0.
We now discuss the case where this does not hold, i.e., g > y2 = 0. In this case λ∗is still obtained as the rightmost eigenvalue of M̃ (λ) by Theorem 3.2, but finding p∗
requires more work.
We show that the condition g > y2 = 0 can happen only if λ∗ = µn .
Proposition 4.1. In the setting of Theorem 3.3, g > y2 = 0 only if λ∗ = µn .
Proof. Recall from (20) that g > y2 = 0 implies (A + λ∗ B)y1 = 0. From the first
block of (20), we obtain(23)−By1 + (A + λ∗ B)y2 = 0,

and hence (A + λ∗ B)B −1 (A + λ∗ B)y2 = 0. Using the diagonalization W > (A, B)W =(Λ, I), this yields W > (λ∗ I+Λ)2 W y2 = 0, from which we obtain W > (λ∗ I+Λ)W y2 = 0.
This is equivalent to (A + λ∗ B)y2 = 0, which, by (23), also implies y1 = 0. Hence, y2
is a nonzero eigenvector of A + λB corresponding to the eigenvalue λ∗ = µn .
Note that the proof also shows that if g > y2 = 0, then y1 , which we usually
use to extract the solution p∗ , is zero. Also note from (3) that λ∗ = µn impliesg ⊥ N (A + λ∗ B). Hence, g > y2 = 0 implies λ∗ = µn , which in turn implies g ⊥
N (A + λ∗ B). This is the so-called “hard case,” a difficulty that is known to arise in
the standard B = I case, and it is of course also present when B  0 is a general
positive definite matrix. We restate the definition for general B  0.
Definition 4.2. A TRS (1) is said to be in the “hard case” if λ∗ = µn , thelargest eigenvalue of the pencil A + λB.
We repeat that λ∗ = µn implies g ⊥ N (A + λ∗ B), which some papers take to be
the definition.
It is worth noting that the last condition pertains to the orthogonality between gand the eigenvectors corresponding to the largest eigenvalue of A+λB in the standard
inner product, not B-orthogonal wi> Bg = 0 as one might expect since the TRS (1)employs the B-norm. We also note that a TRS solution lies on the boundary in thehard case, since we have either λ∗ > 0 (trivially a boundary solution), or A  0 with
A having null vectors, so again a TRS solution is on the boundary.
Although mathematically the hard case represents only a set of TRS instancesof measure zero, it can happen for matrices with special structures, and numericallythere are “nearly hard cases,” in which g > y2 ≈ 0 and hence λ∗ ≈ µn . These can beequally challenging. Indeed a number of studies such as [25, 32] have focused on thehard case with B = I.
In our approach, the reason the hard case is difficult is that recalling (18), theeigenvectors of M or M̃ for λ∗ do not give us the TRS solution. Indeed the linear
system (A + λ∗ B)x = −g has infinitely many solutions. The challenge is, therefore,

279

TRS VIA GEP

to find the solution p∗ such that (A + λ∗ B)p∗ = −g and kp∗ kB = ∆. We manage to
do so by modifying an approach in [8] to adapt to a general positive definite B.
4.1. Solution for the hard case. We form a nonsingular linear system using
N (A + λ∗ B) such that we can obtain p∗ from its solution.
Theorem 4.3. Suppose the TRS problem (1) belongs to the “hard case” and(λ∗ , p∗ ) satisfies (3)–(5) and ||p∗ ||B = ∆ with λ∗ = µn . Let d = dim(N (A + λ∗ B))
and V := [v1 , . . . , vd ] be a basis of N (A+λ∗ B) that is B-orthogonal, i.e., V > BV = Id .
For an arbitrary α > 0, defineH :=A + λ∗ B + αdXBvi vi> B .

i=1

Then H is positive definite. Moreover, q := −H −1 g is the minimum-norm solution
to the linear system (A + λ∗ B)p = −g in the B-norm, that is,
q = argminp {kpkB |(A + λ∗ B)p = −g}.

(24)

Furthermore, for any nonzero v ∈ N (A + λ∗ B) there exists a scalar η ∈ R such thatp∗ = q + ηv is a TRS solution.
Proof. First, we prove that H  0. Let W be the matrix that simultaneously diagonalize A, B as in (9). Then W > (A+λ∗ B)W = diag(λ∗ −µ1 , . . . , λ∗ −µn−d , 0, . . . , 0)
is a diagonal matrix and W > BWP= I.
d
We claim that defining G = i=1 Bvi (Bvi )> = BV V > B, we haveOn−d
W > GW =.
Id
To see this, writing W = [W1 W2 ] we have (A + λ∗ B)W2 = (A + λ∗ B)V = On×d andW2> BW2 = V > BV = Id . Since V, W2 are of the same size, these two equalities imply
that V, W2 are both B-orthonormal matrices that span the same subspace, namely
the space spanned by the eigenvectors of A + λB corresponding to the eigenvalue λ∗ .
Hence there exists an orthogonal matrix Q ∈ Rd×d such that V = W2 Q. Therefore,
W > (BV V > B)W = [W1 W2 ]> BW2 Q(W2 Q)> B[W1 W2 ]
= [W1 W2 ]> BW2 W2> B[W1 W2 ]O(n−d)×d On−dOd×(n−d) Id =
=
Id


Id

,

where we have used the fact W > BW = [W1 W2 ]> B[W1 W2 ] = In . Therefore,>

W HW =>W AW + λ∗ In + α

d
X

W

>

Bvi vi> BW

i=1= diag(λ∗ − µ1 , . . . , λ∗ − µn−d , α, . . . , α)is positive definite, so by Sylvester’s law of inertia it follows that H is also positive
definite.
We next show that q := −H −1 g is a solution to the singular linear system (A +λ∗ B)p = −g, which necessarily has infinitely many solutions. To see that (A+λ∗ B)q =280

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA−g, define Λ̃ = −diag(µ1 , . . . , µn−d ) (recall Λ in (9)) and note that −g = (A + λ∗ B)p
implies (for the remainder of the proof, I = In−d and O = Od )−> Λ̃ + λ∗ I
−g = WW −1 p,O
and also thatq := −H −1 g = −W

(25)

−1


Λ̃ + λ∗ I
αIW >gfrom which we obtain(A + λ∗ B)q = −W −>
= −W

−>


Λ̃ + λ∗ I


O


I



W −1 W

>

αI

I

−>

−1


Λ̃ + λ∗ I


>

W g=W
W W
O


Λ̃ + λ∗ I= W −>
W −1 p = −g.
O

W >g

−> Λ̃ + λ∗ IOOW −1 p

To prove that q = −H −1 g is the minimum B-norm solution to (A + λ∗ B)p = −g,note that we can write a general solution to (A + λ∗ B)p = −g as p = q + v wherev ∈ N (A + λ∗ B). We shall show that kqkB ≤ kq + vkB for any such v. We havekq + vk2B = (q + v)> B(q + v) = q > Bq + 2v > Bq + v > Bv.
Since B is positive definite the third term is nonnegative, it suffices to show thatv > Bq = 0 for all v ∈ N (A + λ∗ B). To see this, using q = −H −1 g = H −1 (A + λ∗ B)q
we obtain Bq = BH −1 (A + λ∗ B)q, and note that
BH −1 (A + λ∗ B) = (A + λ∗ B)H −1 B,
which we can also verify using the decomposition with respect to W (essentiallybecause diagonal matrices commute):
BH −1 (A + λ∗ B)
=W

−>

=W

−>

W

−1

W


I


O

−1


Λ̃ + λ∗ I

W

αId

>

!
W

−>




Λ̃ + λ∗ I
O

W

−1



W −1



Λ̃ + λ∗ I= W −>


O

W

−1


W

−1


Λ̃ + λ∗ I
αId

!
W

>

W −> W −1= (A + λ∗ B)H −1 B.
Therefore, we obtain for all v ∈ N (A + λ∗ B)v > Bq = v > (BH −1 (A+λ∗ B)q) = v > ((A+λ∗ B)H −1 Bq) = (v > (A+λ∗ B))H −1 Bq = 0,
as required. Note from (25) that kqk2B =of Theorem 3.2.

(wj> g)2
j=1 (µn −µj )2 , which we used in the proof

Pn−1

TRS VIA GEP

281

It remains to establish the final statement. Since in the hard case the solution
is on the boundary with Lagrange multiplier λ∗ equal to µn , we are able to obtain avector q + ηv on the boundary by finding the scalar η by solving a scalar quadratic
equation kq + ηvk2B = ∆2 , to obtain a global solution to the TRS satisfying (2)–(5).
In the last paragraph of the proof, it may appear that kqkB > ∆ is possible, inwhich case there is no η such that q+ηv is on the boundary. However, by Theorem 1.1,
there must be a boundary solution in the hard case. In other words, if even theminimum-norm solution had norm kqkB > ∆, then this implies the situation was not
the hard case, and the process in the previous section solves the TRS.
The proof above shows that in the hard case the TRS solution p∗ is generally notunique; the goal is to find one solution, which we denote simply by p∗ .
Theorem 4.3 shows that p∗ can be computed by finding the null vectors of A+λ∗ B,
forming H, and solving the nonsingular linear system Hx + g = 0. Finding η is thenan easyP
task of solving a scalar quadratic equation. Note that due to the low-rankdterm α i=1 Bvi vi> B, H can be dense even when A, B are sparse. Nonetheless, thelinear system Hx + g = 0 can be solved efficiently by the CG algorithm (employing anappropriate preconditioner if available), which (as with any Krylov-type algorithm)
only requires a routine for computing matrix-vector multiplications with respect tothe coefficient matrix.
4.2. Detecting the hard case in our algorithm. As noted after Proposition 4.1, the hard case is indicated by y1 = 0, that is, the vector y1 that we usually obtain the solution from becomes zero; indeed, otherwise we obtain the solutionby (21).
This suggests that we can detect the hard case by looking at the first half elements
y1 of the computed eigenvector y = yy12 , which we take to have unit norm kyk = 1:
the hard case is when y1 = 0. In practice, due to roundoff errors the computed
vector ŷ1 (recall that quantities wearing a hat represent computed approximations)has nonzero but small elements in the (near) hard case. To distinguish the hard andgeneric cases, we set a threshold τ < 1 such that TRS with ky1 k < τ is treated as
hard case. To set an appropriate value for τ , note that if the top part ŷ1 is small innorm, denoting it by , we have
# 
  "
−B
A + λ̂∗ B


.
(26)
M̃ (λ̂∗ )
=
gg >
ŷ
ŷ2
2
A + λ̂∗ B− ∆2
For this to be numerically zero, by the first block row we have
(A + λ̂∗ B)ŷ2 = B.
Since the right-hand side is O(), it follows that (λ̂∗ , ŷ2 ) can be regarded as an approximate eigenpair for (A, B).
To choose an appropriate threshold we analyze the accuracy
 of the computed ŷ1
as an approximation to p∗ . Let the computed eigenvector ŷ = ŷŷ12 be normalized tohave unit norm. Assuming the matrix of eigenvectors is well conditioned, the accuracyof ŷ is known to be O( residualgap ) [37, Chap. 5]. Here the residual is kM (λ̂∗ )ŷk, whichis generally O(u) with a numerically stable algorithm,1 where u is the unit roundoff,
1 To avoid unnecessary cluttering we assume kAk, kBk are O(1). This causes no loss of generality
as we can scale the matrices without changing the TRS solution: A ← c1 A by taking g ← c1 g, andB ← c2 B by taking ∆ ← c2 ∆.

282

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

and gap is the smallest distance between λ∗ and the rest of the eigenvalues of M (λ).
Moreover, the loss of accuracy in extracting a vector of norm ky1 k as a part of aunit-norm vector is a factor O( ky11 k ). Since the solution p∗ is obtained by normalizing
y1 as in (21), overall the accuracy in p∗ is estimated to be O( ky1 kugap ).
On the other hand, if we treat the problem as the hard case, recalling Theorem 4.3,
we need to compute the null vectors N (A+ λ̂∗ B). Numerically, these are the vectors v
for which k(A+ λ̂∗ B)vk is negligible. With the tolerance τ for detecting the hard case,recalling (26), we expect the vectors we consider will have k(A + λ̂∗ B)vk = O(ky1 k),suggesting this entails an error of size O(ky1 k).
We suggest choosing the threshold τ for ky1 k based on which is likely to give theu= τ , that is,more accurate solution, according to the above discussion: τ gap
r
(27)

τ=

u
.
gap

p
In double precision arithmetic u ≈ 10−16 , and we choose τ to be about 10−8 1/gap.
We can estimate the gap by computing two rightmost eigenvalues of M (λ) instead ofone. Alternatively and more efficiently, since a rough estimate for gap would suffice,we can use the inverse power method for M (λ) = M0 +λM1 by starting with a random
vector w that is M1 -orthogonal to the eigenvector v corresponding to λ∗ , and taking
a few iterations of power method with the matrix (M0 + λ∗ M1 )−1 M1 , executed bysolving the linear system (M0 + λ∗ M1 )wnew = M1 w for wnew . Here we solve thelinear system only crudely using an iterative solver such as MINRES [28].
Below is the process to deal with the hard case TRS in pseudocode.
Algorithm 1 Algorithm for hard-case TRS.
1: Compute the eigenvectors of (A, B) corresponding to λ∗ (i.e., the null vectors of
A + λ∗ B).
2: Solve Hq + g = 0 for q by the CG method.
3: Take an eigenvector v ∈ N (A + λ∗ B) computed above, and find η ∈ R such thatkq + ηvkB = ∆ by a quadratic equation.
4: Return q + ηv as a candidate for the global TRS solution.
5. Pseudocode. We summarize our TRS algorithm in pseudocode (Algorithm 2).
The core part of Algorithm 2 is (28). We note that some steps can be skipped if additional information is available. For example, we can conclude that p0 is the (unique)optimal TRS solution if (i) it is feasible kp0 kB ≤ ∆, and (ii) A  0. Therefore, when
the positive definiteness of A is known or easily verifiable, if kp0 kB ≤ ∆ after step 1
in Algorithm 2, then we can dismiss the remaining steps in Algorithm 2 and take p0
as the TRS solution.
Note that the CG algorithm is originally designed for positive definite linear systems, and if CG does not converge, then this implies that A is numerically indefinite.
However, in practice CG often converges even when A is indefinite [27], so we cannot
conclude that A is positive definite just because the linear system Ap0 = −g was
solved by CG.
Another situation worth mentioning is when A is known (or easily verifiable) to
be indefinite with one or more negative eigenvalues. In this case a TRS solution mustlie on the boundary, and so there is no point in executing step 1; we directly proceed
to step 2.

283

TRS VIA GEP

Algorithm 2 Solve the TRS (1).
1: (Consider the case λ∗ = 0.) Solve Ap0 = −g by the CG algorithm, and keep p0 if
it satisfies kp0 kB < ∆.
y 
2: Compute the rightmost eigenvalue λ∗ of M̃ (λ) and an eigenvector y12 such that

−B
A

(28)

 

y1
0
=
−λ
∗
y
B
− gg
2
2
∆
A

>

B
0

 
y1
.
y2

3: Estimate
qthe gap between λ∗ and the other eigenvalues via a power method.

ky1 k ≤

If
,
then
treat
as
hard
case:
run
Algorithm
1
toobtainp.
Otherwise,1
gapu

obtain p1 by p1 := −sign(g > y2 )∆ kyy11kB .
4: The solution p∗ is either p1 or p0 (if it exists), whichever gives the smallest objective value.

6. Comparison with existing methods. Here we compare our Algorithm 2
with previous methods and argue that ours has attractive properties in terms ofefficiency and simplicity. Moreover, numerical experiments illustrate the excellent
accuracy as we see in section 7.
6.1. Efficiency. The computational cost of our algorithm is O(n3 ) flops for
dense A, B. When A, B are large and sparse, it is essentially the cost of an Arnolditype method for computing one rightmost eigenpair, which is typically in the order ofa constant times the cost of a matrix-vector product, or a shifted-and-inverted linear
system (M0 + σM1 )x = b for some σ ∈ C; the cost also depends on the separationof eigenvalues etc. This is in the same ballpark as the cost of the existing algorithms[13, 25, 35] when B = I, both in the dense and large-sparse cases.
Let us compare the cost in more detail. Recall that most conventional algorithms
for an accurate TRS solution involve n × n linear systems or eigenvalue problems
iteratively, whereas ours solves a double-sized 2n × 2n eigenproblem once. The advantage of a noniterative approach becomes significant especially when the eigenproblem
scales mildly with n. For example, in the dense case where eigensolvers (and linear
systems) require O(n3 ) cost, our approach is roughly comparable to a conventionalalgorithm that iterates eight times. If the eigensolver needs O(np ) cost with p < 3,then this number reduces to 2p .
We note that (when B = I) it is possibly less than twice as expensive to solvethe 2n × 2n eigenvalue problem M̃ rather than the size-(n + 1) eigenproblem with
 >respect to the matrix N (α) = αg gA , which is used in [35, 32, 33]. This is because
the dominant cost in the Arnoldi iteration  is in matrix-vector multiplication, andmultiplying the matrix (22) to a vector xx12 can be done by computing A[x1 x2 ] anda vector-vector multiplication g > x1 . Now, computing A[x1 , x2 ] is usually faster than
two independent matrix-vector multiplications Ax1 , Ax2 due to the use of a higher
level BLAS routine [12].
6.2. Ease of implementation. As mentioned before, the main feature of ourapproach is that the TRS can be solved essentially by one generalized eigenvalue
problem.
Besides the aesthetic advantage of directly giving a solution without iterations,another advantage of our approach is its ease of implementation. For example, in284

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

MATLAB we can execute (28), the main part of Algorithm 2, in just four lines (hereM̃ (λ) = M0 + λM1 ):M0 = [-B A;A -g*g’/Delta^2];
M1 = [zeros(n) B;B zeros(n)];lam = max(eig(M0,-M1));
p = -(A+lam*B)\g.
This is strikingly simple when compared with those of the existing algorithms, suchas the codes of [9, 14, 15, 33]. When A, B are large-scale and sparse, it is advised toreplace eig with eigs and M̃ with M . Specifically, in the large-scale case, after thesecond line the code should be[v,lam] = eigs(@(x)M0x(x),2*n,-M1,1,’lr’);p = v(n+1:end);p = p/sqrt(p’*B*p)*Delta.
The last line is a normalization to force kpkB = ∆, because the computed eigenvector
is normalized so that kvk2 = 1. Here M0x(x) is a function handle that left-multiplies
M0 to an input vector x. This saves memory over storing the matrix M0 because thisway we essentially need only store the matrices A and B, along with the vectors g, v.
This way, despite their doubled size, storing M, M̃ requires no more memory thanA, B, and g.
6.3. TRS as a subproblem. We argue that Algorithm 2 is the first algorithm
to solve the TRS accurately and is suited to the large-scale sparse case with B 6= I
without requiring a change of variables or outer iterations. However, its performancefor efficiently solving the overall nonlinear optimization problem using TRS as a subproblem is not easily predictable, in which solving the TRS exactly is not necessary,and it suffices if the TRS solution results in sufficient reduction in the original nonlinear function f .
To gain an idea of the performance in such cases, we present experiments where
we obtain approximate TRS solutions by stopping the Arnoldi iteration early when
computing the eigenpair (28).
7. Numerical experiments. We now turn to experiments to examine the performance of the proposed methods for a variety of TRS instances. All experiments
were carried out in MATLAB 2013A on a Blade server machine with Xeon CPU and
64 GB memory.
In the figures below, Algorithm 2 is shown as GEP (standing for generalized
eigenvalue problem). We compare GEP with three other codes that are publicallyavailable, implementing existing algorithms: (i) the code by Fortin–Wolkowicz [9]which is based on the algorithm by Rendl and Wolkowicz [30] (shown as FRW),(ii) Rojas, Santos, and Sorensen [33], based on [35] (shown as RSS), and (iii) thegalahad gltr in the Galahad library [14] (shown as GLTR), based on the algorithm
by Gould et al. [13]; this code is written in Fortran, and here we used its MATLAB
interface. We used the default parameters for FRW and RSS, and for GLTR, weused the options extra vectors=100 and stop relative=1e-15 (suggested by Nick
Gould, chosen to obtain high-accuracy solutions).
We examine the performance in the following cases: (a) B = I and A is sparse,(b) B = I and A is dense, (c) B 6= I and A, B are sparse, and (d) the hard case with
B = I.
In each set of experiments, we ran the codes 20 times and reported the averageruntime and accuracy. To measure the accuracy, we have computed the relative

285

TRS VIA GEP

objective function difference as follows:f (p̂∗ ) − f (p̂best ).
|f (p̂best )|(29)

Here p̂∗ denotes the computed solution of each method and p̂best is the solution withthe smallest objective value among the four algorithms. The reason the accuracy
measure is usually positive in the plots below (exceptions include Figure 1 (right) for
GLTR) is that the algorithm that achieves p̂best varies from problem to problem, andwe report the average of 20 runs. Throughout, unless otherwise specified, we take2∆ = 1 and the vector g is randomly generated by randn(n,1). To ensure the computed solutions are always feasible, after each algorithm we applied the normalizationp̂ := ∆p̂/kp̂kB if the computed p̂ violated the constraint kp̂kB ≤ ∆.
When B = I and A is sparse. We first set B = I and let A be large-sparse
matrices, varying the matrix size n from 103 (moderate) to 107 (large).
In Figure 1 we show the results with A = L − 5I, where L is the discrete twodimensional (2D) Laplacian (Figure 1); this follows the test matrices in [33]. Forn > 106 we ran only GEP and GLTR, which complete the computation in a reasonableamount of time.
100

104

10-5
accuracytime(s)

102

FRW
RSS
GEP
GLTR

100

FRW
RSS
GEP
GLTR

10-10

10-2

10-15
10-4
103

104

105
matrix size

106

107

103

104

105
matrix size

106

107

Fig. 1. Runtime (left) and accuracy (right) for A = L − 5I and B = I.

We also experiment with a random symmetric sparse matrix (Figure 2), generatedby(30)

A = sprandsym(n, density),

where we took density=1e-2.
The left plots of Figures 1 and 2 show the runtime in seconds. We see that while
GLTR was the fastest in these examples—our algorithm was notably faster than FRW
and RSS, the codes natively programmed in MATLAB. Figure 1 indicates that GEP
scales like O(n) for this problem. Recalling the discussion in section 6.1, this isfavorable for our noniterative algorithm, because the iterative algorithms often needto solve more than ten eigenproblems of half the size.
The right plots of Figures 1 and 2, which show the difference in the objective
values, illustrate that our algorithm GEP obtained solutions within about 10−15 of
2 The choice of ∆ is usually made adaptively during the overall trust-region method [4, sect.
6.1]. It is worth noting that the volume of the unit ball in Rn shrinks factorially with n [18], so
a geometrically appropriate choice may be far from ∆ = 1. Our experiments suggest that theperformance of GEP depends only mildly on the value of ∆.

286

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA
100
FRW
RSS
GEP
GLTR

10-5
accuracytime(s)

102

100

FRW
RSS
GEP
GLTR

10-10

10-2

10-15
10

-4

102

103

104

102

105

103

104

105

matrix size

matrix size

Fig. 2. Runtime (left) and accuracy (right) for A = sprandsym(n, 1e − 2) and B = I.

the optimal for every problem, which are accurate enough to be regarded as exact
solutions in double precision arithmetic. These results suggest that with the eigensolvers available today, when A, B are sparse, our algorithm is applicable and effectivefor very large problems.
When B = I and A is dense. We next examine the dense case. We generate
A by forming n random real numbers µi via randn(n,1), then generating a random
orthogonal matrix Q, formed by orth(randn(n)) and setting A = Q> diag(µi )Q.
Clearly, we are limited to much smaller matrix size n than in the previous sparse
case; here we take n ≤ 104 .
The results are shown in Figure 3. The accuracy behaved much the same as in thesparse case. For the speed, the difference is more benign than in the sparse case. We
can explain this qualitatively as follows: in the dense case all the algorithms requireO(n3 ) operations, and, recalling the discussion in section 6, our algorithm is expected
to be fast especially when the generalized eigenproblem can be solved efficiently. In
the sparse case its cost is often much less than O(n3 ) as we saw above, and this is afavorable situation for our algorithm. In the dense case, for which all algorithms are
expected to scale like O(n3 ), our algorithm still has efficiency comparable with otherapproaches. GLTR appears to scale slightly worse than the other algorithms here.
Another alternative for dense problems is GALAHAD’s galahad trs, which implements the algorithm in [15], although it did not perform better in the experiments
here.
100
FRW
RSS
GEP
GLTR

10

-2

10

-4

FRW
RSS
GEP
GLTR

10-5

100

accuracytime(s)

102

10-10

10-15
101

102

103
matrix size

104

101

102

103
matrix size

Fig. 3. Dense A, B = I. Runtime (left) and accuracy (right).

104

287

TRS VIA GEPWhen B  0 and A is sparse. We now examine problems with B 6= I. We let A
be sparse as in and B be a tridiagonal matrix, defined as A = sprandsym(n,0.01)
as in (30) and B = tridiag(1,3,1). g is a random vector as before.
Since FRW and RSS are not directly applicable when B 6= I, to invoke these algorithms we first compute the Cholesky factorization B = LL> , then form L−1 AL−>and apply the codes to A ← L−1 AL−> , g ← L−1 g, and B ← I. RSS is a matrix-freealgorithm that allows the input for the matrix A to be a routine for matrix-vector
multiplication, so we can use RSS without forming the (dense) matrix L−1 AL−> bysolving a linear system with L, L> for each matrix-vector multiply, thus RSS allows
us to work only with sparse matrices; in the above example, L here has bandwidth
2. The code FRW does not allow such input, which makes the code too costly inmemory and hence inapplicable for large n. GEP and GLTR do not require modified
treatment when B 6= I; they handle the B 6= I case exactly the same way as when
B = I.
Figure 4 shows the results. As before, our algorithm outperforms FRW and
RSS. GLTR is the most efficient here, although the difference shrinks with n as GEP
appears to scale better.
100

100

10

-2

10

-4

FRW
RSS
GEP
GLTR

10-5
accuracytime(s)

102

FRW
RSS
GEP
GLTR

10-10

10-15
102

103

104
matrix size

105

102

103

104

105

matrix size

Fig. 4. Sparse A, B, B 6= I. Runtime (left) and accuracy (right) for A = sprandsym, B
tridiagonal.

Hard case. We now turn to the hard case. The standard process is then insufficient, and algorithms usually employ a remedy for the hard case, such as Algorithm 1for GEP. An exception is RSS, in which the hard case is treated within the main
iteration of the algorithm.
We let B = I and form a tridiagonal A, with 2 on the diagonal and random N (0, 1)
elements on the off-diagonals; the Hessian A is tridiagonal if the objective function in
the nonlinear optimization problem depends only on adjacent variables xi−1 , xi , xi+1 .
To generate a “hard case,” we first set g to be a random vector and calculate thelargest eigenvalue λ∗ of the pencil A + λB with a corresponding eigenvector v. We
then update g as g ← g − (v > g/||v||2 )v. To enforce the hard case we set ∆ to belarge: ∆ = 103 . Indeed we verified that λ∗ = µn in all the examples thus generated
indicating the hard case.
We set the convergence criteria of the CG method as follows: the relative residualis less than 10−12 , or the maximum allowed number of iterations 3000 is reached.
Since galahad gltr appeared to return an incorrect output in the hard case, here we
do not present its accuracy. The results are shown in Figure 5.
Observe that, especially in the hard case, our algorithm is fast and reliable, giving

288

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA

10

10

2

10

0

FRW
RSS
GEP
GLTR

FRW
RSS
GEP

10-5
accuracytime(s)

100
4

10-10

10-15
10

-2

103

104

105

106

103

matrix size

104

105

106

matrix size

Fig. 5. Hard case. Runtime (left) and accuracy (right), B = I and A sparse.
Table 1
Error in computed objective value.

n = 100
n = 1000
n = 10000

FRW
1.50e-12
1.49e-12
4.39e-14

RSS
9.37e-14
8.90e-14
8.32e-14

GEP
1.44e-15−6.22e-15
3.87e-14

the most accurate solutions. We note that while GALAHAD’s GLTR underperformedin this example, the latest prereleased version 3.0 has addressed the inefficiency andaccuracy issue for the hard case.
Problems whose exact solution is known. We can generate a TRS problem in thehard case with known exact solutions as follows: first, set
A = diag(−1, 2, 3, . . . , n),g = (0, −3α∆, 0, . . . , 0)> ,B = I, ∆ = 1, and α = 10−2 . In this case, the optimal value of TRS (1) is −(1 +3α2 )∆2 /2. We then generate an orthogonal matrix Q by the MATLAB commandqr(rand(n)) and update A, g as A ← QAQ> , g ← Qg, which does not change theoptimal objective value.
We measured the difference between the computed optimal value of each method
and the exact optimal value. Note that we sometimes obtained objective value slightly
smaller than the exact optimal value; this is caused by roundoff error (recall that we
impose that the obtained solution is feasible to O(u)). The result is shown in Table
1, which confirms that our algorithm indeed computes solutions with high accuracy,
and the previous accuracy plots give reliable estimates for the errors.
Approximate solutions. Recalling the discussion in section 6.3, here we consider
using our approach to efficiently obtain an approximate TRS solution by solving thegeneralized eigenvalue problem (28) approximately instead of to high accuracy O(u).
A natural strategy would be to terminate the iteration for computing the eigenpair
before convergence to full precision is attained.
To illustrate this idea, we formed an n = 1000 example and took B = I and let
A be a random matrix defined by A=randn(n);A=A’+A, and we apply k steps of thegg > 
Arnoldi process [1, Chap. 7] to −M1−1 M0 = −A ∆2 and approximate the eigenvecB−A

tor in (28) by the Ritz vector (e.g., [1, Chap. 3]) corresponding to the rightmost Ritz
value, from which we obtain approximate solutions xk for each k.
Figure 6 shows the accuracy of the resulting approximate solutions pk as the289

TRS VIA GEP

100

10-5

(pk , p∗ )||M0 xk − λM1 xk ||f (pk ) − f (p∗ )
6

10

-10

10-15
20

40
60
Krylov dimension

80

100

Fig. 6. Accuracy in solution ∠(pk , p∗ ) and objective value f (pk ) − f (p∗ ) as Krylov subspacedimension k varies. Also shown is the residual kM0 x − λM1 xk2 of the approximate eigenpair.

dimension of the Arnoldi subspace k varies. Here we quantify the accuracy from three
aspects: denoting by p∗ the solution (obtained by Algorithm 2) the angular error|p> p |∠(pk , p∗ ) = arccos kpkkkkp∗∗ k , the distance from the optimal objective value f (pk ) −f (p∗ ), and the eigenpair residual kM0 x − λM1 xk2 , which measures the quality of theapproximate eigenpair.
Observe in Figure 6 that the computed TRS solution pk improves its accuracy
as the quality of the Ritz vector improves. Since Ritz vectors are known to typicallyconverge geometrically, which we can see in the figure, this suggests that our approachis also attractive for computing an approximate TRS solution.
Summary of experiments. The results of our experiments can be summarized
as follows.
• Algorithm 2, GEP, based on one generalized eigenproblem is consistently
reliable and its accuracy is often among the highest, including the hard case.
• Algorithm 2 has speed comparable to the state-of-the-art methods, largely
independent of the problem: dense/sparse and/or B 6= I.
While GLTR was often the most efficient in the nonhard case, our algorithm is thesimplest to implement (largely due to the sophisticated but black-box eigensolvers
available via eig and eigs), robust and reliable without being much slower, tends toscale well as n grows and outperforms existing MATLAB implementations.
8. Discussion. As we have seen, the real eigenvalues of M (λ) correspond to the
KKT points for the TRS, and the largest eigenvalue provides a solution that minimizes
the objective function g > p + 12 p> Ap. In fact, more can be said: as shown by Forsythe
and Golub [7], the objective function value is an increasing function of λ, so we canalso maximize the objective function by finding the smallest real eigenvalue of M (λ).
Further analysis of the KKT points is given in [21]. The fact that we can bothmaximize and minimize the objective function is perhaps unsurprising as we impose
no positive definiteness assumption on A, so the objective function is nonconvex andthere is no fundamental difference between minimizing and maximizing it.
Future directions include performance benchmarking on parallel systems in thecontext of solving the TRS approximately in the trust-region method for nonlinear
optimization problems, comparing in particular with [13, 36]. Also worth considering
is extending the eigenvalue-based approach to solve other trust-region type problems

290

S. ADACHI, S. IWATA, Y. NAKATSUKASA, AND A. TAKEDA[3, 29], and dealing with a general quadratically constrained quadratic programming(QCQP) with one constraint. We note that an eigenvalue-based algorithm for QCQP
with two constraints is developed in [34], which also mentions in its appendix analgorithm for one constraint. However, that algorithm involves the computation of all
the eigenvalues, and thus requires O(n3 ) flops in all cases. It remains open to develop
a more efficient algorithm that exploits structure such as sparsity.

Acknowledgments. We are grateful to Nick Gould for providing comments on
the manuscript and project, and helping us with the Galahad library. We thank thereferees for their useful comments and suggestions. We thank Bill Hager and Henry
Wolkowicz for discussions and testing our codes.
REFERENCES[1] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, Templates for the
Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia, 2000,
https://doi.org/10.1137/1.9780898719581.
[2] S. P. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge, 2004.
[3] S. Burer and K. M. Anstreicher, Second-order-cone constraints for extended trust-region
subproblems, SIAM J. Optim., 23 (2013), pp. 432–451, https://doi.org/10.1137/110826862.
[4] A. R. Conn, N. I. M. Gould, and P. L. Toint, Trust Region Methods, Volume 1. SIAM,
Philadelphia, 2000, https://doi.org/10.1137/1.9780898719857.
[5] J. B. Erway and P. E. Gill, A subspace minimization method for the trust-region step, SIAM
J. Optim., 20 (2009), pp. 1439–1461, https://doi.org/10.1137/08072440X.
[6] J. B. Erway, P. E. Gill, and J. D. Griffin, Iterative methods for finding a trust-region step,
SIAM J. Optim., 20 (2009), pp. 1110–1131, https://doi.org/10.1137/070708494.
[7] G. E. Forsythe and G. H. Golub, On the stationary values of a second-degree polynomial
on the unit sphere, J. Soc. Indust. Appl. Math., 13 (1965), pp. 1050–1068, https://doi.org/
10.1137/0113073.
[8] C. Fortin and H. Wolkowicz, The trust region subproblem and semidefinite programming, Optim. Methods Softw., 19 (2004), pp. 41–67, https://doi.org/10.1080/
10556780410001647186.
[9] C. Fortin and H. Wolkowicz, Trust region subroutine algorithm: Algorithm and Documentation, 2010, http://www.math.uwaterloo.ca/∼hwolkowi/henry/software/trustreg.d.
[10] W. Gander, G. H. Golub, and U. von Matt, A constrained eigenvalue problem, Linear
Algebra Appl., 114 (1989), pp. 815–839, https://doi.org/10.1016/0024-3795(89)90494-1.
[11] D. M. Gay, Computing optimal locally constrained steps, SIAM J. Sci. and Stat. Comput., 2(1981), pp. 186–197, https://doi.org/10.1137/0902016.
[12] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., The Johns Hopkins University Press, Baltimore, MD, 2013.
[13] N. I. M. Gould, S. Lucidi, M. Roma, and P. L. Toint, Solving the trust-region subproblem
using the Lanczos method, SIAM J. Optim., 9 (1999), pp. 504–525, https://doi.org/10.
1137/S1052623497322735.
[14] N. I. M. Gould, D. Orban, and P. L. Toint, GALAHAD, a library of thread-safe fortran 90
packages for large-scale nonlinear optimization, ACM Trans. Math. Software, 29 (2003),
pp. 353–372, https://doi.org/10.1145/962437.962438.
[15] N. I. M. Gould, D. P. Robinson, and H. S. Thorne, On solving trust-region and otherregularised subproblems in optimization, Math. Program. Comput., 2 (2010), pp. 21–57.
[16] W. W. Hager, Minimizing a quadratic over a sphere, SIAM J. Optim., 12 (2001), pp. 188–208,
https://doi.org/10.1137/S1052623499356071.
[17] N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, 2008,https://doi.org/10.1137/1.9780898717778.
[18] G. Huber, Gamma function derivation of n-sphere volumes, Amer. Math. Monthly, 89 (1982),
pp. 301–302.
[19] S. Iwata, Y. Nakatsukasa, and A. Takeda, Global optimization methods for extended Fisher
discriminant analysis, in Proceedings of the Seventh AISTATS, JMLR W&CP 33, 2014,
pp. 411–419.

TRS VIA GEP

291[20] R. Lehoucq, D. C. Sorensen, and C. Yang, Arpack User’s Guide: Solution of Large-Scale
Eigenvalue Problems With Implicitly Restorted Arnoldi Methods, Vol. 6, SIAM, Philadelphia, 1998, https://doi.org/10.1137/1.9780898719628.
[21] S. Lucidi, L. Palagi, and M. Roma, On some properties of quadratic programs with a convex
quadratic constraint, SIAM J. Optim., 8 (1998), pp. 105–122, https://doi.org/10.1137/
S1052623494278049.
[22] J. M. Martı́nez, Local minimizers of quadratic functions on Euclidean balls and spheres, SIAM
J. Optim., 4 (1994), pp. 159–176, https://doi.org/10.1137/0804009.
[23] J. J. Moré, Recent developments in algorithms and software for trust region methods, Mathematical Programming: the State of the Art, Springer, Berlin, 1983, pp. 258–287.
[24] J. J. Moré, Generalizations of the trust region problem, Optimization Methods and Software,
2 (1993), pp. 189–209, https://doi.org/10.1080/10556789308805542.
[25] J. J. Moré and D. C. Sorensen, Computing a trust region step, SIAM J. Sci. and Stat.
Comput., 4 (1983), pp. 553–572, https://doi.org/10.1137/0904038.
[26] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer-Verlag, New York,1999.
[27] C. C. Paige, B. N. Parlett, and H. A. Van der Vorst, Approximate solutions and eigenvalue
bounds from Krylov subspaces, Numer. Linear Algebra Appl., 2 (1995), pp. 115–133, https://doi.org/10.1002/nla.1680020205.
[28] C. C. Paige and M. A. Saunders, Solution of sparse indefinite systems of linear equations,
SIAM J. Numer. Anal., 12 (1975), pp. 617–629, https://doi.org/10.1137/0712047.
[29] T. K. Pong and H. Wolkowicz, The generalized trust region subproblem, Comput. Optim.
Appl., 58 (2014), pp. 273–322, https://doi.org/10.1007/s10589-013-9635-7.
[30] F. Rendl and H. Wolkowicz, A semidefinite framework for trust region subproblems withapplications to large scale minimization, Math. Programming, 77 (1997), pp. 273–299,
https://doi.org/10.1007/BF02614438.
[31] E. Rimon and S. P. Boyd, Obstacle collision detection using best ellipsoid fit, J. Intell. Robot.
Syst., 18 (1997), pp. 105–126, https://doi.org/10.1023/A:1007960531949.
[32] M. Rojas, S. A. Santos, and D. C. Sorensen, A new matrix-free algorithm for the largescale trust-region subproblem, SIAM J. Optim., 11 (2001), pp. 611–646, https://doi.org/
10.1137/S105262349928887X.
[33] M. Rojas, S. A. Santos, and D. C. Sorensen, Algorithm 873: LSTRS: MATLAB software
for large-scale trust-region subproblems and regularization, ACM Trans. Math. Softw., 34(2008), 11, https://doi.org/10.1145/1326548.1326553.
[34] S. Sakaue, Y. Nakatsukasa, A. Takeda, and S. Iwata, Solving generalized CDT problems
via two-parameter eigenvalues, SIAM J. Optim., 26 (2016), pp. 1669–1694, https://doi.
org/10.1137/15100624X.
[35] D. C. Sorensen, Minimization of a large-scale quadratic function subject to a spherical constraint, SIAM J. Optim., 7 (1997), pp. 141–161, https://doi.org/10.1137/
S1052623494274374.
[36] T. Steihaug, The conjugate gradient method and trust regions in large scale optimization,
SIAM J. Numer. Anal., 20 (1983), pp. 626–637, https://doi.org/10.1137/0720042.
[37] G. W. Stewart and J.-G. Sun, Matrix Perturbation Theory, Computer Science and Scientific
Computing, Academic Press, Boston, 1990.
[38] P. D. Tao and L. T. H. An, A DC optimization algorithm for solving the trust-region subproblem, SIAM J. Optim., 8 (1998), pp. 476–505, https://doi.org/10.1137/S1052623494274313.
[39] P. L. Toint, Towards an efficient sparsity exploiting Newton method for minimization, Sparse
Matrices and Their Uses, I. S. Duff, ed., Academic Press, London, 1981, pp. 57–88.