A modification of the Jacobi-Davidson method

TL;DR

This paper introduces a modified Jacobi-Davidson method for large sparse eigenvalue problems, incorporating a new correction equation inspired by least squares, with analysis and numerical validation of its convergence and efficiency.

Contribution

It proposes a novel modification to the Jacobi-Davidson method using a least squares motivated correction equation, enhancing convergence analysis for symmetric matrices.

Findings

Convergence properties are established for symmetric matrices.

Numerical experiments demonstrate the method's computational viability.

The modified method offers potential improvements over existing approaches.

Abstract

Each iteration in Jacobi-Davidson method for solving large sparse eigenvalue problems involves two phases, called subspace expansion and eigen pair extraction. The subspace expansion phase involves solving a correction equation. We propose a modification to this by introducing a related correction equation, motivated by the least squares. We call the proposed method as the Modified Jacobi-Davidson Method. When the subspace expansion is ignored as in the Simplified Jacobi- Davidson Method, the modified method is called as Modified Simplified Jacobi-Davidson Method. We analyze the convergence properties of the proposed method for Symmetric matrices. Numerical experiments have been carried out to check whether the method is computationally viable or not.

Equations

$$(I-u{u}^{\ast })(A-\theta I)(I-u{u}^{\ast })t=-r,\text{leavevmode}\text{leavevmode}t\perp u$$

$$v_{k+1}=\frac{(I-V_k{V}_k^{\ast })t}{\Vert (I-V_k{V}_k^{\ast })t\Vert }$$

$$\frac{u+(I-u{u}^{\ast })t}{\Vert u+(I-u{u}^{\ast })t\Vert }$$

$$A=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 2& 0\\ 0& 2& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

$$B:=(I-u{u}^{\ast })(A-\theta I)(I-u{u}^{\ast })=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

$$\Vert (A-\theta I)u+(A-\theta I)(I-u{u}^{\ast })t{\Vert }^2$$

$$(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)(I-u{u}^{\ast })t=-(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)u$$

$$(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)(I-u{u}^{\ast })t+h^{2}t=-(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)u$$

$$\Vert (A-\theta I)u+(A-\theta I)(I-u{u}^{\ast })t{\Vert }^2+h^2\Vert t{\Vert }^2$$

$$(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)(u+(I-u{u}^{\ast })t)=0$$

$$(I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}(u+(I-uu^{\ast})t)=0$$

$$(I-u{u}^{\ast })E(u+(I-u{u}^{\ast })t=0$$

$$\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}(u+(I-uu^{\ast})t)=k_{1}u$$

$$\Vert (A-\theta I)(u+(I-u{u}^{\ast })t){\Vert }^2+(u+(I-u{u}^{\ast })t{)}^{\ast }E(u+(I-u{u}^{\ast })t)$$

$$=k_1(u+(I-u{u}^{\ast })t{)}^{\ast }u$$

$$\Vert (A-\theta I)(u+(I-u{u}^{\ast })t){\Vert }^2=-(u+(I-u{u}^{\ast })t{)}^{\ast }E(u+(I-u{u}^{\ast })t)+k_1$$

