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

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.

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 (𝐴,𝐵). 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 (𝐴,𝐵) is to find a scalar 𝜆 and nonzero vector 𝑥∈ℜ^𝑛 to satisfy 𝐴𝑥=𝜆𝐵𝑥, (1) where 𝐴 and 𝐵 are 𝑛×𝑛 real symmetric matrices and 𝐵 is positive definite. Problem (1) can be converted to an equivalent equality-constrained optimization problem [4]: min_{𝑥∈ℜ𝑛}𝑥^𝑇𝐴𝑥,s.t.𝑥^𝑇𝐵𝑥=1. (2) If 𝑥^∗ 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 (𝑥^∗,𝜆^∗): 𝐴𝑥^∗−𝜆^∗𝐵𝑥^∗𝑥=0,_∗𝑇𝐵𝑥^∗=1. (3) Therefore, the solution 𝑥^∗ of problem (2) is a generalized eigenvector of (𝐴,𝐵).

Property 1 (Gao et al., 2008 [4]). If 𝑥^∗ is an optimal solution of problem (2) and 𝜆^∗ is the corresponding Lagrange multiplier, then they satisfy the second-order optimality condition 𝐴−𝜆^∗𝐵≽0. (4) Furthermore, a local minimizer of (2) is also its global minimizer. If 𝑥^∗ 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) on the generalized spherical surface, 𝑥^𝑇𝐵𝑥=1. (6) 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 constraint 𝑥^𝑇𝐵𝑥=1 along the trajectory (𝑥(𝑡),𝜆(𝑡)) of (5) and (6), we obtain 𝑑𝑥^𝑇𝐵𝑥=∇𝑑𝑡^𝑥𝑥^𝑇𝐵𝑥^𝑇𝑑𝑥𝑑𝑡=2𝑥^𝑇𝐵(−𝐴𝑥+𝜆𝐵𝑥)=2−𝑥^𝑇𝐵𝐴𝑥+𝜆𝑥^𝑇𝐵²𝑥=0. (7) Consequently, we obtain 𝑥𝜆=^𝑇𝐵𝐴𝑥𝑥^𝑇𝐵²𝑥. (8) Substituting (8) into (5), we also obtain the following generalized gradient flow: 𝑑𝑥(𝑡)𝑥𝑑𝑡=−𝐴−𝐵^𝑇𝐵𝐴𝑥𝑥^𝑇𝐵²𝑥𝑥,𝑥(0)^𝑇𝐵𝑥(0)=1. (9) Conversely, if 𝑥(𝑡) is a solution of (9), then 𝑑(𝑥^𝑇𝐵𝑥)/𝑑𝑡=0; that is, 𝑥(𝑡) satisfies 𝑥^𝑇𝐵𝑥=1, and consequently 𝑥(𝑡) 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 𝑥(𝑡) for any initial point 𝑥(0)=𝑥₀. We can verify that the objective function 𝑓(𝑥)=𝑥^𝑇𝐴𝑥 is monotonically decreasing along the trajectory 𝑥(𝑡) of (5)-(6). Actually, from (5)–(8), using the Cauchy-Schwartz inequality |𝑦^𝑇𝑧|≤‖𝑦‖‖𝑧‖, we have 𝑑𝑓(𝑥)=∇𝑑𝑡^𝑥𝑓(𝑥)^𝑇𝑑𝑥𝑑𝑡=2−𝑥^𝑇𝐴²𝑥𝑥+^𝑇𝐴𝐵𝑥𝑥^𝑇𝐵𝐴𝑥𝑥^𝑇𝐵²𝑥≤2𝑥^𝑇𝐵²𝑥−‖𝐴𝑥‖²‖𝐵𝑥‖²+(𝐴𝑥)^𝑇(𝐵𝑥)(𝐵𝑥)^𝑇(𝐴𝑥)≤0. (10) Furthermore, the solution 𝑥(𝑡) of (9) converges to an eigenvector 𝑥^∗ associated with the smallest generalized eigenvalue 𝜆^∗ of (𝐴,𝐵), for almost all initial condition 𝑥(0)=𝑥₀, where 𝑥₀ satisfies ‖𝑥₀‖=1, except for a set of measure zero, that is, the case of a multiple generalized eigenvalue 𝜆^∗ [5]. Therefore, by following the trajectory 𝑥(𝑡) of (5)-(6), we can compute a global minimum point 𝑥^∗ of problem (2), and consequently obtain the smallest generalized eigenvalue 𝜆^∗ of (𝐴,𝐵).

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 (𝐴,𝐵). 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 𝑑𝑥/𝑑𝑡=−𝜇𝑥, where 𝜇>0 [7], and consequently the implicit Euler method can take large steps at the steady-state phase so that the iteration sequence {𝑥_𝑘} converges rapidly to a stationary point 𝑥^∗ 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: 𝑥_𝑘+1=𝑥_𝑘−Δ𝑡_𝑘𝐴−𝜆_𝑘+1𝐵𝑥_𝑘+1,𝑥^𝑇_𝑘+1𝐵𝑥_𝑘+1=1, (11)(12) where Δ𝑡_𝑘 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 𝜆_𝑘+1 with 𝜆_𝑘 in (11) to obtain the predicted point 𝑥^𝑃_𝑘+1 as follows: 𝐼+Δ𝑡_𝑘𝐴−𝜆_𝑘𝐵𝑥^𝑃_𝑘+1=𝑥_𝑘. (13) Since the predicted point 𝑥^𝑃_𝑘+1 moves away from the elliptic spherical surface, we project this predicted point 𝑥^𝑃_𝑘+1 onto the unit elliptic spherical surface, so that the iterated point 𝑥_𝑘+1 satisfies the algebraic constraint (6); namely, we find the shortest distance between 𝑥^𝑃_𝑘+1 and the unit elliptic spherical surface. We achieve this aim by solving the following equality-constrained problem: min_{𝑥∈𝔑𝑛}𝑥−𝑥^𝑃_𝑘+1,s.t.𝑥^𝑇𝐵𝑥=1. (14)

