Abstract and Applied Analysis

Volume 2013 (2013), Article ID 845459, 8 pages

http://dx.doi.org/10.1155/2013/845459

Research Article

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

Xin-long Luo, Jia-ru Lin, and Wei-ling Wu

School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, P.O. Box 101, Beijing 100876, China

Received 12 May 2013; Revised 18 July 2013; Accepted 1 August 2013

Academic Editor: Chang-Hua Lien

Copyright © 2013 Xin-long Luo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

where and are real symmetric matrices and is positive definite. Problem (1) can be converted to an equivalent equality-constrained optimization problem [4]:

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 :

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

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:

on the generalized spherical surface,

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 along the trajectory of (5) and (6), we obtain

Consequently, we obtain

Substituting (8) into (5), we also obtain the following generalized gradient flow:

Conversely, if is a solution of (9), then ; that is, satisfies , 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 . 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

Furthermore, the solution of (9) converges to an eigenvector associated with the smallest generalized eigenvalue of , for almost all initial condition , where satisfies , 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 [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:

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 with in (11) to obtain the predicted point as follows:

Since the predicted point moves away from the elliptic spherical surface, we project this predicted point onto the unit elliptic spherical surface, so that the iterated point satisfies the algebraic constraint (6); namely, we find the shortest distance between and the unit elliptic spherical surface. We achieve this aim by solving the following equality-constrained problem:

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:

Since , we obtain the suboptimal solution by solving the following problem:

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 . Then, the suboptimal solution point of problem (14) is obtained as

Now, we use the correction point to update , so that the steady-state equation is satisfied as possible. Namely, we find the minimum point of the following problem:

Since matrix is positive, it can be decomposed as . Let . Then, the previous minimum problem is equivalent to the following problem:

Using the property , we obtain the suboptimal solution of problem (18) by solving the following problem:

It is not difficult to obtain the minimum solution of the previous problem as

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 is also near the elliptic spherical surface ; otherwise, the time step is reduced. In order to measure the distance between and the unit elliptic spherical surface, we construct the following approximate model:

Then,

where constants , , , and satisfy

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

Compute and the residual .

Step  2. Repeat for until :

(a)solve (13) to obtain the predicted point ;(b)evaluate (17) to obtain the corrected point ;(c)compute (21) to update the Lagrange multiplier ;(d)compute the residual

(e)update the time step according to the time-stepping scheme (22)-(23);(f)accept the trial vector. If , let , and ;(g)try again starting at (a) with .

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 . 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

where and . From (27), we obtain

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

It gives

Thus, we choose a small enough positive , such as , and a large enough to satisfy

Then, combining inequality (30), we have

Consequently, we obtain

We can verify that there exists a positive constant such that all the time steps are bounded below by ; that is, . Actually, from (13), we know that will be close to , when is small enough. Consequently, for small enough , we have

According to the time-stepping strategy (23), will be enlarged. Therefore, are bounded below by a positive number .

On the other hand, from (13) and (27), we obtain

where and . Noticing and , from (35), we obtain

Since is orthogonal to and ,   can be decomposed as , where . Let be the eigenvalues of in ascending order. Thus, from the spectral decomposition of , for large enough , we obtain

It gives

From (36) and (38), for large enough , we have

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

Consequently, the error angles tend to zero cubically, if we adopt a suitable time-stepping strategy such as 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 , , , and and compute an initial time-step length as follows:

For the GGL method, we use the solver ODE45 [12] for the continuous projected dynamical method and set and , 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

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 . The stiffness and mass matrices are given as follows:

It is not difficult to verify that matrix is positive definite. We choose an initial point .

Example 6 (see [4, 21]). The stiffness and mass matrices are given as follows:

and with , , and for .

Example 7 (see [22]). This test problem is from fluid flow generalized eigenvalues (see the Bcsstruc1 set in [22]).

Consider the following

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:

Then, the smallest generalized eigenvalue of is .

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 ; (2)let ; (3)let ; (4)let ; (5)let .

Here, sparse.m, diag.m, and ones.m are the Matlab functions, and is the dimension of . In our numerical experiments, we set 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.

Examples

PCM

GGL
EIGS
JDQZ

Iter (CPU)

Iter (CPU)
CPU
CPU

Example 5

12 (0.021)

77 (0.178)
FAIL
FAIL

Example 6
0.01337
18 (0.004)
0.01344
3389 (0.36)
0.01
0.49

Example 7
0.4998
49 (2.556)
0.5
33 (6.643)
0.24
14.17

Example 8

28 (0.024)
0.0122
514 (169.3)
0.02
0.686

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 and . We choose four groups of parameters and as following:

(I), ;(II), ;(III), ;(IV), .

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

Table 2: Numerical results for different parameters of Algorithm 1.

Examples
(I)
(II)
(III)
(IV)

Iter (CPU)
Iter (CPU)
Iter (CPU)
Iter (CPU)

Example 5
12
12
12
12

Example 6
19
19
19
17

Example 7
49
49
50
45

Example 8
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 (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.

References

T. Auckenthaler, T. Auckenthaler, H.-J. Bungartz, et al., “Parallel solution of partial symmetric eigenvalue problems from electronic structure calculations,” Parallel Computing, vol. 37, no. 11, pp. 783–794, 2011.

Y. Saad, Numerical Methods for Large Eigenvalue Problems, Algorithms and Architectures for Advanced Scientific Computing, Manchester University Press, Manchester, England, 1992.

G. H. Golub and L.-Z. Liao, “Continuous methods for extreme and interior eigenvalue problems,” Linear Algebra and its Applications, vol. 415, no. 1, pp. 31–51, 2006.

X.-B. Gao, G. H. Golub, and L.-Z. Liao, “Continuous methods for symmetric generalized eigenvalue problems,” Linear Algebra and its Applications, vol. 428, no. 2-3, pp. 676–696, 2008.

U. Helmke and J. B. Moore, Optimization and Dynamical Systems, Springer, 2nd edition, 1996.

R. Mahony and P.-A. Absil, “The continuous-time Rayleigh quotient flow on the sphere,” Linear Algebra and its Applications, vol. 368, pp. 343–357, 2003.

C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, NJ, USA, 1971.

B. He, X. Yuan, and J. J. Z. Zhang, “Comparison of two kinds of prediction-correction methods for monotone variational inequalities,” Computational Optimization and Applications, vol. 27, no. 3, pp. 247–267, 2004.

B. He and X. Yuan, “On the convergence rate of the Douglas-Rachford alternating direction method,” SIAM Journal on Numerical Analysis, vol. 50, no. 2, pp. 700–709, 2012.

J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research, Springer, New York, NY, USA, 1999.

U. M. Ascher and L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1998.

L. F. Shampine and M. W. Reichelt, “The Matlab ODE suite,” SIAM Journal on Scientific Computing, vol. 18, no. 1, pp. 1–22, 1997.

S.-T. Liu and X.-L. Luo, “A method based on Rayleigh quotient gradient flow for extreme and interior eigenvalue problems,” Linear Algebra and its Applications, vol. 432, no. 7, pp. 1851–1863, 2010.

M. K. Ng, F. Wang, and X. Yuan, “Inexact alternating direction methods for image recovery,” SIAM Journal on Scientific Computing, vol. 33, no. 4, pp. 1643–1668, 2011.

C. T. Kelley, L.-Z. Liao, L. Qi, M. T. Chu, J. P. Reese, and C. Winton, “Projected pseudotransient continuation,” SIAM Journal on Numerical Analysis, vol. 46, no. 6, pp. 3071–3083, 2008.

X. L. Luo, “A dynamical method of DAEs for the smallest eigenvalue problem,” Journal of Computational Science, vol. 3, no. 3, pp. 113–119, 2012.

B. N. Parlett, The Symmetric Eigenvalue Problem, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1998, Corrected reprint of the 1980 original.

MATLAB R2009a, http://www.mathworks.com

D. R. Fokkema, G. L. G. Sleijpen, and H. A. Van der Vorst, “Jacobi-Davidson style QR and QZ algorithms for the reduction of matrix pencils,” SIAM Journal on Scientific Computing, vol. 20, no. 1, pp. 94–125, 1998.

http://www.staff.science.uu.nl/?sleij101/index.html

T. Zhang, K. H. Law, and G. H. Golub, “On the homotopy method for perturbed symmetric generalized eigenvalue problems,” SIAM Journal on Scientific Computing, vol. 19, no. 5, pp. 1625–1645, 1998.

Matrix Market, http://math.nist.gov/MatrixMarket

X.-l. Luo, “A second-order pseudo-transient method for steady-state problems,” Applied Mathematics and Computation, vol. 216, no. 6, pp. 1752–1762, 2010.