$$-2(I-uu^{\ast})E(u+(I-uu^{\ast})t+(I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}u=0$$

$$(I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}(I-uu^{\ast})t=0$$

$$\Vert (A-\theta I)(I-u{u}^{\ast })t\Vert =0$$

$$-(I-uu^{\ast})(E+E^{\ast})\big{(}u+(I-uu^{\ast})t\big{)}+(I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}u=0$$

$$-(I-uu^{\ast})E^{\ast}\big{(}u+(I-uu^{\ast})t\big{)}+(I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}u=0$$

$$-k(I-u{u}^{\ast })w+(I-u{u}^{\ast })(A-\theta I{)}^{\ast }(A-\theta I)u=0$$

$$-k(I-u{u}^{\ast })w+(A-\theta I{)}^{\ast }(A-\theta I)u=\Vert (A-\theta I)u{\Vert }^{2}u$$

$$\big{(}\leavevmode\nobreak\ (A-\theta I)^{\ast}(A-\theta I)\big{)}(u+(I-uu^{\ast})t)=k_{1}u$$

$$\frac{(A-\theta I{)}^{-2}u}{\Vert (A-\theta I{)}^{-2}u\Vert }$$

$$\Vert (A-{\rho }_{k}I)(u_k+(I-u_k{u}_k^{\ast })t_k){\Vert }^2\le \Vert (A-{\rho }_{k}I)u_k{\Vert }^2$$

$$u_{k+1}=\frac{u_k+(I-u_k{u}_k^{\ast })t_k}{(1+\Vert (I-u_k{u}_k^{\ast })t_k{\Vert }^2{)}^{\frac{1}{2}}}$$

$$\Vert (A-{\rho }_{k}I)u_{k+1}{\Vert }^2=\frac{\Vert (A-{\rho }_{k}I)(u_k+(I-u_k{u}_k^{\ast })t_k){\Vert }^2}{1+\Vert (I-u_k{u}_k^{\ast })t_k{\Vert }^2}$$

$$\le \Vert (A-{\rho }_{k}I)(u_k+(I-u_k{u}_k^{\ast })t_k){\Vert }^2$$

$$\displaystyle\|(A-\rho_{k+1}I)u_{k+1}\|^{2}=\|\big{(}(A-\rho_{k}I)-(\rho_{k+1}-\rho_{k})I\big{)}u_{k+1}\|^{2}$$

$$\displaystyle=\|(A-\rho_{k}I)u_{k+1}\|^{2}-|\rho_{k+1}-\rho_{k}|^{2}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\Vert (A-{\rho }_{k+1}I)u_{k+1}{\Vert }^2\le \Vert (A-{\rho }_{k}I)u_{k+1}{\Vert }^2\le \Vert (A-{\rho }_{k}I)u_k{\Vert }^2$$

Full text

A modification of the Jacobi-Davidson method

††thanks: Research supported by CSIR India(09/084(0563)/2010 EMR-I) and NBHM India(2/40(3)/2016/R&D-II/9602)

Mashetti Ravibabu

Department of Computational and Data Sciences

Indian Institute of Science

Bengaluru

India-560012 mashettiravibabu2@gmail.com

Each iteration in Jacobi-Davidson method for solving large sparse eigenvalue problems involves two phases, called subspace expansion and eigen pair extraction. The subspace expansion phase involves solving a correction equation. We propose a modification to this by introducing a related correction equation, motivated by the least squares. We call the proposed method as the Modified Jacobi-Davidson Method. When the subspace expansion is ignored as in the Simplified Jacobi-Davidson Method, the modified method is called as Modified Simplified Jacobi-Davidson Method. We analyze the convergence properties of the proposed method for Symmetric matrices. Numerical experiments have been carried out to check whether the method is computationally viable or not.

**keywords:**Jacobi-Davidson method, Subspace expansion, Eigen pair extraction

1 Introduction

Projection methods are widely used for solving large sparse eigenvalue problems. Jacobi-Davidson method is a quite well known projection method that approximates the smallest eigenvalue or eigenvalue near a given shift of a symmetric matrix. This method starts with an arbitrary initial nonzero vector and involves two phases. In the first phase, the eigenvector and eigenvalue approximations are obtained from the existing subspace by applying a projection to the given matrix. From these approximate eigenvalues, we select the one with desired properties and its corresponding eigenvector approximation. The second phase involves solving the correction equation to obtain a vector from the orthogonal complement of the selected eigenvector approximation. An existing subspace is expanded by adding this vector after orthogonalizing it against the existing subspace. In the proposed method, we modify the correction equation using least squares. A vector is determined from the orthogonal complement of a selected eigenvector approximation so that the norm of the residual associated with the resultant of the selected eigenvector approximation and the vector determined from its orthogonal complement is minimum.

In Section 2, we briefly review Jacobi-Davidson method. In Section 3, a new correction equation is introduced based on least squares heuristics, and convergence properties of the proposed method are discussed. In Section 4, we report the results of numerical experiments which were carried out to understand the viability of the method developed in Section 3. In Section 5, we consider the method used with restarting, and report the results of numerical examples. Section 6 concludes the paper.

2 Jacobi-Davidson and simplified Jacobi-Davidson method

We briefly describe Jacobi-Davidson method. Let $A$ be a given matrix. Assume that we have already computed a matrix $V_{k}$ with $k$ orthonormal columns, which span the existing subspace of dimension $k$, starting with an initial nonzero vector $v_{1}$. Approximations to eigenvalues and eigenvectors of $A$ are obtained from those of the matrix $H_{k}=V_{k}^{\ast}AV_{k}$ as follows. If $(\theta,y)$ is an eigen pair of $H_{k}$, then $\theta$ is called a Ritz value and $u=V_{k}y$ is called a Ritz vector, and the associated residual is $r=Au-\theta u$. Since $V_{k}$ has orthonormal columns, $V_{k}^{\ast}r=V_{k}^{\ast}AV_{k}y-\theta V_{k}y=0$. Thus, the residual is orthogonal to the existing subspace. We then select one of the Ritz values with desired properties and its corresponding Ritz vector. Using this vector, we solve the following correction equation, for the vector $t$:

[TABLE]

The vector $t$ is orthogonalized with respect to the existing subspace of dimension $k$ spanned by $k$ orthonormal columns of the matrix $V_{k}$ and then normalized to get the vector

[TABLE]

Then $v_{k+1}$ is added to the existing subspace of dimension $k$ and $V_{k}$ is updated to $V_{k+1}=[V_{k}\leavevmode\nobreak\ \leavevmode\nobreak\ ,v_{k+1}]$. An eigen pair of $H_{k+1}$ is then taken as an improved eigen pair. Since $A$ is very large and sparse, Equation (1) is solved using a Krylov subspace method for linear systems, such as GMRES or FOM. The algorithm is continued until the norm of the residual satisfies some fixed tolerance, that is, $\|r|\|_{2}\leq tol$.

In this process, if we ignore the subspace expansion, then the method is called the Simplified Jacobi-Davidson method. Here, for an eigenvector approximation $u$ and the solution $t$ of the correction equation (1), the vector $\frac{u+t}{\|u+t\|}$ is considered as the new eigenvector approximation, The process is continued by replacing $u$ and $\theta$ in the correction equation (1), with $\frac{u+t}{\|u+t\|}$ and the Rayleigh quotient of the vector $\frac{u+t}{\|u+t\|},$ respectively. The process terminates when norm of the residual vector satisfies a pre-set tolerance. It differs from Jacobi-Davidson method by ignoring the information from previous vectors, where as in Jacobi-Davidson method, a subspace is formed from these vectors.

3 Modification to the Simplified Jacobi-Davidson method

As we have seen in Section 2, Simplified Jacobi Davidson method uses at each step, the solution of the correction equation (1), where $u$ is a known eigenvector approximation, $\theta$ is its Rayleigh quotient with respect to the matrix $A,$ and the vector

[TABLE]

is the new eigenvector approximation. It is a general belief that the matrix $(I-uu^{\ast})(A-\theta I)(I-uu^{\ast})$ is well-conditioned compared to the matrix $(A-\theta I);$ but this is not correct as shown in the following Example.

Example 1**.**

Consider the matrix

[TABLE]

Take the vector $u$ as $(0\leavevmode\nobreak\ 0\leavevmode\nobreak\ 1\leavevmode\nobreak\ 0)^{t}$. Then $\theta=u^{\ast}Au=0,$ $r=(A-\theta I)u=(0\leavevmode\nobreak\ 2\leavevmode\nobreak\ 0\leavevmode\nobreak\ 0)^{t},$ and

[TABLE]

Here, $A-\theta I$ is non-singular but the matrix $B$ is singular. Further, $-r$ is not in the range space of $B$. So there exists no vector $t$ satisfying the correction equation (1).

Example 1 illustrates the situation where the correction equation has no solution. To avoid this, we propose the following method of choosing a new eigenvector approximation. We call this method as the Modified Simplified Jacobi Davidson method (MSJD). It has the following two steps.

Algorithm 1**.**

MSJDStep 1: For a given eigenvector approximation $u$ of a matrix $A$, find a vector $t$ that minimizes

[TABLE]

*where $\theta$ is the Rayleigh quotient of $u$ with respect to the matrix $A$.

Step 2: Take the new eigenvector approximation as   $\dfrac{u+(I-uu^{\ast})t}{\|u+(I-uu^{\ast})t\|}$.*

In Step 1, the vector $t$ is determined by solving the following normal equation:

[TABLE]

Notice that, if $(A-\theta I)$ is symmetric, then (2) can be obtained from (1) by replacing the matrix $(A-\theta I)$ with $(A-\theta I)^{2},$ and then applying the orthogonal projection $(I-uu^{\ast})$ to both the sides.

The linear system  (2) may be solved by well known methods such as Gaussian elimination, Conjugate gradients, and etc. When $\theta$ is very close to an eigenvalue $\lambda$ of $A$, the matrix $(I-uu^{\ast})(A-\theta I)^{2}(I-uu^{\ast})$ is expected to be more ill-conditioned than $(I-uu^{\ast})(A-\theta I)(I-uu^{\ast})$. Such ill-conditioned problems can be solved by using well known regularization techniques, for example, Tikhonov Regularization. Tikhonov regularization method, instead of solving the normal equations, requires the solution of the following equation:

[TABLE]

where $h^{2}$ is a parameter chosen in such a way that as $h$ tends to zero, the solution of the above equation converges to the solution of the least squares problem in Step 1. This is equivalent to finding the vector that minimizes the functional

[TABLE]

Here we propose a new way, which avoids problems associated with parameter selection. For this purpose, we rewrite Equation (2) as

[TABLE]

When $\theta$ is close to an eigenvalue of $A$, the system matrix in Equation (3) is ill-conditioned. To avoid ill-conditioning, perturb the matrix $(A-\theta I)^{\ast}(A-\theta I)$ by a matrix $E$, and write the perturbed equation as

[TABLE]

We require that $\|(A-\theta I)u+(A-\theta I)(I-uu^{\ast})t\|$ be minimum. Using Equation (2), we simplify (4) to obtain

[TABLE]

Matrices that transform the vector $u+(I-uu^{\ast})t$ to a vector parallel to $u$ are possible candidates for $E$ in Equation (5). Also, $E\leavevmode\nobreak\ :=\leavevmode\nobreak\ kuw^{\ast}$ for any scalar $k$ and any vector $w,$ satisfies (5). For the choice $w=u$, the matrix $E$ becomes symmetric and it may have computational advantages. However this simple choice is not possible, in general.

Observation 1**.**

The choice $w=u$ is not possible unless $(\theta,(I-uu^{\ast})t)$ is an eigen pair of $A$ or $u$ is a right singular vector of $(A-\theta I)$.

Proof.

Rewrite Equation (4) as

[TABLE]

where k_{1}:=u^{\ast}\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}(u+(I-uu^{\ast})t) is a scalar. Taking the inner product with $u+(I-uu^{\ast})t,$ we obtain