Since the optimal solution of the previously problem is not apparently obtained, we turn to achieve its suboptimal solution. Since matrix 𝐵 is positive, it can be decomposed as 𝐵=𝐿^𝑇𝐿. Let 𝑦=𝐿𝑥. Then, problem (14) is equivalent to the following problem: min_{𝑦∈𝔑𝑛}𝐿⁻¹𝑦−𝐿𝑥^𝑃_𝑘+1,s.t.𝑦=1. (15) Since ‖𝐿⁻¹(𝑦−𝐿𝑥^𝑃_𝑘+1)‖≤‖𝐿⁻¹‖‖𝑦−𝐿𝑥^𝑃_𝑘+1‖, we obtain the suboptimal solution by solving the following problem: min_{𝑦∈𝔑𝑛}𝑦−𝐿𝑥^𝑃_𝑘+1,s.t.𝑦=1. (16) 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 𝑦_𝑘+1=𝐿𝑥^𝑇_𝑘+1/‖𝐿𝑥^𝑃_𝑘+1‖. Then, the suboptimal solution point of problem (14) is obtained as 𝑥_𝑘+1=𝐿⁻¹𝑦_𝑘+1=𝑥^𝑃_𝑘+1𝐿𝑥^𝑃_𝑘+1. (17)

Now, we use the correction point 𝑥_𝑘+1 to update 𝜆_𝑘+1, so that the steady-state equation 𝐴𝑥−𝜆𝐵𝑥=0 is satisfied as possible. Namely, we find the minimum point of the following problem: min_{𝜆∈𝔑𝑛}𝐴𝑥_𝑘+1−𝜆𝐵𝑥_𝑘+1. (18) Since matrix 𝐵 is positive, it can be decomposed as 𝐵=𝐿^𝑇𝐿. Let 𝑧_𝑘+1=𝐿𝑥_𝑘+1. Then, the previous minimum problem is equivalent to the following problem: min_{𝜆∈𝔑𝑛}𝐿^𝑇𝐿^−𝑇𝐴𝑥_𝑘+1−𝜆𝐿𝑥_𝑘+1. (19) Using the property ‖𝐿^𝑇(𝐿^−𝑇𝐴𝑥_𝑘+1−𝜆𝐿𝑥_𝑘+1)‖≤‖𝐿^𝑇‖‖(𝐿^−𝑇𝐴𝑥_𝑘+1−𝜆𝐿𝑥_𝑘+1)‖, we obtain the suboptimal solution of problem (18) by solving the following problem: min𝐿^−𝑇𝐴𝑥_𝑘+1−𝜆𝐿𝑥_𝑘+1. (20) It is not difficult to obtain the minimum solution of the previous problem as 𝜆_𝑘+1=𝑥^𝑇_𝑘+1𝐴𝑥_𝑘+1𝑥^𝑇_𝑘+1𝐵𝑥_𝑘+1, (21) which is the suboptimal solution of problem (18).

