# Behavior of the Correction Equations in the Jacobi–Davidson Method

Authors: Yuan Kong; Yong Fang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The Jacobi–Davidson iteration method is efficient for computing several eigenpairs of Hermitian matrices. Although the involved correction equation in the Jacobi–Davidson method has many developed variants, the behaviors of them are not clear for us. In this paper, we aim to explore, theoretically, the convergence property of the Jacobi–Davidson method influenced by different types of correction equations. As a by-product, we derive the optimal expansion vector, which imposed a shift-and-invert transform on a vector located in the prescribed subspace, to expand the current subspace.


---

## Body

## 1. Introduction

The complexn-dimensional vector space and n-by-n matrix space are denoted by ℂn and ℂn×n, respectively; and these notations can be straightforwardly carried over to real space. Many iterative methods for computing a few extreme eigenpairs of the standard eigenvalue problem(1)Ax=λx,withx=1,where A∈ℂn×n is a large, sparse Hermitian matrix, λ∈ℝ is an eigenvalue, and x∈ℂn is the associated eigenvector of matrix A, have been developed, such as the Krylov subspace methods [1, 2], the locally optimal block preconditioned conjugate gradient method [3], the Davidson method [4], and their variants [5–8]. Here, we use ⋅ to denote the Euclidean norm of a vector, which will be defined in the next section.In the Jacobi–Davidson iteration method, the projection onto the orthogonal complement of the current approximation elegantly avoids the singularity of the shifted matrix and ensures the positive definite of the resulting coefficient matrix when the current approximation is near to the desired eigenvalue. In [9], the authors proposed several alternative variants of the correction equations and presented some predictions reasonably on the convergence property influenced by these correction equations. Also, in [10], the authors explored the existence and uniqueness on solutions of different correction equations. However, to our best knowledge, there is no analysis to compare the convergence affected by different types of correction equations [11, 12]. In this paper, we aim to analyse, theoretically, the convergence property of the Jacobi–Davidson method influenced by different types of correction equations.The organization of this paper is as follows. In Section2, we firstly introduce some notations and present some preliminaries. In Section 3, we then give the main convergence results. And, finally, in Section 5, we end the paper with some concluding remarks.


## 2. Preliminaries

For convenience of statements, we firstly introduce some notations. In this paper,In is used to represent the n×n identity matrix, and, for short, it can be simplified as I. And, ei with a suitable dimension is used to denote the i-th column of I. We use ⋅,⋅ to denote the inner product defined in the complex vector space ℂn, and the corresponding induced Euclidean norm of either a vector or a matrix is denoted by ⋅. For two given vectors v,w∈ℂn, cosϕ=v,w/vw indicates the angle between v and w, and, thus, the distance between them, with v=1 and w=1, can be expressed, equivalently, as sin∠v,w=I−vv∗w. In addition, for a given matrix G∈ℂn×n, GT and G∗ denote its transpose and conjugate transpose, respectively, and this notation can be straightforwardly carried over to vectors.The eigenpairs of the Hermitian eigenvalue problem are denoted byλi,xii=1n, with xi, i=1,2,…,n, being orthonormal, (1). For a given matrix pair A,B and a nonzero vector v∈ℂn, the function,(2)RA,B;v=v∗Avv∗Bv,is defined as the generalized Rayleigh quotient of v associated with the matrix pair A,B.In [6], Sleijpen and Van der Vorst proposed the efficient Jacobi–Davidson iteration method for computing a few eigenpairs of the Hermitian eigenvalue problem (1). In each iteration step, say, the k-th, of the Jacobi–Davidson method, the so-called correction equation,(3)I−u1u1∗A−θII−u1u1∗t=−r,fort⊥u1,with r=A−θIu1 being the residual with respect to the current unit approximation u1 and its Rayleigh quotient θ=u1∗Au1, requires to be solved to expand the current projection subspace.It should be noticed that the projection onto the orthogonal complement ofu1 guarantees positive definite and well-conditioning of the coefficient matrix of the correction equation in (3). In addition, the Jacobi–Davidson iteration method attains, elegantly, the cubic convergence rate locally [6, 13–15].In [9], Genseberger and Sleijpen proposed a more restrictive alternative correction equation:(4)I−UU∗A−θII−UU∗t=−r,fort⊥U,where U is an n-by-m matrix of which the columns expand an orthonormal basis of the current projection subspace, say, subspace U.Although the projection in (4) is more complicated and expensive than that in (3), projecting onto the orthogonal complement of the current approximate subspace may further improve the conditioning of the correction equation, especially in the case of clustering eigenvalues. Obviously, if U is chosen as U=u1, with u1 being the current Ritz vector, the correction equation in (4) is reduced to that in (3).Generally speaking, the solutions of the correction equations in (3) and (4) do not lead to the same expansion of the projection subspace, and, as emphasized in [9], it is unclear which of them is more effective in actual computations. Therefore, different convergence behaviors of the Jacobi–Davidson method are expected.In the Jacobi–Davidson iteration method, once we obtain the current search subspaceU, we firstly project matrix A onto it and obtain the corresponding projected matrix T=U∗AU; we then compute the eigenpairs θi,sii=1m of T with smaller size, that is, Tsi=θisi with si=1; and, finally, use the Ritz pairs θi,Usii=1m to approximate the desired eigenpairs of matrix A, see [6] for more details. Note that the Ritz vectors ui=Usi, i=1,2,…,m, can also span subspace U, and, thus, in the subsequent discussions, we assign U=u1,u2,⋯,um in the correction equation in (4).