[TABLE]

Since $\|u\|=1$, this equation is simplified to

[TABLE]

In order that the left hand side of Equation (6) is minimum over all vectors $t$, the gradient of the right hand side of Equation (6) with respect to the vector $t$ vanishes. Therefore,

[TABLE]

Using Equation (5), we have (I-uu^{\ast})\big{(}(A-\theta I)^{\ast}(A-\theta I)+E\big{)}u=0. This together with Equation (4) gives

[TABLE]

If $E$ is of the form $kuu^{\ast}$, then $(I-uu^{\ast})(A-\theta I)^{\ast}(A-\theta I)(I-uu^{\ast})t=0$. This implies

[TABLE]

Therefore, the choice $w=u$ is not possible, unless $(\theta,(I-uu^{\ast})t)$ is an eigen pair of $A$ or $(I-uu^{\ast})t=0$. From Equation (3), it is clear that if $(I-uu^{\ast})t=0,$ then $u$ is a right singular vector of $(A-\theta I)$. ∎

Observation 2**.**

If $E=kuw^{\ast}$ satisfies Equation (5), then the component of $w$ orthogonal to $u$ is parallel to the vector \big{(}(A-\theta I)^{\ast}(A-\theta I)-\|(A-\theta I)u\|^{2}I\big{)}u.

Proof.

Vanishing of the gradient of right hand side of Equation (6) gives

[TABLE]

Using Equation (5), we have

[TABLE]

With $E=kuw^{\ast}$ and $\|u\|=1$, we have

[TABLE]

It gives

[TABLE]

Therefore, k(I-uu^{\ast})w=\big{(}(A-\theta I)^{\ast}(A-\theta I)-\|(A-\theta I)u\|^{2}I\big{)}u. ∎