Another issue is how to choose the time step Δ𝑡_𝑘 at every iteration. We adopt an analogous trust-region technique to adaptively adjust the time step Δ𝑡_𝑘. Our intuition is that the time step Δ𝑡_𝑘 can be enlarged when the predicted point 𝑥^𝑃_𝑘+1 is also near the elliptic spherical surface 𝑥^𝑇𝐵𝑥=1; otherwise, the time step Δ𝑡_𝑘 is reduced. In order to measure the distance between 𝑥^𝑃_𝑘+1 and the unit elliptic spherical surface, we construct the following approximate model: 𝜌_𝑘=1−𝑥^𝑃_𝑘+1^𝑇𝐵𝑥^𝑃_𝑘+1^1/2. (22) Then, Δ𝑡_𝑘+1=𝛾₁Δ𝑡_𝑘,if𝜌_𝑘<𝜂₁,Δ𝑡_𝑘,if𝜂₁≤𝜌_𝑘𝜂₂,𝛾₂Δ𝑡_𝑘,if𝜌_𝑘≥𝜂₂, (23) where constants 𝜂₁, 𝜂₂, 𝛾₁, and 𝛾₂ satisfy 0<𝜂₁<𝜂₂<1,0<𝛾₂<1<𝛾₁. (24)

According to the previous discussion, we give the following discrete dynamical method for the smallest generalized eigenvalue of (𝐴,𝐵).

Algorithm 1. Prediction-correction method for the smallest generalized eigenvalue.
Step  1. Initialize the parameters. Specify a tolerated error TOL and an initial point 𝑥₀ to satisfy 𝑥^𝑇₀𝐵𝑥₀=1. (25) Compute 𝜆₀=𝑥^𝑇₀𝐴𝑥₀/𝑥^𝑇₀𝐵𝑥₀ and the residual 𝑟₀=‖𝐴𝑥₀−𝜆₀𝐵𝑥₀‖.
Step  2. Repeat for 𝑘=0,1,2,… until 𝑟_𝑘=‖𝐴𝑥_𝑘−𝜆_𝑘𝐵𝑥_𝑘‖≤TOL: (a)solve (13) to obtain the predicted point 𝑥^𝑃_𝑘+1;(b)evaluate (17) to obtain the corrected point 𝑥_𝑘+1;(c)compute (21) to update the Lagrange multiplier 𝜆_𝑘+1;(d)compute the residual 𝑟_𝑘+1=𝐴𝑥_𝑘+1−𝜆_𝑘+1𝐵𝑥_𝑘+1; (26)(e)update the time step Δ𝑡_𝑘+1 according to the time-stepping scheme (22)-(23);(f)accept the trial vector. If 𝜌_𝑘<𝜂₂, let 𝑥_𝑘=𝑥_𝑘+1, 𝜆_𝑘=𝜆_𝑘+1 and 𝑟_𝑘=𝑟_𝑘+1;(g)try again starting at (a) with 𝑘←𝑘+1.

