# A Prediction-Correction Dynamic Method for Large-Scale Generalized Eigenvalue Problems

Authors: Xin-long Luo; Jia-ru Lin; Wei-ling Wu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper gives a new prediction-correction method based on the dynamical system of differential-algebraic equations for the smallest generalized eigenvalue problem. First, the smallest generalized eigenvalue problem is converted into an equivalent-constrained optimization problem. Second, according to the Karush-Kuhn-Tucker conditions of this special equality-constrained problem, a special continuous dynamical system of differential-algebraic equations is obtained. Third, based on the implicit Euler method and an analogous trust-region technique, a prediction-correction method is constructed to follow this system of differential-algebraic equations to compute its steady-state solution. Consequently, the smallest generalized eigenvalue of the original problem is obtained. The local superlinear convergence property for this new algorithm is also established. Finally, in comparison with other methods, some promising numerical experiments are presented.


---

## Body

## 1. Introduction

In this paper, we consider the smallest generalized eigenvalue problem, which is often encountered in engineering applications such as automatic control, dynamical analysis of structure, electronic structure calculations, and quantum chemistry (see [1, 2] and references therein). For this old and active problem, recently, Gao et al. gave an interesting continuous projected method [3, 4]. This article follows this line and gives a new prediction-correction method based on the dynamical system of differential-algebraic equations (DAEs) for this problem.First, we convert the smallest eigenvalue problem into an equivalent equality-constrained optimization problem. Second, from the Karush-Kuhn-Tucker conditions of this optimization problem, we obtain a variant of the Rayleigh quotient gradient flow [5, 6], which is formulated by a system of DAEs. Third, applying the implicit Euler method [7], a projected technique [4, 8, 9], and a new time-stepping strategy to that dynamical system of DAEs, we construct a new prediction-correction method to compute a steady-state solution of that special dynamical system. Consequently, we obtain the smallest generalized eigenvalue of (A,B). We also establish the local superlinear convergence property for this new method. Last, in comparison with other methods, some promising numerical experiments are presented. Throughout this article, ∥·∥ denotes the Euclidean vector norm and the corresponding induced matrix norm.


## 2. Continuous Dynamical Model

The generalized eigenvalue problem(A,B) is to find a scalar λ and nonzero vector x∈ℜn to satisfy
(1)Ax=λBx,
where A and B are n×n real symmetric matrices and B is positive definite. Problem (1) can be converted to an equivalent equality-constrained optimization problem [4]:
(2)minx∈ℜnxTAx,s.t.xTBx=1.
If x* is an optimal solution of problem (2), according to the Karush-Kuhn-Tucker conditions [10, p. 328], then there exists a Lagrange multiplier λ* such that the following conditions are satisfied at (x*,λ*):
(3)Ax*-λ*Bx*=0,(x*)TBx*=1.
Therefore, the solution x* of problem (2) is a generalized eigenvector of (A,B).Property 1 (Gao et al., 2008 [4]).
Ifx* is an optimal solution of problem (2) and λ* is the corresponding Lagrange multiplier, then they satisfy the second-order optimality condition
(4)(A-λ*B)≽0.
Furthermore, a local minimizer of (2) is also its global minimizer. If x* is a global minimizer of (2), then it is an eigenvector associated with the smallest generalized eigenvalue of (1).From the first-order optimality conditions (3) of problem (2), we construct the following negative gradient flow:
(5)dx(t)dt=-(Ax-λBx),
on the generalized spherical surface,
(6)xTBx=1.
The system of DAEs (5)-(6) is a special dynamical system which is derived from the equality-constrained optimization problem (2). According to the definition of the differential index of DAEs [11], we can verify that the index of DAEs (5)-(6) is two. Therefore, the solution of DAEs (5)-(6) is not easy to be obtained by the general software package such as bdf15s [7, 11, 12]. Thus, we need to devise a special technique for its steady-state solution.Differentiating the algebraic constraintxTBx=1 along the trajectory (x(t),λ(t)) of (5) and (6), we obtain
(7)d(xTBx)dt=(∇x(xTBx))Tdxdt=2xTB(-Ax+λBx)=2(-xTBAx+λxTB2x)=0.
Consequently, we obtain
(8)λ=xTBAxxTB2x.
Substituting (8) into (5), we also obtain the following generalized gradient flow:
(9)dx(t)dt=-(A-BxTBAxxTB2x)x,x(0)TBx(0)=1.
Conversely, if x(t) is a solution of (9), then d(xTBx)/dt=0; that is, x(t) satisfies xTBx=1, and consequently x(t) also satisfies (5)-(6). Therefore, the continuous dynamical system of (5) and (6) is equivalent to the generalized gradient flow (9).Since the right-hand-side function of (9) is Lipschitz continuous, it has a unique solution x(t) for any initial point x(0)=x0. We can verify that the objective function f(x)=xTAx is monotonically decreasing along the trajectory x(t) of (5)-(6). Actually, from (5)–(8), using the Cauchy-Schwartz inequality |yTz|≤∥y∥∥z∥, we have
(10)df(x)dt=(∇xf(x))Tdxdt=2(-xTA2x+xTABxxTBAxxTB2x)≤2xTB2x(-∥Ax∥2∥Bx∥2+(Ax)T(Bx)(Bx)T(Ax))≤0.
Furthermore, the solution x(t) of (9) converges to an eigenvector x* associated with the smallest generalized eigenvalue λ* of (A,B), for almost all initial condition x(0)=x0, where x0 satisfies ∥x0∥=1, except for a set of measure zero, that is, the case of a multiple generalized eigenvalue λ* [5]. Therefore, by following the trajectory x(t) of (5)-(6), we can compute a global minimum point x* of problem (2), and consequently obtain the smallest generalized eigenvalue λ* of (A,B).