Observation-2 gives the direction of the component of $w$ orthogonal to the vector $u.$ Moreover, it is clear that the component of $w$ in the direction of $u$ can be chosen as $\alpha u$, for some scalar $\alpha$. Therefore, the vector $w$ in Observation-2 is of the form w=\alpha u+\big{(}(A-\theta I)^{\ast}(A-\theta I)-\|(A-\theta I)u\|^{2}I\big{)}u. A scalar $\alpha$ can be chosen by imposing an extra condition on the vector $w$ such as $\|w\|_{2}=1.$ Further, notice that when $E=kuw^{\ast},$ solutions of Equations (3) and (4) coincide. Then Equation (3) can be rewritten as

[TABLE]

with $k_{1}=u^{\ast}(A-\theta I)^{\ast}(A-\theta I)(u+(I-uu^{\ast})t)$. If $A$ is symmetric, then from Equation (7) we have $u+(I-uu^{\ast})t=k_{1}(A-\theta I)^{-2}u$. In this case, the new eigenvector approximation is

[TABLE]

Next, we deal with the convergence of the norms of residual vectors in MSJD method.

Theorem 3.1**.**

Let $u_{k}$ denote an approximation to an eigenvector $x$ corresponding to the eigenvalue $\lambda$ of $A$ at the $k$th iteration of the MSJD method. Let $\rho_{k}:=\rho(u_{k})$ be the Rayleigh quotient of $u_{k}$ with respect to the matrix $A$. Then the sequence of residual norms $\{\|(A-\rho_{k}I)u_{k}\|^{2}\}$ and the sequence $\{\|(A-\rho_{k}I)u_{k+1}\|^{2}\}$ converge to the same limit. Further the sequence of absolute differences $|\rho_{k+1}-\rho_{k}|$ between Rayleigh quotients in two consecutive iterations of MSJD method converges to [math] as $k\rightarrow\infty.$

Proof.

From Step 1 of MSJD method, it follows that

[TABLE]

And from Step 2, we have

[TABLE]

Notice that $\|u_{k}\|=1$ and $\|u_{k}+(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}=1+\|(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}\geq 1.$ Then Equations (9) and (10) give

[TABLE]

As $\|u_{k+1}\|=1$ and $\rho_{k+1}$ is a Rayleigh quotient of the vector $u_{k+1}$ with respect to $A$, we have $\langle(A-\rho_{k}I)u_{k+1},(\rho_{k+1}-\rho_{k})u_{k+1}\rangle=(\rho_{k+1}-\rho_{k})^{2}.$ Thus,

[TABLE]

Therefore $\|(A-\rho_{k+1}I)u_{k+1}\|^{2}\leq\|(A-\rho_{k}I)u_{k+1}\|^{2}$. Combining this inequality with Equations (3) and (3), we obtain

[TABLE]

It shows that $\{\|(A-\rho_{k}I)u_{k}\|^{2}\}$ is a monotonically decreasing sequence of non-negative real numbers. Suppose it converges to $\alpha$. By Sandwich theorem, the sequence $\{\|(A-\rho_{k}I)u_{k+1}\|^{2}\}$ also converges to $\alpha$. Using Equation (3), we conclude that $|\rho_{k+1}-\rho_{k}|\rightarrow 0$ as $k\rightarrow\infty.$ ∎

Though $|\rho_{k+1}-\rho_{k}|\rightarrow 0$ is enough for computational purposes, it does not imply that $\{\rho_{k}\}$ is a convergent sequence. In order to prove a result on convergence alternatives, we use the following two lemmas.

Lemma 3.2**.**

Let $u$ be unit vector. Let $\alpha$ be the Rayleigh quotient of $u$ with respect to matrix $A$. Let $s:=u+\tau(I-uu^{\ast})t$, where $\tau$ is chosen such that Rayleigh quotient of $s$ is minimum over the space of vectors spanned by $u$ and $(I-uu^{\ast})t.$ Write $J_{u,s}:=(I-uu^{\ast})(A-\rho(s)I)(I-uu^{\ast}.)$ Then the following relationship holds:

[TABLE]

For the proof of Lemma 3.2, see Theorem 2.6 in [3].

Lemma 3.3**.**

Let $\rho_{k}$, $u_{k}$ and $t_{k}$ be as in Theorem (3.4). Let $\tau$ be such that for the vector $s:=u_{k}+\tau(I-u_{k}u_{k}^{\ast})t_{k}$ norm of the residual $\frac{\|(A-\rho_{k}I)s\|}{\|s\|}$ is minimum over the subspace spanned by $u_{k}$ and $(I-u_{k}u_{k}^{\ast})t_{k}$. Then

[TABLE]

and

[TABLE]

Proof.

For the vectors $u_{k}$ and $s,$ we have

[TABLE]

and $u_{k}-s=-\tau(I-u_{k}u_{k}^{\ast})t_{k}$. Using Equation (13) for the matrix $(A-\rho_{k}I)^{\ast}(A-\rho_{k}I)$ and the vectors $u_{k}$ and $s,$ we obtain

[TABLE]

It follows that

[TABLE]

Since $\frac{\|(A-\rho_{k}I)s\|^{2}}{\|s\|^{2}}\leq\|(A-\rho_{k}I)u_{k}\|^{2}$, Equation (16) yields the inequality in Equation (14). From Equation (9), we have

[TABLE]

It gives $\|(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}\leq 1$. ∎

Theorem 3.4**.**