4. Local Convergence Analysis

The sequence (𝑥_𝑘,𝜆_𝑘) 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 𝑥_𝑘 is 𝑂(Δ𝑡_𝑘). According to the standard analysis of numerical ordinary differential equations [7, 11], we know that 𝑥_𝑘 is close to 𝑥(𝑡_𝑘) when Δ𝑡_𝑘 is small enough in the finite time interval [0,𝑇]. Therefore, we only need to estimate the local convergence rate of 𝑥_𝑘 in the steady-state phase.

For the convenience of analysis, we decompose matrix 𝐵 as 𝐵=𝐿^𝑇𝐿 and denote 𝐶=𝐿^−𝑇𝐴𝐿⁻¹. Then, the smallest generalized eigenvalue 𝜆^∗ of (𝐴,𝐵) and associated with the eigenvector 𝑥^∗ is equivalent to the smallest eigenvalue of symmetric matrix 𝐶 and associated with the eigenvector 𝑦^∗=𝐿𝑥^∗, respectively. Therefore, we only need to consider the convergence property of (𝑦_𝑘,𝜆_𝑘), where 𝑦_𝑘=𝐿𝑥_𝑘. 𝑥_𝑘 is generated by Algorithm 1.

We denote the error angle between 𝑦_𝑘 and 𝑦^∗ by 𝜙_𝑘=∠(𝑦_𝑘,𝑦^∗). Then, the current iteration 𝑦_𝑘 can be decomposed as 𝑦_𝑘=𝑦^∗cos𝜙_𝑘+𝑢_𝑘sin𝜙_𝑘, (27) where 𝑢^𝑇_𝑘𝑦^∗=0 and ‖𝑢_𝑘‖=‖𝑦^∗‖=1. From (27), we obtain 𝜆_min−𝜆_𝑘=𝜆_min−𝑦^𝑇_𝑘𝐶𝑦_𝑘=𝜆_min−𝑦^∗cos𝜙_𝑘+𝑢_𝑘sin𝜙_𝑘^𝑇×𝐶𝑦^∗cos𝜙_𝑘+𝑢_𝑘sin𝜙_𝑘=𝜆_min−𝜆_mincos²𝜙_𝑘−𝑢^𝑇_𝑘𝐶𝑢_𝑘sin²𝜙_𝑘=𝜆_min−𝑢^𝑇_𝑘𝐶𝑢_𝑘sin²𝜙_𝑘. (28)

According to the previous discussion, we give the following local convergence property for Algorithm 1.

Theorem 2. Assume that the sequence (𝑥_𝑘,𝜆_𝑘) is generated by Algorithm 1. Then, the sequence of error angles 𝜙_𝑘 tends to zero superlinearly.