## 3. The Main Results

In this section, we exhibit the main result which depicts different convergence behaviors of the Jacobi–Davidson method affected by different types of correction equations (3) and (4).If the correction equations (3) and (4) are solved exactly, we, frequently, obtain the solutions:(5)to=−u1+A−θI−1u1u1∗A−θI−1u1,tg=−u1+A−θI−1UU∗A−θI−1U−1U∗u1,respectively. Thus, expanding the current search subspace U by to and tg in the above equations is equivalent to(6)to=A−θI−1Uyo,tg=A−θI−1Uyg,where yo=e1, yg=U∗A−θI−1U−1e1, and e1=1,0,⋯,0∈ℝm, as u1 is included in the current search subspace.Next, we first present a result which depicts the distancesin∠x1,ty between ty=A−θI−1Uy and the desired eigenvector x1.Theorem 1.
ConsiderU=u1,u2,⋯,um whose columns ui, i=1,2,…,m, are Ritz vectors of A, and θ=u1∗Au1, the Rayleigh quotient of u1. Assume that ui admits the orthogonal direct sum decomposition:(7)ui=c1ix1+c2ix2+…+cnixn,with∑k=1ncki2=1,i=1,2,…,m,and let ty=A−θI−1Uy, then ye=C∗e1 solves the following minimum problem:(8)miny∈ℂm×1sin∠x1,ty=miny∈ℂm×1I−x1x1∗A−θI−1UyA−θI−1Uy,where C=cij∈ℂn×m with its elements being cij and e1=1,0,⋯,0T∈ℝn.Proof.
Based on the orthogonal direct sum decomposition (7), the matrix U can be represented as U=XC with X=x1,x2,…,xn, which indicates that(9)A−θI−1U=A−θI−1XC=XΣ−1C,where Σ=diagλ1−θ,λ2−θ,…,λn−θ. Then, it follows from I−x1x1∗X=0,x2,…,xn that(10)I−x1x1∗A−θI−1Uy2=y∗C∗Σ−10x2∗⋮xn∗0,x2,…,xnΣ−1Cy=y∗C∗00⋯00λ2−θ−2⋯0⋮⋮⋱⋮00⋯λn−θ−2Cy,A−θI−1Uy2=y∗C∗Σ−2Cy.
According to the property of the generalized Rayleigh quotient, we know that vectory such that Cy=e1 minimizes(11)sin∠x1,ty=I−x1x1∗A−θI−1Uy2A−θI−1Uy2,with the minimum being zero, and, thus, we obtain the least squares solution ye=C∗e1 uniquely.
In [16], the author considered the problem of optimal expansion of subspaces for eigenvector approximations. Specifically, starting from a prescribed subspace, the author derived the optimal resulting vector, obtained from multiplying matrix A with a vector located in the prescribed subspace, to expand the current prescribed subspace. Likewise, Theorem 1 provides the optimal vector, obtained by imposing a shift-and-invert transform A−θI−1 on a vector located in the prescribed subspace, to expand the current subspace.Theorem 2.
Assume that the current search subspace spanned by U obtained by the Jacobi–Davidson method satisfies(12)C1Tyg>c11yg,and then, we have the following estimate(13)cos∠yg,ye>cos∠yo,ye,where ye, yo, and yg are defined as above and C1T=c11,c12,…,c1m is the first row of matrix C.Proof.
It follows from the orthogonal direct sum decomposition in (7) that(14)yg=U∗A−θI−1U−1e1=C∗Σ−1C−1e1.
As a result, under assumption (12) and by straightforward computations, we obtain(15)cos∠yg,ye>cos∠yo,ye,which completes the proof.Theorem 2 indicates that vectoryg approximates the least squares solution ye better than that of yo under the assumption in (12), which straightforwardly results in that the exact solution of (4) approximates the desired eigenvector x1 better than that of the exact solution of (3).