Let the scalar $\rho_{k}$ and the vector $u_{k}$ be as in Theorem 3.1. Then the sequence of eigenvector approximations in MSJD method converges to either an eigenvector of the matrix $A$ corresponding to the eigenvalue $\rho,$ or to a vector in an invariant subspace corresponding to eigenvalue $\rho,$ or to a vector in an invariant subspace spanned by eigenvectors corresponding to eigenvalues $\rho\pm\alpha.$

Proof.

We see that

[TABLE]

Taking inner product with $t_{k}$ on both sides of equation (2), we get

[TABLE]

Therefore,

[TABLE]

Using Equation (10), we obtain

[TABLE]

It implies that

[TABLE]

Due to Theorem 3.1, the left hand side of Equation (3) converges to zero as $k\rightarrow\infty$. Therefore

[TABLE]

From Equation (14), we have

[TABLE]

Therefore, $\|(A-\rho_{k}I)u_{k}\|^{2}\|(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}\rightarrow 0.$ Since $\|(A-\rho_{k}I)u_{k}\|^{2}$ is a monotonically decreasing sequence of nonnegative terms, it converges to a non-negative real number. Suppose it converges to $\alpha^{2}\geq 0.$ We have two cases:

[TABLE]

Note that in both the cases, Equation (8) can be written as

[TABLE]

Case 1: $\|(I-u_{k}u_{k}^{\ast})t_{k}\|\rightarrow 0$

From Step 2 in Algorithm 1, we have

[TABLE]

Since $u_{k}$ is orthogonal to $(I-u_{k}u_{k}^{\ast})t_{k}$ and $\|u_{k}\|=1$, we have

[TABLE]

Writing $\alpha_{k}=\sqrt{1+\|(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}},$ we see that $u_{k+1}-\frac{u_{k}}{\alpha_{k}}=\frac{(I-u_{k}u_{k}^{\ast})t_{k}}{\alpha_{k}}.$ Since $\|(I-u_{k}u_{k}^{\ast})t_{k}\|\rightarrow 0,\leavevmode\nobreak\ \alpha_{k}\rightarrow 1$ as $k\to\infty.$ And then $u_{k+1}-u_{k}\rightarrow 0.$ From Equation (18), we thus obtain

[TABLE]

As $u_{k+1}-u_{k}\rightarrow 0,$ $(A-\rho_{k}I)^{2}u_{k}=\alpha^{2}u_{k}\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{for large}\leavevmode\nobreak\ \leavevmode\nobreak\ k.$ Then

[TABLE]

Therefore, In this case also, $\rho_{k}$ is equal to either $\lambda_{i}+\alpha$ or $\lambda_{i}-\alpha$ for one or more eigenvalues of $A.$

Case 2: $\|(A-\rho_{k}I)u_{k}\|^{2}\rightarrow 0$ and the sequence $\{u_{k}\}$ converges to $u.$ Then $u$ is an eigenvector corresponding to $\rho.$

Case 3: $\|(A-\rho_{k}I)u_{k}\|^{2}\rightarrow 0$ and the sequence $\{u_{k}\}$ is not convergent.

Then Equation (10) gives

[TABLE]

From Equation (15), we have $\|(I-u_{k}u_{k}^{\ast})t_{k}\|^{2}\leq 1$. This together with Equation  (18) gives

[TABLE]

which implies that $\|(A-\rho_{k}I)^{2}u_{k+1}\|\rightarrow 0$. Then, $\|(A-\rho_{k}I)u_{k}\|^{2}\to 0$ and $\|(A-\rho_{k}I)^{2}u_{k+1}\|\to 0$ as $k\rightarrow\infty.$ It follows that, as $k\rightarrow\infty,$ $u_{k}$ is in the invariant subspace corresponding to the eigenvalue $\rho$ . ∎

We show that if $\{u_{k}\}$ converges to an eigenvector corresponding to the eigenvalue $\lambda$, then the order of convergence is $5$. The proof follows the line of proof of Theorem 4.7.1 in [5].

Theorem 3.5**.**

Let the sequence $\{u_{k}\}$ generated by MSJD method converge to an eigenvector $x$ of $A$ corresponding to the eigenvalue $\lambda$. Let $\Phi_{k}$ denote the angle between vectors $u_{k}$ and $x$. Then $|\Phi_{k+1}|\leq|\Phi_{k}|^{5}$ for large $k.$

Proof.

Since $\Phi_{k}$ is the angle between the unit vectors $u_{k}$ and $x$, $u_{k}$ can be written as

[TABLE]

Where $v_{k}^{\ast}x=0$ and $\|v_{k}\|=\|x\|=1.$ Let $\rho_{k}:=\rho(u_{k})$ be the Rayleigh quotient of $A$ with respect to $u_{k}$. Pre-multiply Equation (19) with $(A-\rho_{k}I)^{-2}$ and use Spectral mapping theorem to obtain

[TABLE]

From Equation (8), the next eigenvector approximation $u_{k+1}$ in MSJD method satisfies the following:

[TABLE]

Using Equation (20), we have,

[TABLE]

Let $v_{k+1}$ be a unit vector such that $u_{k+1}=x\cos\Phi_{k+1}+v_{k+1}\sin\Phi_{k+1}$ and $v_{k+1}^{\ast}x=0.$ Then

[TABLE]

[TABLE]

[TABLE]

We have

[TABLE]

From equation (19), we have

[TABLE]

Since $Ax=\lambda x$ and $v_{k}^{\ast}x=0,$ it follows that $\rho_{k}=\lambda\cos^{2}\Phi_{k}+\rho(v_{k})\sin^{2}\Phi_{k}$. Therefore

[TABLE]

Using this in Equation (21), we obtain

[TABLE]

Since $v_{k}$ is orthogonal to $x,$