Proof. From (21) and 𝐶𝑦^∗=𝜆^∗𝑦^∗, we have 𝜆_𝑘−𝜆^∗=𝑦^𝑇_𝑘𝐶−𝜆^∗𝐼𝑦_𝑘=𝑦_𝑘−𝑦^∗𝑇𝐶−𝜆^∗𝐼𝑦_𝑘−𝑦^∗. (29) It gives 𝜆_𝑘−𝜆^∗≤𝜆_max−𝜆^∗𝑦^𝑘−𝑦^∗2. (30) Thus, we choose a small enough positive 𝜀₁>0, such as 𝜀₁=3/4|𝜆^∗−𝜇₂|, and a large enough 𝑘 to satisfy 𝜆_𝑘−𝜇₂>𝜀₁,𝜆_max−𝜆^∗𝑦^𝑘−𝑦^∗2<12𝜀₁. (31) Then, combining inequality (30), we have 𝜆_𝑘−𝜆^∗≤12𝜀₁. (32)
Consequently, we obtain 𝑞_𝑘Δ𝑡_𝑘∶=1/Δ𝑡_𝑘+𝜆^∗−𝜆_𝑘1/Δ𝑡_𝑘+𝜇₂−𝜆_𝑘≤1/Δ𝑡_𝑘+1/2𝜀₁1/Δ𝑡_𝑘+𝜀₁≤1/Δ𝑡_min+1/2𝜀₁1/Δ𝑡_min+𝜀₁<1. (33)
We can verify that there exists a positive constant 𝜀_𝑡 such that all the time steps Δ𝑡_𝑘 are bounded below by Δ𝑡_min; that is, Δ𝑡_𝑘≥Δ𝑡_min. Actually, from (13), we know that ‖𝑥^𝑃_𝑘+1‖ will be close to ‖𝑥_𝑘‖, when Δ𝑡_𝑘 is small enough. Consequently, for small enough Δ𝑡_𝑘, we have 𝑥^𝑃_𝑘+1𝐵𝑥_𝑘+1−1=𝑥^𝑃_𝑘+1𝐵𝑥_𝑘+1−𝑥_𝑘𝐵𝑥_𝑘=𝑥^𝑃_𝑘+1−𝑥_𝑘^𝑇𝐵𝑥^𝑃_𝑘+1−𝑥_𝑘<𝜂₁. (34) According to the time-stepping strategy (23), Δ𝑡_𝑘+1 will be enlarged. Therefore, Δ𝑡_𝑘 are bounded below by a positive number Δ𝑡_min.
On the other hand, from (13) and (27), we obtain 𝑦^𝑃_𝑘+1=𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑦_𝑘=𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹×𝑦^∗cos𝜙_𝑘+𝑢_𝑘sin𝜙_𝑘=𝑦^∗cos𝜙_𝑘1+Δ𝑡_𝑘𝜆_min−𝜆_𝑘+𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑢_𝑘sin𝜙_𝑘=𝑦^∗cos𝜙_𝑘1+Δ𝑡_𝑘𝜆_min−𝜆_𝑘+𝑢_𝑘+1sin𝜙_𝑘𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑢_𝑘, (35) where 𝑢_𝑘+1=(𝐼+Δ𝑡_𝑘(𝐶−𝜆_𝑘𝐼))⁻¹𝑢_𝑘/‖(𝐼+Δ𝑡_𝑘(𝐶−𝜆_𝑘𝐼))⁻¹𝑢_𝑘‖ and 𝑢^𝑇_𝑘+1𝑦^∗=0. Noticing 𝑦_𝑘+1=𝑦^𝑃_𝑘+1/‖𝑦^𝑃_𝑘+1‖ and 𝑦_𝑘+1=𝑦^∗cos𝜙_𝑘+1+𝑢_𝑘+1sin𝜙_𝑘+1, from (35), we obtain tan𝜙_𝑘+1=tan𝜙_𝑘𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑢_𝑘×1+Δ𝑡_𝑘𝜆_min−𝜆_𝑘. (36)
Since 𝑢_𝑘 is orthogonal to 𝑦^∗ and ‖𝑢_𝑘‖=1,  𝑢_𝑘 can be decomposed as 𝑢_𝑘=𝛼₂𝑧₂+⋯+𝛼_𝑛𝑧_𝑛, where 𝛼²₂+𝛼²₃+⋯+𝛼²_𝑛=1. Let 𝜇_𝑖(𝑖=1,2,…,𝑛) be the eigenvalues of 𝐶=𝐿^−𝑇𝐴𝐿⁻¹ in ascending order. Thus, from the spectral decomposition of 𝐶, for large enough 𝑘, we obtain 𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑢_𝑘²=^𝑛_𝑖=21+Δ𝑡_𝑘𝜇_𝑖−𝜆_𝑘⁻¹𝛼_𝑖𝑧_𝑖²=^𝑛_𝑖=21+Δ𝑡_𝑘𝜇_𝑖−𝜆_𝑘⁻²𝛼²_𝑖≤^𝑛_𝑖=21+Δ𝑡_𝑘𝜇₂−𝜆_𝑘⁻²𝛼²_𝑖=1+Δ𝑡_𝑘𝜇₂−𝜆_𝑘⁻². (37) It gives 𝐼+Δ𝑡_𝑘𝐶−𝜆_𝑘𝐼⁻¹𝑢_𝑘≤11+Δ𝑡_𝑘𝜇₂−𝜆_𝑘. (38) From (36) and (38), for large enough 𝑘, we have tan𝜙_𝑘+1≤tan𝜙_𝑘1/Δ𝑡_𝑘+𝜆^∗−𝜆_𝑘1/Δ𝑡_𝑘+𝜇₂−𝜆_𝑘. (39)
Combining (33) and (39), we know that 𝜙_𝑘 converges linearly to zero; namely, 𝑦_𝑘 converges linearly to 𝑦^∗, when we choose the initial 𝑘 such that (𝑦_𝑘,𝜆_𝑘) is close enough (𝑦^∗,𝜆^∗).
Since (𝑦_𝑘,𝜆_𝑘) converges linearly to (𝑦^∗,𝜆^∗), according to the time-stepping strategy (23), we obtain that the time step Δ𝑡_𝑘 will tend to infinity. Consequently, from inequality (39), we know that 𝜙_𝑘 superlinearly converges to 0.