## 4. Numerical Results

In this section, we use one simple example to examine the numerical behavior of the Jacobi–Davidson method affected by different type of correction equations. We simplify the Jacobi–Davidson method equipped with (3) as the JD method and that equipped with (4) as the JDG method.We would like to compare the two methods JD and JDG in aspects of the number of iteration steps (“IT”) and the computing time in seconds (“CPU”). In addition, for the last approximationθko, we use the absolute error (“ERR”) of the approximation to the desired eigenvalue, namely, ERR=θko−λ1, to check whether the two methods converge to a same target eigenvalue. The two methods are terminated once their residual norms could achieve the stopping criterion rk/rk≤10−8, where rk=A−θkIxk is the residual, xk and θk are the approximate eigenvector, and the Rayleigh quotient of xk, respectively.Example 1.
[1, 5] Let us consider the following partial differential equation:(16)−∂∂xμx,y∂u∂x−∂∂yνx,y∂u∂y=λu,defined on the domain of 0,1×0,1 imposed with the Dirichlet boundary conditions. We could acquire the eigenvalue problem (1) by discretization of the partial differential equation (16) by five-point finite difference scheme on an m×m grid. Here, μx,y and νx,y are chosen with μx,y=νx,y=−exy/1+x+y, and the dimension of the resulting matrix is n=m2.
We use Table1 to list the numerical results of Example 1. From this table, we observe that the JDV method behaves more efficient and robust than the JD method, and that the JDV method outperforms the JD method in both terms of iteration steps and CPU time, although the two methods converge to the desired eigenpair globally.Table 1
Numerical results for Example 1.
Method
Index
m
2
5
2
6
2
7
2
8
JD
IT
85
94
166
175
CPU
1.315
8.124
60.881
249.865
ERR
4.62E − 11
1.36E − 10
9.18E − 10
9.31E − 10
JD
IT
71
85
91
94
CPU
0.936
6.349
26.837
120.923
ERR
1.18E − 11
1.56E − 10
6.11E − 10
5.53E − 09


## 5. Concluding Remarks

In this paper, we analyse, theoretically, the behaviors of the Jacobi–Davidson method affected by correction equations (3) and (4), respectively. Admittedly, the comparisons may be impractically useful for unknowns of the eigenvectors, and, thus, presenting their comparisons in actual computations may be a challenging and interesting problem.



---

*Source: 5169362-2019-08-05.xml*