[TABLE]

Since $\Phi_{k}\rightarrow 0$, Equation (22) implies that $|\lambda-\rho_{k}|\rightarrow 0$ as $k\rightarrow\infty$. As $\rho_{k}$ converges to $\lambda$ as $k\rightarrow\infty$, there exists a real number $M$ such that

[TABLE]

for all $\lambda_{i}\neq\lambda$. Now Equations (23) and (24) together with $\tan\Phi=O(\Phi)$ and $\sin(\Phi)=O(\Phi)$ for small $\Phi$, imply that $|\Phi_{k+1}|\leq|\Phi_{k}|^{5}$ for large $k.$ ∎

4 Numerical experiments

The MSJD method has been tested on many numerical examples for checking whether it really works, using Matlab R2014A on an Intel core 3 processor. Out of these we report three examples for demonstrating various features. We compare the performance of the Jacobi-Davidson method and the proposed MSJD method by solving the corresponding correction equations

[TABLE]

[TABLE]

using Gaussian elimination as well as with the approximate solution obtained after a few steps of GMRES. Since the matrices in the left hand side of the correction equation (25) in Jacobi-Davidson method may become ill-conditioned, especially, when the matrix $A$ has multiple eigenvalues, it has been proposed in [6] to use the following equivalent form of (25):

[TABLE]

In a similar vein, we define the following correction equation in MSJD method:

[TABLE]

which is theoretically equivalent to (26).

In the following examples, we compare the performance of the correction equation (25) in Jacobi-Davidson method with the correction equation (26) in Modified Jacobi-Davidson method. We also compare the numerical results obtained using the correction equation (27) in Jacobi-Davidson method with the correction equation (28) in Modified Jacobi-Davidson method. We also demonstrate in the following examples that the performance of JD method is numerically different for the two theoretically equivalent correction equations (25) and (27). This is because of the presence of rounding errors in floating point arithmetic. In a similar vein, we also demonstrate that the MJD performs numerically different for the two theoretically equivalent correction equations (26) and (28).

In all the figures,“JD" means either the correction equation (25) or (27) in Jacobi-Davidson method is used whereas “MJD" means either the correction equation (26) or (28) is used. In order to check the performance of the proposed method, without restarting, we consider the following algorithm:

Algorithm 2**.**

*Unrestarted algorithm1. Solve the correction equation, either (25) or (27) in Jacobi-Davidson method or either (26) or (28) in Modified Jacobi Davidson method.

2. Expand the subspace with the vector obtained in Step 1, as explained in Section 2.

3. Apply Rayleigh-Ritz projection/Harmonic Rayleigh-Ritz projection for eigenvalues, and calculate the norm of residuals associated with refined Ritz/Harmonic Ritz vectors. If they reach $`tol^{\prime},$ stop. Otherwise, go to Step 1.*

When a good approximation to eigenvalues in the interior of the spectrum is required we use Harmonic Rayleigh-Ritz projection [2] in Algorithm 1. The results on performance of Algorithm 1 for three Examples are reported below.

Example 2**.**

Consider the diagonal matrix $A$ of order $100$ with diagonal elements as $(\frac{j}{100})^{2}-0.8$ for $j=1,2...100$. Take the initial vector as the vector with each entry equal to $1$. It is required to compute an eigenvalue with the smallest absolute value.

This matrix is same as in Example 3 in [6] and is used to show that the proposed method can also be used in computing interior eigenvalues. As the method is not restarted, in this and next example, iteration number denotes the size of subspace used in Harmonic Projection. For this example, we used Harmonic projection in Algorithm 2 and Harmonic Ritz vectors in the correction equations (25) and (26).

With an approximate solution of Equations (25) and (26) obtained using $8$ steps of GMRES, the convergence to the required eigenvalue is observed at iterations $69$ and $66,$ respectively. The $\log_{10}$ of the norm of the residuals are $-12.26$ and $-13.65,$ respectively. The numerical results with the correction equation (25) are shown in Figures 2-2. Figures 4-4 show the results, when the correction equation (26) is used. In Figures 2 and 4, $`\ast^{\prime}$ represents the exact eigenvalues of $A$. From these Figures, it is clear that convergence results of approximate eigenvalues in both the methods are on par.

Similarly, using the solutions of the correction equations (27) and (28) obtained using $8$ steps of GMRES, convergence occurs at iterations $57$ and $75,$ with $\log_{10}$ of residual norms as $-12.15$ and $-13.78,$ respectively. Convergence for these cases are shown in Figures 6-6 and 8-8, respectively. Using the correction equation (27), convergence of norms of residual vectors is faster when compared to MJD using the correction equation (28).

If Gaussian elimination is used instead of GMRES to solve the correction equations (27) and (28), with the correction equation (27), convergence to the required eigenvalue occurred at $57^{th}$ iteration, whereas with the correction equation (28), convergence is achieved at $46^{th}$ iteration. Thus, when the correction equations (27) and (28) are solved by using the Gaussian elimination, we noticed that with the correction equation (28) the convergence to the required eigenvalue is faster than that with the correction equation (27). This is exactly opposite to the scenario that we observed when the correction equations (27) and (28) are solved approximately by using the GMRES.

In Examples 2, even though the matrix is symmetric, we used the GMRES to solve the correction equations approximately as we did not take the advantage of this for generating matrices $H_{k}$ in JD and MJD methods (For details, see Section-2). In the following example, we consider a non-Hermitian matrix. We observe that like Jacobi-Davidson method, the modified method is also useful to approximate non-real eigenvalues.

Example 3**.**

The matrix in this Example is a block diagonal matrix $diag\{A1,A2\}$ where