Remark 3. From (28) and inequality (39), we have tan𝜙_𝑘+1≤tan𝜙_𝑘1/Δ𝑡_𝑘+𝜆_min−𝜆_𝑘1/Δ𝑡_𝑘+𝜇₂−𝜆_𝑘=tan𝜙_𝑘1/Δ𝑡_𝑘+𝜆_min−𝑢^𝑇_𝑘𝐴𝑢_𝑘sin²𝜙_𝑘1/Δ𝑡_𝑘+𝜇₂−𝜆_𝑘. (40) Consequently, the error angles 𝜙_𝑘 tend to zero cubically, if we adopt a suitable time-stepping strategy such as 1/Δ𝑡_𝑘≈sin²𝜙_𝑘 for large enough 𝑘.

Remark 4. The convergence rate of the Rayleigh quotient iteration is cubic; however, its iterate point sequence {𝑥_𝑘} is not guaranteed to converge to the smallest generalized eigenvalue of (𝐴,𝐵) [17]. The convergence rate of Algorithm 1 is less than the convergence rate of the Rayleigh quotient iteration, but the iterate point sequence {𝑥_𝑘} of Algorithm 1 follows the trajectory of the special dynamical system of DAEs (5)-(6), and consequently it converges to a stationary point 𝑥^∗ of the system of DAEs (5)-(6), that is, one of generalized eigenvectors associated with the smallest generalized eigenvalue 𝜆^∗ of (𝐴,𝐵).

5. Numerical Experiments