## 3. A Prediction-Correction Method

In this section, in order to obtain a steady-state solution of DAEs (5)-(6), we construct a prediction-correction method to follow its trajectory. According to the discussion of Section 2, we know that a steady-state solution of DAEs (5)-(6) is a generalized eigenvector associated with the smallest generalized eigenvalue of (A,B). Therefore, by following the discrete trajectory of DAEs (5)-(6), we can compute an approximate eigenvector associated with the smallest generalized eigenvalue and consequently obtain the smallest generalized eigenvalue.Note that we mainly consider the steady state of the system of DAEs (5)-(6) and do not care about an accurate solution at its transient-state phase. In order to avoid consuming unnecessary computing time, we adopt the first-order implicit Euler method to follow its trajectory. The time steps of the implicit Euler method are not restricted by the absolute stability property for the linear test equation dx/dt=-μx, where μ>0 [7], and consequently the implicit Euler method can take large steps at the steady-state phase so that the iteration sequence {xk} converges rapidly to a stationary point x* of the system of DAEs (5)-(6).By applying one implicit Euler iteration to the system of DAEs (5)-(6), we obtain the following iteration formula of the form:
(11)xk+1=xk-Δtk(A-λk+1B)xk+1,(12)xk+1TBxk+1=1,
where Δtk is the time step. Note that the system of (11)-(12) is nonlinear, and generally, its solution cannot explicitly be obtained. Therefore, we use the alternating projection idea [13, 14] to handle it.We replaceλk+1 with λk in (11) to obtain the predicted point xk+1P as follows:
(13)(I+Δtk(A-λkB))xk+1P=xk.
Since the predicted point xk+1P moves away from the elliptic spherical surface, we project this predicted point xk+1P onto the unit elliptic spherical surface, so that the iterated point xk+1 satisfies the algebraic constraint (6); namely, we find the shortest distance between xk+1P and the unit elliptic spherical surface. We achieve this aim by solving the following equality-constrained problem:
(14)minx∈𝔑n∥x-xk+1P∥,s.t.xTBx=1.Since the optimal solution of the previously problem is not apparently obtained, we turn to achieve its suboptimal solution. Since matrixB is positive, it can be decomposed as B=LTL. Let y=Lx. Then, problem (14) is equivalent to the following problem:
(15)miny∈𝔑n∥L-1(y-Lxk+1P)∥,s.t.∥y∥=1.
Since ∥L-1(y-Lxk+1P)∥≤∥L-1∥∥y-Lxk+1P∥, we obtain the suboptimal solution by solving the following problem:
(16)miny∈𝔑n∥(y-Lxk+1P)∥,s.t.∥y∥=1.
From the shortest distance between a point and a unit spherical surface being the intersection point of the straight line through the center of that circle and the point and that unit spherical surface, we obtain the previous minimum point as yk+1=Lxk+1T/∥Lxk+1P∥. Then, the suboptimal solution point of problem (14) is obtained as
(17)xk+1=L-1yk+1=xk+1P∥Lxk+1P∥.Now, we use the correction pointxk+1 to update λk+1, so that the steady-state equation Ax-λBx=0 is satisfied as possible. Namely, we find the minimum point of the following problem:
(18)minλ∈𝔑n∥Axk+1-λBxk+1∥.
Since matrix B is positive, it can be decomposed as B=LTL. Let zk+1=Lxk+1. Then, the previous minimum problem is equivalent to the following problem:
(19)minλ∈𝔑n∥LT(L-TAxk+1-λLxk+1)∥.
Using the property ∥LT(L-TAxk+1-λLxk+1)∥≤∥LT∥∥(L-TAxk+1-λLxk+1)∥, we obtain the suboptimal solution of problem (18) by solving the following problem:
(20)min∥L-TAxk+1-λLxk+1∥.
It is not difficult to obtain the minimum solution of the previous problem as
(21)λk+1=xk+1TAxk+1xk+1TBxk+1,
which is the suboptimal solution of problem (18).Another issue is how to choose the time stepΔtk at every iteration. We adopt an analogous trust-region technique to adaptively adjust the time step Δtk. Our intuition is that the time step Δtk can be enlarged when the predicted point xk+1P is also near the elliptic spherical surface xTBx=1; otherwise, the time step Δtk is reduced. In order to measure the distance between xk+1P and the unit elliptic spherical surface, we construct the following approximate model:
(22)ρk=|1-((xk+1P)TBxk+1P)1/2|.
Then,
(23)Δtk+1={γ1Δtk,ifρk<η1,Δtk,ifη1≤ρkη2,γ2Δtk,ifρk≥η2,
where constants η1, η2, γ1, and γ2 satisfy
(24)0<η1<η2<1,0<γ2<1<γ1.According to the previous discussion, we give the following discrete dynamical method for the smallest generalized eigenvalue of(A,B).Algorithm 1.
Prediction-correction method for the smallest generalized eigenvalue.
Step  1. Initialize the parameters. Specify a tolerated error TOL and an initial point x0 to satisfy
(25)x0TBx0=1.
Compute λ0=x0TAx0/x0TBx0 and the residual r0=∥Ax0-λ0Bx0∥.
Step  2. Repeat for k=0,1,2,… until rk=∥Axk-λkBxk∥≤TOL:  (a)
solve (13) to obtain the predicted point xk+1P;(b)
evaluate (17) to obtain the corrected point xk+1;(c)
compute (21) to update the Lagrange multiplier λk+1;(d)
compute the residual(26)rk+1=∥Axk+1-λk+1Bxk+1∥;(e)
update the time stepΔtk+1 according to the time-stepping scheme (22)-(23);(f)
accept the trial vector. Ifρk<η2, let xk=xk+1, λk=λk+1 and rk=rk+1;(g)
try again starting at (a) withk←k+1.


## 4. Local Convergence Analysis

The sequence(xk,λk) generated by Algorithm 1 follows the continuous trajectory (5)-(6). Adopting the similar estimation technique in [15, 16], we obtain that the local truncation error of xk is O(Δtk). According to the standard analysis of numerical ordinary differential equations [7, 11], we know that xk is close to x(tk) when Δtk is small enough in the finite time interval [0,T]. Therefore, we only need to estimate the local convergence rate of xk in the steady-state phase.For the convenience of analysis, we decompose matrixB as B=LTL and denote C=L-TAL-1. Then, the smallest generalized eigenvalue λ* of (A,B) and associated with the eigenvector x* is equivalent to the smallest eigenvalue of symmetric matrix C and associated with the eigenvector y*=Lx*, respectively. Therefore, we only need to consider the convergence property of (yk,λk), where yk=Lxk. xk is generated by Algorithm 1.We denote the error angle betweenyk and y* by ϕk=∠(yk,y*). Then, the current iteration yk can be decomposed as
(27)yk=y*cosϕk+uksinϕk,
where ukTy*=0 and ∥uk∥=∥y*∥=1. From (27), we obtain
(28)λmin-λk=λmin-ykTCyk=λmin-(y*cosϕk+uksinϕk)T×C(y*cosϕk+uksinϕk)=λmin-λmincos2ϕk-(ukTCuk)sin2ϕk=(λmin-(ukTCuk))sin2ϕk.According to the previous discussion, we give the following local convergence property for Algorithm1.Theorem 2.
Assume that the sequence(xk,λk) is generated by Algorithm 1. Then, the sequence of error angles ϕk tends to zero superlinearly.Proof.
From (21) and Cy*=λ*y*, we have
(29)λk-λ*=ykT(C-λ*I)yk=(yk-y*)T(C-λ*I)(yk-y*).
It gives
(30)|λk-λ*|≤|λmax-λ*|∥yk-y*∥2.
Thus, we choose a small enough positive ε1>0, such as ε1=3/4|λ*-μ2|, and a large enough k to satisfy
(31)|λk-μ2|>ε1,|λmax-λ*|∥yk-y*∥2<12ε1.
Then, combining inequality (30), we have
(32)|λk-λ*|≤12ε1.
Consequently, we obtain(33)qk(Δtk)∶=1/Δtk+|λ*-λk|1/Δtk+|μ2-λk|≤1/Δtk+1/2ε11/Δtk+ε1≤1/Δtmin+1/2ε11/Δtmin+ε1<1.
We can verify that there exists a positive constantεt such that all the time steps Δtk are bounded below by Δtmin; that is, Δtk≥Δtmin. Actually, from (13), we know that ∥xk+1P∥ will be close to ∥xk∥, when Δtk is small enough. Consequently, for small enough Δtk, we have
(34)|xk+1PBxk+1-1|=|xk+1PBxk+1-xkBxk|=|(xk+1P-xk)TB(xk+1P-xk)|<η1.
According to the time-stepping strategy (23), Δtk+1 will be enlarged. Therefore, Δtk are bounded below by a positive number Δtmin.
On the other hand, from (13) and (27), we obtain
(35)yk+1P=(I+Δtk(C-λkI))-1yk=(I+Δtk(C-λkI))-1×(y*cosϕk+uksinϕk)=y*cosϕk(1+Δtk(λmin-λk))+(I+Δtk(C-λkI))-1uksinϕk=y*cosϕk(1+Δtk(λmin-λk))+uk+1sinϕk∥(I+Δtk(C-λkI))-1uk∥,
where uk+1=(I+Δtk(C-λkI))-1uk/∥(I+Δtk(C-λkI))-1uk∥ and uk+1Ty*=0. Noticing yk+1=yk+1P/∥yk+1P∥ and yk+1=y*cosϕk+1+uk+1sinϕk+1, from (35), we obtain
(36)tanϕk+1=tanϕk∥(I+Δtk(C-λkI))-1uk∥×(1+Δtk(λmin-λk)).
Sinceuk is orthogonal to y* and ∥uk∥=1,  uk can be decomposed as uk=α2z2+⋯+αnzn, where α22+α32+⋯+αn2=1. Let μi(i=1,2,…,n) be the eigenvalues of C=L-TAL-1 in ascending order. Thus, from the spectral decomposition of C, for large enough k, we obtain
(37)∥(I+Δtk(C-λkI))-1uk∥2=∥∑i=2n‍(1+Δtk(μi-λk))-1αizi∥2=∑i=2n‍(1+Δtk(μi-λk))-2αi2≤∑i=2n‍(1+Δtk(μ2-λk))-2αi2=(1+Δtk(μ2-λk))-2.
It gives
(38)∥(I+Δtk(C-λkI))-1uk∥≤1(1+Δtk(μ2-λk)).
From (36) and (38), for large enough k, we have
(39)|tanϕk+1|≤|tanϕk|(1/Δtk+|λ*-λk|)(1/Δtk+μ2-λk).
Combining (33) and (39), we know that ϕk converges linearly to zero; namely, yk converges linearly to y*, when we choose the initial k such that (yk,λk) is close enough (y*,λ*).
Since(yk,λk) converges linearly to (y*,λ*), according to the time-stepping strategy (23), we obtain that the time step Δtk will tend to infinity. Consequently, from inequality (39), we know that ϕk superlinearly converges to 0.Remark 3.
From (28) and inequality (39), we have
(40)|tanϕk+1|≤|tanϕk|(1/Δtk+|λmin-λk|)(1/Δtk+μ2-λk)=|tanϕk|(1/Δtk+|(λmin-ukTAuk)sin2ϕk|)(1/Δtk+μ2-λk).
Consequently, the error angles ϕk tend to zero cubically, if we adopt a suitable time-stepping strategy such as 1/Δtk≈sin2ϕk for large enough k.Remark 4.
The convergence rate of the Rayleigh quotient iteration is cubic; however, its iterate point sequence{xk} is not guaranteed to converge to the smallest generalized eigenvalue of (A,B) [17]. The convergence rate of Algorithm 1 is less than the convergence rate of the Rayleigh quotient iteration, but the iterate point sequence {xk} of Algorithm 1 follows the trajectory of the special dynamical system of DAEs (5)-(6), and consequently it converges to a stationary point x* of the system of DAEs (5)-(6), that is, one of generalized eigenvectors associated with the smallest generalized eigenvalue λ* of (A,B).


## 5. Numerical Experiments

We give some numerical experiments for Algorithm1 (which is denoted by the PCM method), in comparison with the recent continuous projected method by Gao et al. (which is denoted by the GGL method [4]) the restarted Arnoldi method (the EIGS method [18]), and the Jacobi-Davidson method (JDQZ [19, 20]). For Algorithm 1, we pick η1=0.25, η2=0.75, γ1=2, and γ2=0.5 and compute an initial time-step length Δt0 as follows:
(41)Δt0=min{0.01,1∥Ax0-λ0Bx0∥2}.For the GGL method, we use the solver ODE45 [12] for the continuous projected dynamical method and set RelTol=10-6 and AbsTol=10-9, which are suggested in [4]. All of our tests were run in Matlab 2009a environment on a Lenovo Thinkpad laptop X200 with 2.4 GHz processor. The compared methods for the test problems are terminated when the condition
(42)∥Axk-λkBxk∥∞≤10-6
is satisfied. The descriptions of test problems and the numerical results are presented by the following.Example 5 (see [4, 21]).
This test problem has four generalized eigenvalues. The first three generalized eigenvalues are 0, 0, and 0, and the other is2. The stiffness and mass matrices are given as follows:
(43)A=[10-100000-10100000],B=[2010020110200102].
It is not difficult to verify that matrixB is positive definite. We choose an initial point x0=[0.2,0.6,-0.8,0.6]T.Example 6 (see [4, 21]).
The stiffness and mass matrices are given as follows:(44)A=[11k1+k2-k2-k2k2+k3⋱⋱⋱⋱⋱kn-1+kn-kn-kn11kn+kn+1],
and B=diag(m1,m2,…,mn) with n=20, ki=4+i, and mi=35-i for i=1,2,…,20.Example 7 (see [22]).
This test problem is from fluid flow generalized eigenvalues (see the Bcsstruc1 set in [22]).
Consider the following(45)Kx=λMx.
Here, the mass matrix M (bcsstm13.mtx) is a real symmetric positive semidefinite matrix; and the stiffness matrix K (bcsstk13.mtx) is a real positive definite matrix. Therefore, its smallest generalized eigenvalue is zero. In order to test the effect of algorithms, we apply a sift to M as follows:
(46)B=M+2K.
Then, the smallest generalized eigenvalue of (K,B) is 1/2.Example 8.
In order to test the performance of these two methods for large-scale problems, we construct a large-sparse matrixA as follows:  (1)
letA1=sparse(diag(ones(n,1)));(2)
letA2=sparse(diag(ones(n-1,1),1));(3)
letA3=sparse(diag(ones(n-1,1),-1));(4)
letA=2A1+A2+A3;(5)
letB=4A1+A2+A3.
Here, sparse.m, diag.m, and ones.m are the Matlab functions, andn is the dimension of A. In our numerical experiments, we set n=1000 for this test problem.Numerical results for these four test problems are reported in Table1. Iter and CPU designate the number of iterations and the computed time (given in seconds), respectively. From the computed results for Examples 5 and 8 in Table 1, we see that the proposed method (the PCM method) manages to find a smaller minimum eigenvalue than the compared method (the GGL method). Furthermore, the computed time of the PCM method is less than the GGL method and the JDQZ method for these test problems. For Example 5, the EIGS method fails since it gives the message: “(A-sigma*B) is singular. The shift is an eigenvalue. Try to use some other shift please,” when we execute it.Table 1
Comparison with different methods.
Examples
PCM
GGL
EIGS
JDQZ
λ
k
Iter (CPU)
λ
k
Iter (CPU)
CPU
CPU
Example5
1.0
e
-
13
12 (0.021)
4.5
e
-
11
77 (0.178)
FAIL
FAIL
Example6
0.01337
18 (0.004)
0.01344
3389 (0.36)
0.01
0.49
Example7
0.4998
49 (2.556)
0.5
33 (6.643)
0.24
14.17
Example8
1.1
e
-
5
28 (0.024)
0.0122
514 (169.3)
0.02
0.686In order to record the convergence history, we draw the iteration trajectory of the smallest eigenvalue computed by the PCM method and by the GGL method for every test problem in Figures1, 2, 3, and 4. From those figures, we find that Algorithm 1 (the PCM method) achieves the steady state of the dynamical system faster than the GGL method.Figure 1
Numerical results of Example5 with PCM and GGL methods.Figure 2
Numerical results of Example6 with PCM and GGL methods.Figure 3
Numerical results of Example7 with PCM and GGL methods.Figure 4
Numerical results of Example8 with PCM and GGL methods.In order to test the sensitivity of Algorithm1, we compute the test problems with Algorithm 1 for different parameters γ1, γ2, η1, and η2. The robustness of Algorithm 1 is mainly determined by the parameters η1 and η2 and the parameters η1 and η2 affect its efficiency. Thus, we fix parameters γ1=2 and γ2=1/2. We choose four groups of parameters η1 and η2 as following:(I)
η1=0.25,  η2=0.75;(II)
η1=0.1, η2=0.9;(III)
η1=0.3, η2=0.7;(IV)
η1=0.4, η2=0.8.The numerical results are reported in Table2. From Table 2, we find that Algorithm 1 is not sensitive for different parameters η1 around 0.25 and η2 around 0.75. Thus, we set the default parameters η1=0.25 and η2=0.75 in Algorithm 1.Table 2
Numerical results for different parameters of Algorithm1.
Examples
(I)
(II)
(III)
(IV)
Iter (CPU)
Iter (CPU)
Iter (CPU)
Iter (CPU)
Example5
12
12
12
12
Example6
19
19
19
17
Example7
49
49
50
45
Example8
26
31
22
22


## 6. Conclusions

This paper discusses the connection between the smallest generalized eigenvalue problem and the continuous dynamical system of DAEs. The trajectory of this special dynamical system is followed by a new prediction-correction method (Algorithm1), which is derived from the implicit Euler formula and an analogous trust-region technique. Consequently, an equilibrium point x* of this continuous dynamical system is obtained. This equilibrium point x* is also a generalized eigenvector associated with the smallest generalized eigenvalue. The local superlinear property for the new prediction-correction method is presented. Numerical experiments indicate that Algorithm 1 is promising for the smallest generalized eigenvalue problem. Another interesting issue is whether or not the second-order, pseudo-transient method [23] obtains the superior numerical property for this special dynamical system of DAEs (5)-(6). We would like to consider this issue in our future works.



---

*Source: 845459-2013-09-11.xml*