[TABLE]

and $A2$ is the matrix of Example 2. Again, each entry of the initial vector is taken as $1$. We apply Harmonic projection with shift $0.81+0.08i$ to find the eigenvalue $0.8+0.1i$.

Here, we compare the results obtained using the correction equations (27) and (28). In both the approaches the resulting linear systems are solved using Gaussian elimination. With these two correction equations the convergence to the desired eigenvalue occurred at $44^{th}$ and $48^{th}$ iterations, respectively. In both the cases, apart from the desired eigenvalue,the algorithm also finds the eigenvalue $0.8-0.1i$ which is far from the shift. The same is observed, using Gaussian elimination for solving the correction equation (25). For the correction equation (27), Figures 10 and 10 show the convergence history of the real parts and imaginary parts of the harmonic Ritz values, respectively. Similarly, for the correction equation (28), the Figures 12 and 12 show the convergence history of the harmonic Ritz values. In all these figures, we have used the same symbols as in the previous example.

For the correction equations (27) and (28), the convergence of norm of residual vectors are shown in the Figures 14 and 16, respectively.For the correction equations (27) and (28), the harmonic Ritz values at the final iteration where the convergence occurred are shown in Figures 14 and 16, respectively. In these Figures, we marked the exact eigenvalues of $A$ with squares. Dots represent the Harmonic Ritz values obtained at the iteration where convergence occurs. It is easy to observe from these two figures that the accurate approximation to larger number of eigenvalues are obtained using the correction equation (28)when compared to the correction equation (27).

In case of solving the correction equation (26) using the Gaussian elimination, the eigenvalue approximations converges to an eigenvalue other than the desired one. When the correction equations (25), (27) and (26), (28) are solved approximately using $10$ steps of GMRES, the approximation converges to the eigenvalue $2.000+0.000i$, which is not the desired eigenvalue.

The next example shows that the new method also works for large sparse matrices.

Example 4**.**

We consider the matrix $A$ as SHERMAN4, a sparse matrix of order $1104$, taken from Harwell-Boeing set of test matrices. The smallest eigenvalue $0.030726$ (accurate upto $5$ decimal places) is required. MATLAB command ‘eigs’ produces the result as $3.072570776499973e-02$. All entries in the initial vector are equal to $1$. Rayleigh-Ritz projection and Refined Ritz vectors are used for approximating the eigen pairs.

In Table 1, we give the numerical results of Jacobi-Davidson and MSJD method for the matrix Sherman4. Fast convergence is observed when Gaussian elimination is used for solving the correction equation (27) compared to using the other correction equations. The comparison of convergence for the correction equations (27) and (28) is done in Figures 18-18.

With Gaussian elimination for solving the correction equation (26), the convergence occurs at the iteration number $11,$ whereas with the correction equation (25), the convergence is reached at the iteration number $10.$ In both the approaches, the obtained eigenvalue approximation is accurate upto $10$ decimal places. The Comparison of convergence of eigenvalue approximations and norms of residual vectors in these two cases is done in Figures 20 and 20, respectively. We observe that with Gaussian elimination for solving the correction equations in the Jacobi-Davidson and the new method, the results are on par

5 With restarting

It is well known that for symmetric matrices, Rayleigh-Ritz projection over large subspaces may give good approximations to an eigen pair. But as the size of a subspace increases, the cost associated with computing an eigen pair also increases. Further, if the size of the given matrix is very large, the space complexity in the computation may become practically unmanageable. For this reason, the method with restarting is favourable. In the following examples, we check the performance of the MSJD method with restarting.

Example 5**.**

Consider the matrix $A$ as the first Example in [6], which is a diagonally dominant tridiagonal matrix of order $200$ with diagonal elements $a_{i,i}=2.4+i/2$ for $i<200$, $a_{200,200}=2.4+200/1.5,$ and with each entry on the super-diagonal and sub-diagonal as $1$. We take the initial vector $v1=(0.03,0.03,....0.03,1)^{\ast}$ as in [6]. The corresponding Rayleigh quotient with respect to $A$ is $1.632770531196111e+02$. Our goal is to approximate the largest eigenvalue. Such a preliminary approximation is obtained by using Matlab command ‘eig’, which computes it as $2.561474561181774e+02$.

We compare the performance of the proposed method with the Jacobi-Davidson method by restarting the algorithm after the size of subspace becomes $3$. Table 2 shows a summary of numerical results.

Using $5$ steps of GMRES to approximate the solution of the correction equation (26) in MJD method gives an eigenvalue approximation near the desired eigenvalue. The Ritz values converge to a spurious eigenvalue from $5^{th}$ restart onwards, but norm of the residual vectors reach the tolerance. The same behaviour is observed with its theoretically equivalent correction equation (28). With correction equations (25) and (27) in Jacobi-Davidson method, from $1^{st}$ restart onwards, the Ritz values stagnated near the desired eigenvalue. Convergence of Ritz values using correction equations (25) and (26) are shown in Figure 22 . For correction equations (27) and (28), they are shown in Figure 24. The dependence of $\log_{10}$ of norms of corresponding residual vectors versus restart numbers are shown in Figures 22 and 24.

We also checked the performance of the proposed method by restarting the algorithm when the size of the subspace for extracting an eigen pair reached $4.$ The approximate solution of correction equations is obtained after $5$ steps of GMRES. Using the correction equation (25), the eigenvalue approximation $2.561474561181778e+02$ is obtained at first restart, that is, without restart. After $4^{th}$ restart, it starts giving spurious eigenvalues whereas with the correction equation (26) the approximate eigenvalue obtained at $2^{nd}$ restart is found to be $2.561474561181784e+02$. Comparison of convergence of approximate eigenvalues in these cases is done in Figure 26.Figure 26 reports the $\log_{10}$ of residual norms.

