Minimizing a Quadratic Over a Sphere.
A new method, the sequential subspace method (SSM), is developed for the problem of minimizing a quadratic over a sphere. In our scheme, the quadratic is minimized over a subspace which is adjusted in successive iterations to ensure convergence to an optimum. When a sequential quadratic programming iterate is included in the subspace, convergence is locally quadratic. Numerical comparisons with other recent methods are given.
Introduction
In this paper we consider the problem of minimizing a quadratic over a sphere:
subject to kxk
where A is a symmetric n  n matrix, b 2 R n , T denotes transpose, and k  k
is the Euclidean norm. This minimization problem is often called the trust region
subproblem since it must be solved in each step of a trust region algorithm [1, 2, 3, 15,
19]. Problems of this form arise in many other applications including regularization
methods for ill-posed problems [14, 26], and graph partitioning problems [10].
Although the solution to (1) can be expressed in terms of a diagonalization of A,
this representation is practical only when n is small. In this paper, we focus on the
large-scale case. One approach to the large-scale case, developed by Golub and von
Matt in [5] (also see [4]), is to (partially) tridiagonalize A using the Lanczos process
and then solve tridiagonal problems to obtain an approximate solution to (1). For
further developments of this approach, including preconditioning and a Fortran 90
implementation HSL VF05 in the Harwell Subroutine Library, see Gould et al. [7].
For the method developed in this paper, we use an approach in the spirit of the
Golub/von Matt/Gould et al. scheme to obtain a starting guess.
Parametric eigenvalue approaches to the sphere constrained problem (1) are developed
by Sorensen [24] and by Rendl and Wolkowicz [20]. The relationship between
these two approaches is discussed in detail in [20]. Roughly, Sorensen's approach
involves constructing an approximation to the solution of (1) from the solution to
a related eigenvalue problem. Since this approximation may not satisfy the bound
on the norm of the solution, a series of eigenvalue problems are solved, and in the
limit, the bound on the norm of the solution is fullled. In the approach of Rendl
and Wolkowicz, the same eigenvalue problem is solved in each iteration, however, the
bound on the norm of the solution is satised by maximizing a related dual func-
tion. The eigenvalue problems arising in either approach can be solved using Arnoldi
techniques such as those developed in [13]. In the \hard case" (see [16]) where b is orthogonal
to the eigenvectors associated with the smallest eigenvalue of A, Sorensen's
approach needs to be modied. An e-cient algorithm for the hard case is developed
by Rojas in her thesis [21]. She also uses this algorithm to solve some di-cult ill-posed
problems of Hansen [11, 12]. The approach of Rendl and Wolkowicz does not
need modication in the hard case, however, the convergence of algorithms for the
eigenvalue problem may be slower when the computed eigevalue is not simple.
The approach in this paper, which we call the sequential subspace method (SSM),
involves solving (1) with the additional constraint that x is contained in a subspace.
We show that convergence is locally quadratic (locally cubic when if the
subspace contains the iterate generated by one step of the sequential quadratic programming
(SQP) algorithm applied to (1). The convergence is quadratic even when
the original problem is degenerate with multiple solutions, and with a singular Jacobian
for the rst-order optimality system. Descent of the cost at a nonoptimal point
can be ensured by including in the subspace either the cost gradient or an eigenvector
associated with the smallest eigenvalue of A. We observe in numerical experiments
that appropriate small dimensional subspaces are generated by preconditioned
Krylov space and minimum residual techniques. Comparisons with the algorithms of
Sorensen [24], of Rendl and Wolkowicz [20], and of Gould, Lucidi, Roma, and Toint
are given in Section 5.
A solution of the problem
subject to
is any eigenvector associated with the smallest eigenvalue of A. In comparing the SSM
approach to algorithms for solving the eigenproblem, it follows from the discussion
of Sleijpen and Van der Vorst in [22] that an SQP iterate for (2) is closely connected
with the Rayleigh quotient iteration [18, p. 70], which is cubically convergent [18, p.
73]. In [22] approximate solutions to the SQP system are used to build up subspaces
containing the approximation to the eigenvector. In this paper, we solve the SQP
system relatively precisely, and we form a small dimensional subspace containing the
SQP iterate. After computing the new approximation in the subspace, the previous
information is discarded; hence, the computer memory requirements are relatively
small.
Complete diagonalization
If there exists a solution y of (1) with kyk < r, then A is positive semidenite and
y is the global minimizer of the quadratic x T Ax 2b T x. Thus, when a minimizer
of (1) lies in the interior of the constraining sphere, the constraint can be ignored
and the optimization problem can be approached using techniques for unconstrained
optimization. Consequently, we restrict our attention to the following equality constrained
problem:
subject to
The solutions to (3) are characterized by the following result (see [23, Lem. 2.4, Lem.
Lemma 1. The vector x is a solution of (3) if and only if r and there
exists  such that A + I is positive semidenite, and
The solution to (3) can be expressed in terms of the eigenpairs of A. Let
T be a diagonalization of A where  is a diagonal matrix with diagonal elements
and  is the matrix whose columns  1 ,
are
orthonormal eigenvectors of A. Dening
Lemma 2. The vector
is a solution of (3) if and only if c is chosen
in the following way:
(a) Degenerate case: If
then
are arbitrary scalars satisfying the condition
(b) Nondegenerate case: If (a) does not hold, then c
where  >  1 is chosen so that
Proof. Simply check that the su-cient optimality conditions of Lemma 1 are
satised. The degenerate case, where the Jacobian of the rst-order optimality system
may be singular, coincides with the \hard case" of More and Sorensen [16] where b
is orthogonal to the eigenspace associated with the smallest eigenvalue of A and the
multiplier  is equal to  1 . In the nondegenerate case, the multiplier  is chosen so
that A+ I is positive denite and the solution
the constraint x T
In the nondegenerate case, equation (5) leads to upper and lower bounds for the
multiplier . Since  i
kbk
r
To obtain a lower bound, observe that
which yields the relation
Utilizing the upper and lower bound  u and  l and the strict convexity of the left
side of (5) on the interval ( l ;  u
it is easy to devise e-cient algorithms to compute
a solution  of (5).
3 Incomplete diagonalization, local convergence
At iteration k in the sequential subspace method (SSM) for (3), we impose the additional
constraint that x lies in a subspace S k of R n . Hence, the new iterate x k+1 is a
solution of the problem
subject to
We show that the convergence is locally quadratic, even when the original problem
(3) is degenerate, if we include an SQP iterate associated with x k
in S k
If V is an n  l matrix with orthonormal columns that span S k , then (8) is
equivalent to the problem
subject to
After substituting for x, (9) reduces to the following problem in R l :
subject to
l is small, then (10) can be solved by complete
diagonalization as in Section 2. Or if B has a sparse factorization, then (10) can be
solved quickly using the Newton approach developed in [16].
In theory, a tridiagonal B is generated using the Lanczos process [6]. In particular,
if v 1 is a unit vector and v i
is the i-th column of V, then the Lanczos process can be
expressed as follows:
Algorithm 1 (Lanczos)
ksk
Here d is the diagonal and u is the superdiagonal of the tridiagonal matrix B. If
then the Lanczos process is terminated and the column space of
V and AV coincide.
It is well-known that the columns of V generated by this process may deviate
signicantly from orthogonality due to the propagation of rounding errors. And when
this happens, (9) is no longer equivalent to (10). Nonetheless, Gould et al. observe
in [7] that the solution to (10) often provides a good approximation to the solution of
despite the loss of orthogonality. The Lanczos process can be repaired, in order to
restore orthogonality, by using a Householder process to generate the columns of V.
This process, however, requires products between a vector and each of the previously
computed columns of V. Thus the overhead needed to maintain orthogonality grows
like nl 2 in the number of
ops and like nl in storage. This overhead can be signicant
when n or l is large. On the other hand, to compute a high accuracy solution, we
need to maintain orthogonality in order to obtain an equivalent problem (10). This
leads us to focus on approaches that involve subspaces where l is much smaller then
n. In particular, for an implementation (Algorithm 4) of the SSM proposed later, l
is either 4 or 5.
Since sequential quadratic programming (SQP) techniques often converge rapidly,
with a good starting guess, we always include the SQP approximation in the subspace
. If x k
is the current iterate, which we assume satises the constraint
and  k is the current approximation to the multiplier associated with the constraint,
then the SQP iterate can be expressed in the following way:
and
z and  are solutions of the following linear system:
When the coe-cient matrix in (11){(12) is singular, we let (z; ) be the minimum
residual/minimum norm solution; that is, (z; ) is gotten (in theory) by multiplying
the right side by the pseudoinverse of the coe-cient matrix (see [8]). The SQP method
is equivalent to Newton's method applied to the nonlinear system
A solution x k+1 to the subspace problem (8) is an approximation to the solution of
(3). To obtain an estimate for the multiplier of Lemma 1, we minimize the Euclidean
norm of the residual b Ax k+1 x k+1 over the scalar . This works out to give
(b Ax) T x
We now examine the local convergence of a solution x
of (8) and the multiplier
estimate (14) under the assumption that S k contains
is a solution to (11). Let S  denote the set of minimizers of (3), and let   be the
multiplier given by Lemma 1. In the nondegenerate setting where A+  I is positive
denite, we show that the iteration is locally, quadratically convergent to the unique
solution of (3). In the degenerate case   =  1 where S  has more than one element,
we obtain local quadratic convergence to S  , where distance is measured in the usual
way:
In the nondegenerate degenerate-case where  contains a single
element, we obtain local quadratic convergence for a \safe-guarded" choice of  k .
Our convergence result in the special nondegenerate degenerate-case is given later in
Lemma 5, while our local convergence result in either the nondegenerate case or the
degenerate case with multiple solutions is the following:
Theorem 1. Let   be the multiplier of Lemma 1 associated with the set of
solutions S  of (3), and suppose that either A+   I is positive denite, or
with (4) a strict inequality. Then there exist positive constants  and C with the
property that for any (x k ;  k ) such that
and for any subspace S k that contains the SQP iterate xSQP associated with (11){
(12), any solution x k+1 of (8) and associated multiplier  k+1 given by (14) satisfy the
following estimate:
The eigenvalue problem (2), corresponding to b = 0, is always degenerate (with
multiple solution) and the error has the following special form:
When the multiplier is estimated using (14), it can be shown, when that the
error in the multiplier is bounded by a constant times the error in the solution vector
squared (see the remark at the end of Section 3.1). It follows that for some constant
the same as the convergence result for the Rayleigh quotient iteration.
3.1 Nondegenerate problems
We begin the derivation of Theorem 1 with the nondegenerate case:
Lemma 3. If (3) has a solution x  and an associated multiplier   with   >  1 ,
then there exist a neighborhood N of (x  ;   ) and a constant C with the property that
for any and for any subspace S k that contains the
SQP iterate xSQP associated with (11){(12), any solution x k+1 of (8) and associated
multiplier  k+1 given by (14) satisfy the following estimate:
Proof. Since   >  1 , the matrix A+  I is positive denite, and the Jacobian
of the nonlinear system (13) is nonsingular at By the standard convergence
theorem for Newton's method applied to a smooth system of equations, there exist a
neighborhood N of
) and a constant c such that
whenever
Let  and  be positive scalars chosen so that
for all x 2 R n , let f be the cost function in (3), let L be
the Lagrangian dened by
A Taylor expansion around x  yields the following relation:
for any x 2 B r
rg. Combining this with (16) gives
for any x 2 B r
If p is the projection of xSQP onto B r , then
Hence, we have
xSQP for some
, it follows that p 2 S k and f(x k+1 )  f(p). Combining
this inequality with (17) and (18) gives
which implies that
Making this substitution gives
r
Combining (21) with (19), the proof is complete.
Remark. For the eigenvalue problem (2), we have x
In this case, (20) yields
and (21) becomes
3.2 Degenerate problems
Now consider local convergence in the degenerate case where  Referring to
Lemma 2, the degenerate case can only happen when  i
Any solution to (3) in the degenerate case can be expressed x
and  1 is any linear combination of the vectors  i
satisfying the relation
Initially, we suppose that k 1 0, in which case the projection of S  on the
eigenspace associated with contains a sphere of radius -. Our convergence result
is the following:
Lemma 4. Suppose that the multiplier   of Lemma 1 associated with the set of
solutions S  of (3) is given by  is the
component of an element of S  in the eigenspace associated with E 1 . Then there exist
positive constants  and C with the property that for any (x k ;  k ) such that
and for any subspace S k that contains the SQP iterate xSQP associated with (11){
(12), any solution x k+1 of (8) and associated multiplier  k+1 given by (14) satisfy the
following estimate:
Proof. Initially, let us assume that  k is near  1 , but  k 6=  1 . In this case,
the linear system (11){(12) is nonsingular, and there exists a unique solution (z; ).
We expand z and x k in terms of the eigenvectors of A writing
and x
. Utilizing (11), we obtain
Substituting this in (12) gives
Let us dene
I)x k and
. For
If x  2 S  , then since
we have
Let  k be the error at step k dened by
By (26), we have
since
is near - 2 > 0 when x k is near S  . From (23), we have
Let x  be the closest element of S  to x k and dene
. Then we have
By (30) the  i
component of in error by O( k ) since  i ,
the  i
component of x k
, is in error by O( k
by (31).
implies that
. Combining this with (28) and
Hence, for i 2 E+ the  i
component of xSQP is in error by O( 2
be the seminorm associated with projection into the eigenspace associated
with
Then we have
for all x ng. Proceeding as we did
earlier, but replacing norms by seminorms,
where p is the projection of xSQP onto the ball B r , and
xSQP for some
by (30) and (32), and z is perpendicular to x k
by (12), we have
This implies that
Consequently, kp x
which combines with (34) to give
By the triangle inequality,
Let k  k 1 be the seminorm dened by
and recall that kx  k . By the Pythagorean theorem and the fact
that x k+1
has length r, we have
which implies that
The distance from x k+1 to S  is given by
where x  is any element of S  . Relations (35){(38) yield dist(x
Combining these estimates, we have  k+1
This analysis was given under the assumption that  k 6=  1 . In the special
case  k
show how the analysis should be modied. With the change
of variables z =
and the substitution x
, the SQP system
(11){(12) is equivalent, by orthogonal transformation, to
D#
where D is a diagonal matrix with diagonal elements
then the rst s diagonal elements of D and the rst s components of  D vanish.
Hence, the rst s equations in (39) imply that The next n s equations give
while the last equation in (39) gives
The minimum norm solution to this last equation is
By (40),
Combining these bounds, we have these relations, all the analysis
from (33) onward can be applied, leading us to the estimate  k+1
Lemmas 3 and 4 yield Theorem 1.
3.3 Nondegenerate degenerate-problems
Finally, let us consider the nondegenerate degenerate-case where
and the  1 component of x  in the eigenspace associated with the smallest
eigenvalue of A vanishes. Our convergence result is the following:
Lemma 5. If (3) has a solution x is given by (22), then there
exist a neighborhood N of (x  ;  1 ) and a constant C with the property that for any
and for any subspace S k that contains the SQP iterate xSQP associated with (11){
(12), the solution x k+1 of (8) and associated multiplier  k+1 given by (14) satisfy the
following estimate:
In the case that  k
k, C can be chosen so that
Proof. Focusing on the numerator in (24), and substituting
, we
With this substitution for the numerator of  in (24), we obtain
the denominator terms in (43) have the following lower bound:
Another lower bound is gotten by neglecting terms corresponding to indices i
where the seminorm k  k 1 is dened in (36). Combining (43){(45) yields:
Returning to our previous analysis of the degenerate case, it follows from (29) and
(46) that for

Here we exploit the fact that for . In order to analyze (47), we
consider two separate cases: (i) kx k
x  k+ and (ii) kx k
x  k+ ,
where  is any xed constant satisfying
In case (i),
x  k+
x  k+
We now derive a similar bound for the left side of (49) in case (ii). In this case, it
follows from (42) that
x  k+
for any x 2 R n , where a + subscript on a vector is used to denote its projection on
the eigenspace associated with E+ . After substituting for  using (20), we obtain
for any x 2 B r
. Assuming x k
is a unit vector (note that when kx k
only if x
We will establish a uniform bound for the expression (51) when x k
is near x  ,
To facilitate this analysis, we rst consider
whether the equation
has a solution of the form
y of this form, the Schwarz inequality gives
jy T
Since the unit vector y is orthogonal to the eigenspace associated with  1 ,
Multiplying (52) by y T and using both (53) and (54) gives
For any x 2 B r , we have
which implies that
since yields the relation:
Referring to (55), we have a contradiction when kx x  k+  .
In summary, the equation (52) has no solution over the set Y consisting of those
y that satisfy the following conditions:
lies in the closure of Y, then by (57),
since any solution of (52) satises (55), y cannot be a solution of (52).
Since (52) has no solution over the closure of Y, the following constant - is strictly
positive:
Since
lim
min
(51) is bounded uniformly over all x k near x  with
x  k+ .
Thus in either case (i) or (ii), the left side of (49) is bounded, and by (47), we have
the same as relation (30) in the degenerate case.
To establish the analogue of (32) for indices i 2 E+ , we need a dierent bound for
the next to last term in (43). From the identity
Hence, we have
Since
implies that
It follows that
This estimate along with the lower bound (44) for the denominator in (43) yields the
relation
The reminder of the analysis is identical to that given for the degenerate case (Lemma
4), starting with (32). Since S it follows from the analysis of Lemma 4 that
In the special case  k
Hence, the in (60) can be absorbed in the kx k
completes the proof.
Implementation
In our experimentation with the SSM, we put the following four vectors in S k in each
iteration: xSQP , x k
, and an estimate for an eigenvector of A associated with
the smallest eigenvalue. By including x k in S k , the value of the cost function can only
decrease in consecutive iterations. The multiple b Ax k
of the cost function gradient
ensures descent if the current iterate does not satisfy the rst-order optimality
conditions. The eigenvector associated with the smallest eigenvalue will dislodge the
iterates from a nonoptimal stationary point. We also use this vector in a \safe-guard"
strategy designed to keep A positive denite.
4.1 The SQP system
Now consider the SQP system (11){(12). According to (12), z is orthogonal to the
prior iterate x k . Let P be the matrix that projects a vector into the space perpendicular
to
Multiplying (11) by P yields
according to (12), we have
We haved found that preconditioned Krylov space methods, such as the Gauss-Seidel
scheme in [9], converge very quickly when applied to (61). As a small illustra-
tion, let us consider the second test problem from [24] with
where  is a 1000  1000 diagonal matrix with diagonal elements selected randomly
from a uniform distribution on ( :5; :5) and I 2qq T where q is gotten by
rst generating random numbers on ( :5; :5) and then scaling the resulting vector
to have unit length. The vector b is generated in the same way as q. The solid
curve in Figure 1 gives the convergence when a Lanczos type process (Algorithm 1,
with starting vector v is used to generate the matrix V used in (9). The
Lanczos process was modied to ensure orthogonality of the columns of V. For each
value of l in Algorithm 1, we solve the l  l tridiagonal problem (10) to obtain an
approximate solution x and associated multiplier for the original problem
(3). In the solid curve of Figure 1, we plot the base 10 logarithm of the norm of the
residual kb (A+I)xk. According to Lemma 1, the residual vanishes at an optimal
solution.
The dashed curve of Figure 1, based on the SSM approach, is gotten in the following
way: Taking Algorithm 1, we generate a V with 40 orthonormal
columns. Solving (10), we obtain starting guess x 0 . In iteration k of the SSM phase,
we start with the vector v we use the Gauss-Seidel/Krylov space
approach of [9] to generate a matrix V, with orthonormal columns, that approximately
contains a solution of (61) in its range. Using the V generated in this way, we
solve (9) to obtain the next iterate x k+1
. The associated multiplier is estimated using
(14). Each kink in the dashed curve of Figure 1 corresponds to the number of iterations
needed to obtain an approximate solution of (61). In this example, roughly 15
multiplications by the elements of the matrix A are used to solve (61). The quadratic
convergence of SSM is re
ected in the rapid decay of the residual norm.
This approach for generating V, using a nonsymmetric Gauss-Seidel matrix,
Krylov spaces, and orthogonalization, can become expensive when n is really large
Matrix-vector products

Figure

1: Convergence of the tridiagonalization approach (solid) and SSM (dashed)
for the second test problem from [24].
since each of the columns of V should be stored in memory. Hence, in the remainder
of this paper, we focus on low-storage symmetric techniques for solving (61), which
we compare to other approaches.
We solve (61) using a preconditioned version of Paige and Saunders' MINRES
algorithm [17]. More precisely, we use Algorithms 3 and 3a in [9] and three dierent
choices for the symmetrizing preconditioner W in that paper: (i) corresponding
to unconditioned iterations, (ii) is the diagonal matrix
whose diagonal matches that of
L is the strictly lower triangular matrix whose lower triangle matches that of C. The
implementations of SSM associated with the latter two preconditioners are denoted
and SSM l
respectively.
Typically, the L matrix associated with
I)P is dense, even when A
is sparse, since P is often dense. Nonetheless, linear systems of the form
can be solved in time proportional to the number of nonzero elements in the lower
triangle of A due to the special structure of C. In terms of the vectors w, q and p
dened by
I)w and
the diagonal d of C can be expressed
while the o-diagonal elements of C are
Exploiting this structure, it can be shown that the solution to (L can
be computed in the following way:
Algorithm
The statement y a i+1:n;i of Algorithm 2 only requires the nonzero
elements in column i of A beneath the diagonal. Hence, the number of
oating point
operations for Algorithm 2 is O(n) plus the number of nonzero elements in the lower
triangle of A.
The analogous procedure for the transposed system is the following:
Algorithm 3
4.2 Positive deniteness
In theory, the MINRES algorithm we use to solve (61) can be applied to any symmetric
matrix. In practice convergence can be extremely slow when C is indenite. For this
reason, we try to choose  k so that A+  k I is positive denite. If e is an eigenvector
of the matrix B in (10) associated with the smallest eigenvalue , then the pair (v; ),
approximates an eigenpair of A corresponding to the smallest
eigenvalue. The error in  can be estimated in the following way: If  is closer to  1
than the other eigenvalues of A, then after substituting
in the residual r = Av v, we have
since
1. Thus j  1 j  krk, which implies that
With this insight, we replace the least squares estimate (14) by the following safeguarded
estimate:
When the approximate eigenpair (v; ) is not very accurate, then the safe-guarded
step (62) is a safe, but poor approximation to   . Hence, whenever
we apply one iteration of SSM to the quadratic eigenvalue problem (2) in order to
compute a more accurate eigenpair. Due to the third and sixth order estimates in
(15), simply one iteration of SSM for the eigenproblem often yields a highly accurate
eigenpair.
4.3 The algorithm
We now collect our observations and present the algorithm that was used to generate
the numerical results of the next section. To simplify the presentation, we introduce
the following subroutines:
This routine applies Algorithm 1 to the matrix A,
starting from the vector v 1 , to generate a matrix V with columns
l
This routine solves the problem (8) generating a
solution denoted x, and an associated multiplier is a matrix
whose columns are an orthonormal basis for S k , then an estimate (v; ) for the
smallest eigenvalue of A and an associated eigenvector is gotten by computing
the smallest eigenvalue  and an associated eigenvector e for
setting
This routine computes a (minimum residual, minimum
solution (z; ) of the following linear system:
Our implementation of the sequential subspace method combines these three routines
and the safe-guarded step (62):
Algorithm 4 (Safe-guarded SSM with Lanczos startup)
while (  ==  &
it
while ( kb
Algorithm 4
For the computational results reported in the next section, we took
:01ng. The \rand" function appearing at the start of Algorithm 4 generates
a vector with components uniformly distributed on [0; 1].
5 Computational results
In this section we compare the performance of SSM to the performance of the algorithms
in [7], [20], and [24], denoted GLRT, RW, and S respectively, using the three
test problems presented in [24]. The results that we report for S were extracted from
[24], while the results reported for GLRT and RW were obtained using codes provided
by the authors. We thank the authors for providing access to their codes. Each of
these codes used dierent stopping criteria. GLRT stopped when kb (A+I)xk=kbk
was bounded by a given tolerance, while RW stopped when the gap between the value
of the primal and dual problem, and hence the error in the primal cost function, was
smaller than a given tolerance. In order to ensure that each code computed a solution
with the same accuracy, we adjusted the error tolerance parameter of each code until
the value of kb for the computed solution was smaller than a given
tolerance (specied below).
In the rst test problem of [24], is the standard 2-D
discrete Laplacian on the unit square based on a 5-point stencil with equally-spaced
mesh points. Taking a series of 20 problems were
generated where b was a vector with elements uniformly distributed on [0; 1]. Each
of these problems was solved using three dierent tolerances,
In

Table

1 we give the average number of matrix-vector products involving A for each
algorithm. Each iteration of the preconditioned MINRES algorithm with lower trian-
Tolerance S RW GLRT SSM SSM d
SSM l

Table

1: Problem 1, average number of matrix-vector products versus tolerance.
gular preconditioner involves roughly twice as many
ops as an iteration of either the
identity or the diagonal preconditioned schemes. Hence, in doing the bookkeeping,
we charged for two matrix-vector products in each iteration of the triangular preconditioned
scheme. As seen in Table 1, SSM l converges more than twice as fast as the
identity and diagonal preconditioned schemes, and overall, SSM l
uses the smallest
number of matrix-vector products for this test problem. Since the parametric eigenvalue
algorithms S and RW compute an extreme eigenvalue for a series of matrices,
we also list in parentheses in Table 1 the number of these eigenproblems that are
solved. Hence, RW is very economical in terms of the number of these eigenproblems
that are solved.
The second suite of test problems in [24] utilizes the matrix described earlier in
Section 3. In these problems, the radius of the sphere is varied and the number of
matrix-vector products is tabulated. For radii of one or smaller, solutions can be
computed extremely quickly, so we focused on and an error
tolerance of 10 7 . In Table 2 we see that for
had the fewest matrix-
Radius S RW GLRT SSM SSM d SSM l

Table

2: Problem 2, average number of matrix-vector products versus radius.
vector products, while GLRT had the fewest for
The nal problem of [24] again employed the discrete Laplacian matrix, but with
100. The vector b was designed to make the problem degenerate;
rst a random b was generated, then its  1 component was removed. Table 3 gives
the results for the various algorithms.
SSM l

Table

3: Problem 3, average number of matrix-vector products.
We placed an asterisk by the result in Table 3 for GLRT since this routine reduced
the error to 10 4 , not the 10 7 tolerance used by the other routines. Among
the routines that achieved the error tolerance, SSM l
performed the best relative to the
number of matrix-vector products. Note that the number of matrix-vector products
given in Table 3 for S was taken from [24] while Rojas, in her recent thesis [21], developed
a more e-cient implementation of Sorensen's approach for degenerate problems.
In summary, a Lanczos type process seems to be very eective when the problem
is very nondegenerate (  >>  1 ). As the problem becomes more degenerate,
preconditioned schemes such as SSM d or SSM l appear more eective. The number of
times that RW computes an extreme eigenpair is often around 5. For the numerical
experiments reported in this paper, Matlab's eig routine was used to compute this
extreme eigenpair. If this routine for computing an extreme eigenpair could be sped
possibly using the Jacobi type methods of Sleijpen and Van der Vorst [22] or
the truncated RQ iteration of Sorensen and Yang [25], the number of matrix-vector
operations used in the parametric eigenvalue approach would be reduced.



--R

A trust region algorithm for nonlinearly constrained optimization
A trust region strategy for non-linear equality constrained optimization

Least squares with a quadratic constraint
Quadratically constrained least squares and quadratic problems
Matrix Computations
Solving the trust-region subproblem using the Lanczos Method
Applied Numerical Linear Algebra
Iterative methods for nearly singular linear systems
Graph partitioning and continuous quadratic programming
a MATLAB package for analysis and solution of discrete ill-posed problem


Geophysical Data Analysis: Discrete Inverse Theory


Solution of sparse inde
The Symmetric Eigenvalue Problem
A trust region algorithm for equality constrained optimization
A semide
A Large-scale Trust-region Approach to the Regularization of Discrete Ill-posed Problems
A Jacobi-Davidson iteration method for linear eigenvalue problems
Newton's method with a model trust region modi
Minimization of a large-scale quadratic function subject to a spherical constraint
A truncated RQ iteration for large scale eigenvalue calculations
Inverse Problem Theory
--TR

--CTR
Peter A. Graf , Wesley B. Jones, A projection based multiscale optimization method for eigenvalue problems, Journal of Global Optimization, v.39 n.2, p.235-245, October   2007
Stanislav Busygin, A new trust region technique for the maximum weight clique problem, Discrete Applied Mathematics, v.154 n.15, p.2080-2096, 1 October 2006
symmetric successive overrelaxation;sparse optimization;arnoldi m orthogonalization;gauss-seidel;minimal residual;large-scale optimization;trust region subproblem;quadratic programming;krylov space;preconditioning;quadratic optimization