We give some numerical experiments for Algorithm 1 (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 𝜂₁=0.25, 𝜂₂=0.75, 𝛾₁=2, and 𝛾₂=0.5 and compute an initial time-step length Δ𝑡₀ as follows: Δ𝑡₀1=min0.01,𝐴𝑥₀−𝜆₀𝐵𝑥₀₂. (41)

For the GGL method, we use the solver ODE45 [12] for the continuous projected dynamical method and set RelTol=10⁻⁶ and AbsTol=10⁻⁹, 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 𝐴𝑥_𝑘−𝜆_𝑘𝐵𝑥_𝑘∞≤10⁻⁶ (42) 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 is 2. The stiffness and mass matrices are given as follows: ,.𝐴=10−100000−10100000𝐵=2010020110200102 (43)
It is not difficult to verify that matrix 𝐵 is positive definite. We choose an initial point 𝑥₀=[0.2,0.6,−0.8,0.6]^𝑇.

Example 6 (see [4, 21]). The stiffness and mass matrices are given as follows: 𝐴=11𝑘₁+𝑘₂−𝑘₂−𝑘₂𝑘₂+𝑘₃⋱⋱⋱⋱⋱𝑘_𝑛−1+𝑘_𝑛−𝑘_𝑛−𝑘_𝑛11𝑘_𝑛+𝑘_𝑛+1, (44) and 𝐵=diag(𝑚₁,𝑚₂,…,𝑚_𝑛) with 𝑛=20, 𝑘_𝑖=4+𝑖, and 𝑚_𝑖=35−𝑖 for 𝑖=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) Here, the mass matrix 𝑀 (bcsstm13.mtx) is a real symmetric positive semidefinite matrix; and the stiffness matrix 𝐾 (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 𝑀 as follows: 𝐵=𝑀+2𝐾. (46) Then, the smallest generalized eigenvalue of (𝐾,𝐵) is 1/2.

Example 8. In order to test the performance of these two methods for large-scale problems, we construct a large-sparse matrix 𝐴 as follows: (1)let 𝐴₁=sparse(diag(ones(𝑛,1))); (2)let 𝐴₂=sparse(diag(ones(𝑛−1,1),1)); (3)let 𝐴₃=sparse(diag(ones(𝑛−1,1),−1)); (4)let 𝐴=2𝐴₁+𝐴₂+𝐴₃; (5)let 𝐵=4𝐴₁+𝐴₂+𝐴₃.
Here, sparse.m, diag.m, and ones.m are the Matlab functions, and 𝑛 is the dimension of 𝐴. In our numerical experiments, we set 𝑛=1000 for this test problem.

Numerical results for these four test problems are reported in Table 1. 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: “(𝐴-sigma∗𝐵) is singular. The shift is an eigenvalue. Try to use some other shift please,” when we execute it.

Table 1 

Comparison with different methods.

In 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 Figures 1, 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 Example 5 with PCM and GGL methods.

Figure 2 

Numerical results of Example 6 with PCM and GGL methods.

Figure 3 

Numerical results of Example 7 with PCM and GGL methods.

Figure 4 

Numerical results of Example 8 with PCM and GGL methods.

In order to test the sensitivity of Algorithm 1, we compute the test problems with Algorithm 1 for different parameters 𝛾₁, 𝛾₂, 𝜂₁, and 𝜂₂. The robustness of Algorithm 1 is mainly determined by the parameters 𝜂₁ and 𝜂₂ and the parameters 𝜂₁ and 𝜂₂ affect its efficiency. Thus, we fix parameters 𝛾₁=2 and 𝛾₂=1/2. We choose four groups of parameters 𝜂₁ and 𝜂₂ as following:(I)𝜂₁=0.25, 𝜂₂=0.75;(II)𝜂₁=0.1, 𝜂₂=0.9;(III)𝜂₁=0.3, 𝜂₂=0.7;(IV)𝜂₁=0.4, 𝜂₂=0.8.

The numerical results are reported in Table 2. From Table 2, we find that Algorithm 1 is not sensitive for different parameters 𝜂₁ around 0.25 and 𝜂₂ around 0.75. Thus, we set the default parameters 𝜂₁=0.25 and 𝜂₂=0.75 in Algorithm 1.

Table 2 

Numerical results for different parameters of Algorithm 1.

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 (Algorithm 1), which is derived from the implicit Euler formula and an analogous trust-region technique. Consequently, an equilibrium point 𝑥^∗ of this continuous dynamical system is obtained. This equilibrium point 𝑥^∗ 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.

Acknowledgments

Xin-long Luo is grateful to Proffessor Li-Zhi Liao for fruitful discussions when he visited Hong Kong Baptist University in winter, 2011. This work was supported in part by Grant 2009CB320401 from the National Basic Research Program of China, Grant YBWL2011085 from Huawei Technologies Co., Ltd., and Grant YJCB2011003HI from the Innovation Research Program of Huawei Technologies Co., Ltd., Grant BUPT2009RC0118 from the Fundamental Research Funds for the Central Universities.