In a similar vein, using $5$ steps of GMRES for solving the correction equation (27) and with subspace size $4,$ the approximation to the desired eigenvalue is obtained without restart, that is, when size of the subspace reaches $4.$ From $5^{th}$ restart, it starts giving a spurious eigenvalue. Whereas with the correction equation (28), norm of residuals reached the tolerance at $2^{nd}$ restart and the approximation to eigenvalue is $2.561474561181783e+02.$ Results of comparison of convergence of Ritz values and norms of residuals with correction equations (27) and (28) are shown in Figure 28 and Figure 28, respectively.

Example 6**.**

Let the matrix $A=QTQ,$ where $T$ is a tridiagonal matrix of order $100$ with super-diagonal and sub-diagonal elements as $-1,$ and diagonal elements $1$, and $Q$ is the Householder transformation of order $100$ with Householder vector $h$ having entries $h_{i}=\sqrt{i+.45}$ for $i=1,2,....100$. The matrix $A$ is not diagonally dominant. We take each entry in the Initial vector as $1$. The largest eigenvalue is required. Matlab command ‘eig’ produces an approximation to the largest eigenvalue as $3.999032564583972e+00$. See Example 2 in [6].

We first computed eigenvalue approximations by applying Rayleigh-Ritz projection over subspaces of dimension ranging from $3$ to $20$. Using Gaussian elimination for solving the correction equation (27), a good approximation to the desired eigenvalue is obtained at $8^{th}$ restart for the subspace of dimension $15,$ where the Ritz value and residual norm are $3.999032566618813e+00$ and $2.034836893782876e-09,$ respectively. The same accuracy to the desired eigenvalue is also obtained for the subspace of dimension $15$, when the correction equation (28) is solved using Gaussian elimination. Figures  30 and 30, respectively, show the comparison of Ritz values and residual norms in these two cases.

For subspaces of dimensions $11$ and $14$ in Rayleigh-Ritz projection, an accuracy upto machine precision is obtained for the desired eigenvalue at $13^{th}$ and $10^{th}$ iterations, respectively. The results are shown in Figure 32 and the convergence of residual norms associated with this is shown in Figure 32. But, eigenvalue approximations obtained using the correction equation (27) are not accurate upto machine precision for a subspace of any size.

As mentioned earlier, the stagnation phenomenon occurs with the correction equation (28) when the subspace is of size $5,$ from $45^{th}$ iteration onwards with Ritz value $\theta=3.999039854574475e+00.$ This is accurate upto $5$ decimal places with the residual norm as $7.289990498570303e-06$. In this case

[TABLE]

where $x$ is a refined Ritz vector. A right singular vector of $A-\theta I$ corresponding to the singular value $(7.289990498570303e-06)^{2}$ is obtained. Figures 34 and 34 show the convergence of Ritz values and residual norms.

Matrices in the correction equation (25) are found to be ill-conditioned during computation, when exact solutions are required except for subspaces of sizes $3,11,16$ and $17.$ In these exceptional cases, the Ritz values are found to be accurate up to $4$ decimal places. Comparison of convergence behaviour of Ritz values and residual norms are shown in Figures 36-36, for subspaces of sizes $3,11,16$ and $17$.

Using Gaussian elimination for solving the correction equation (26), the exact eigenvalue is obtained at $16^{th}$ iteration, when the subspace of size $10$ is used in the Rayleigh-Ritz projection. Matrices in the correction equation (26) become almost singular only for subspaces of size $6$ and $11$. Convergence of Ritz values and residual norms with various subspace sizes, using the correction equation (26) are shown in Figures 38-38.

When an approximate solution of Equations (25), (27) and (26), (28) are obtained using $5$ steps of GMRES method, an eigenvalue approximation obtained is found to be accurate upto $4$ decimal places. It may be explained as follows. A possible explanation is that due to the correction equations in Jacobi Davidson method, the difference between Ritz values in two consecutive iterations is very high whereas with new correction equations, the difference is low leading to slow convergence.

6 Conclusion and future work

In this paper, we have proposed a modification to the subspace expansion phase in Jacobi-Davidson method for computing approximate eigenvalues of a large sparse matrix. The modification uses the heuristic of least squares. Theoretically, the modification has the advantage that it is still applicable to the cases when the correction equation obtained in Jacobi-Davidson method results in a singular system matrix. Further, the modified method is theoretically equivalent to the Alternating Rayleigh quotient iteration, which converges globally and proposed by B.N. Parlett in [4]. To check whether the modification performs well computationally, we have considered many bench mark examples. It is observed that the over all performance of the modified algorithm is well comparable with the Jacobi-Davidson method. Along with the required eigenvalue, approximations to other eigenvalues are also obtained. In case the proposed modified method exhibits slow convergence, compared to Jacobi-Davidson method, it gives good approximation to many eigenvalues including the desired one. The slow convergence is attributed to the clustering of eigenvalues near the current approximation. While Jacobi-Davidson method jumps away from this cluster resulting in an approximation to a different eigenvalue than the desired one, the proposed method approaches slowly towards the desired eigenvalue. It has been observed that when stagnation occurs in the modified method, an approximation to the right singular vector of a matrix $(A-\theta I)$ is obtained. This is observed in Example 6, where the norm of the residual vector is a singular value. When the norm of the residual coincides with a smallest singular value, the obtained vector is likely to be a good approximation to an eigenvector of the matrix, associated with an approximate eigenvalue $\theta.$ The theory about this coincidence is yet to be developed.