A Block-Preconditioned Inexact Linear Solver for Computing the Complex Eigenpairs of a Large Sparse Matrix ()
1. Introduction
Let
be a large sparse, real n by n nonsymmetric matrix and
a symmetric positive definite matrix. In this paper, we consider the problem of computing the eigenpair
from the following generalized complex eigenvalue problem
(1)
where
is the eigenvalue of the pencil
and
its corresponding complex eigenvector. We assume that the eigenpair of interest
is algebraically simple, with
the corresponding left eigenvector such that [1]
(2)
By adding the normalization
(3)
to (1), the combined system of equations can be expressed in the form
as
(4)
Note that
is real since
is symmetric and positive definite. This results in solving a system of nonlinear n complex and one real equation for the
complex unknowns
. The reason we cannot use Newton’s method to solve (4) is because
in the normalization
is not differentiable.
Recall that for a real eigenpair
, (4) is
real equations for
real unknowns and Newton’s method for solving (4) involves the solution of the
square linear systems
(5)
for the
real unknowns
, and updating
. Secondly, if instead of the normalization (3), we add
, where
is a fixed complex vector (see, for example, [2] ), (1) and
provide
complex equations for
complex unknowns, and the Jacobian of this new system is
The above Jacobian is square and can be easily shown to be nonsingular, using the ABCD Lemma [3] if the eigenvalue of interest is algebraically simple and
. Thirdly, if
is complex, then, as stated earlier, we have n complex and one real equation. Also, if
solves (4), then so does
for any
, such that
.
Our approach for analyzing the solution of (4) for
begins by splitting the eigenpair
into their real and imaginary parts:
,
where
, and
. After expanding (4), we will obtain a real
under-determined system of nonlinear equations in
real unknowns
, and it is natural to use the Gauss-Newton method (see, for example, Deuflhard ( [4] , pp. 222-223)) to obtain a solution (see also, [5] [6] [7] [8] ). By linearizing the under-determined system of nonlinear equations, we obtain an under-determined system of linear equations involving the Jacobian. This paper is structured as follows: in Section 2, we show that the Jacobian has a unique nullvector at the root. This is then followed in Section 3, we present two orthogonality results and Algorithm 1 is given. In Section 4, we present an inexact inverse iteration algorithm with preconditioning for solving the large system of equations encountered.
Algorithm 1. Computing the complex eigenvalues of the pencil
.
The main mathematical tools used in this paper are the LU factorization, inexact inverse iteration and preconditioned Generalized Minimal Residual (GMRES) [9] . The main reason for using inexact inverse iteration is due to the fact that as mentioned in an earlier paper, we do not solve
, in practice but we will use it to show that the solution is possible. Theorem 2.1 shows that the Jacobian has a single nullvector at the root, while Theorem 3.1 gives an important orthogonality result. Algorithms 1-3 are presented. We remark that in the limit, the approximate nullvector converges to the exact. A numerical example is given which supports the validity of the algorithms presented though, as usual, relies on good initial guesses to the desired eigenpair. The classical inverse iteration for the matrix pencil converges slowly for some eigenvalue problems while we present algorithms that converge quadratically. Throughout this paper, unless otherwise stated all norms are the 2-norm.
The following result helps to enforce the validity of the results in this paper.
Lemma 1.1: [10] Let
be of full rank. If
, is an under-determined linear system of equations, then its least squares solution
, is orthogonal to the nullspace of
.
Algorithm 2. Inexact inverse iteration algorithm.
Algorithm 3. Complex eigenvalues of the pencil
using Inexact Inverse Iteration with preconditioning.
In the next section, we will express both
and
as
and
, convert the nonlinear system (4) to a real under-determined system of nonlinear equations and prove some important results.
2. Computation of Complex Eigenpairs by Solving an Under-Determined System of Nonlinear Equations
In this section, we will expand the system of nonlinear n complex and one real equations in
complex unknowns (4) by writing
and
as
and
, respectively. The reason for having an underdetermined system of equations instead of a square system of equations is because, expanding
gives only one real equation, since
is symmetric positive definite, while
results in 2n real equations. This results in a
real under-determined system of nonlinear equations in
real unknowns. This will then be followed by presenting the underdetermined system of nonlinear equations and explicit expression for its Jacobian. Furthermore, we will show in the main result of this section-Theorem 2.1 that if the eigenvalue of interest in
is algebraically simple, then the Jacobian has linearly independent rows. We will find the right nullvector of the Jacobian at the root and proof that it is unique.
If we let
and
, then the square nonlinear system of Equations (4) can be written as
(6)
and
(7)
Hence,
Since
, we equate the real and imaginary parts of (6) to zero and obtain the 2n real equations
and
. This means,
consists of the 2n real equations
arising from (6) and one real equation
;
(8)
where
. The Jacobian
of
which has the following explicit expression
(9)
is a
by
real matrix. From the Jacobian (9) above, we define the real 2n by 2n matrix
as
(10)
Also, we form the 2n by 2 real matrix
(11)
consisting of the product of
and the matrix of right nullvectors of
at the root, where
(12)
and
is the n by n zero matrix. The Jacobian (9) can be rewritten in the following partitioned form
(13)
with
,
are as defined in (10) and (11) respectively. Note that because at the root,
this implies that
or its nonzero scalar multiple is a right nullvector of
. In the same vein, we find
and
or its nonzero scalar multiple is also a right nullvector of
at the
root. Since the eigenvalue
of
is algebraically simple by assumption, then by (2), we need to give explicit expressions for the left nullvector of
in order to prove that the Jacobian has full row rank at the root. Observe that if we define
, where
for all
, then this implies
Hence,
and
. The implication of this is that
which means,
or its nonzero scalar multiple is a left nullvector of
. Similarly,
and it shows that
is also a left nullvector of
.
So we form the matrix
consisting of the 2-dimensional left nullvectors of
at the root (in practice
is not computed), as
(14)
Now, observe that the condition (2), implies
Therefore, at the root, either
or
;
or
or both
and
are nonzero. It excludes the possibility that they are both zero.
Before we continue with the rest of the analysis, we pause a little to present the main result of this section which shows that the Jacobian (9) has a one dimensional nullvector at the root.
Theorem 2.1 Assume that the eigenpairs
of the pencil
are algebraically simple. If
and
are nonzero vectors, then
,
is the unique nonzero nullvector of
at the root.
Proof: See [5] .
After linearizing
, we have the following under-determined linear system of equations
(15)
Let
be the exact nullvector of the Jacobian
. By adding the extra condition
, which stems from Lemma 1.1 to the underdetermined linear system of Equations (15), we obtain the following square linear system of equations
(16)
3. Square System of Equations for the Numerical Computation of the Complex Eigenvalues of a Matrix
In the preceding section, we saw that by adding an extra equation to the under-determined linear system of equations 15, we obtained a square system of equations (16). However, in practice we would never compute
, but Theorem 2.1 guarantees the existence of a unique nullvector
at the root. This is the motivation for the discussion in this section. In this section, we will use
defined by
as an approximation to
in (16) and show that the solution obtained by solving (15) is equivalent to the solution obtained by solving
(17)
in the absence of round off errors. To do this, we will show that
for each k, and this is presented in the main result of this section: Theorem 3.1. Algorithm 1 is given for computing an algebraically simple eigenpair of the pencil
. Note that since
has been shown to be singular at the root in section 2, this section is anchored on the assumption that when
is not at the root,
is nonsingular.
First, we define the 2n by 2n matrix
as (see, also [11] )
(18)
and
(19)
The matrix
satisfies the following properties:
1)
.
2)
, where
is the 2n by 2n identity matrix.
3)
4) The matrix
commutes with
, i.e.,
.
5) For
,
.
6) Let
be an unknown vector that solves
. By premultiplying both sides by
we obtain
and hence
because of the commutativity of
and
. Therefore, if
, then
solves
.
We begin by writing the linear system of Equations (15) explicitly. For ease of notation, we shall drop the superscripts and define
where
, and replace
and
with
and
respectively. This means that (15) can now be rewritten as:
(20)
which is equivalent to the following system of equations
After rearrangement, the first n-equation reduces to
(21)
By multiplying both sides of the
equation by 2, we obtain:
This in turn reduces to
(22)
Observe that since
,
and
. Now,
. Consequently,
(23)
The combined set of Equations (21) and (23), which is the simplified form of (20), can be expressed as:
(24)
We assume that the 2n by 2n matrix
is nonsingular except at the root. This is what forms the basis for the following discussion. That is to say, we want to show that when not at the root,
.
First of all, let the exact nullvector
of
be defined as
, where
,
are real scalars,
and
are defined respectively by (19) and (10). Hence,
then after expanding the matrix-vector multiplication, we obtain
(25)
(26)
We make distinctly clear at this juncture, that the nullvector
is not exactly the same as
because, the later has the form of the exact nullvector at the root, but is evaluated at the kth iterate while the former is the nullvector even when not at the root.
Another way of writing (24) is as follows
(27)
This means that we could solve (24) by solving
(28)
for
. After which the solution of (27) is given by
(29)
With this expression for
, it can be observed that
Which means that
is well defined. Furthermore, from (25)
using the fact that
commutes with
and (28) gives
(30)
Since
is
-orthogonal to
by virtue of Equation (26), taking the
-inner product of both sides of the above with
yields
From which we deduce
(31)
Consider the problem of solving the under-determined linear system of Equations (20) for the
real unknowns
. It was stated in Lemma 1.1 that the minimum norm solution to an under-determined linear system of equations is orthogonal to the nullspace. It is an application of this result that yields the following important relationship:
(32)
If we add the nullvector
to the last row of (24), then
(33)
By expanding the second to the last row,
. But from
(29),
. This implies that, by taking the inner product of both sides with
, yields
Using the definition (31) for
and
, we obtain
(34)
where the unknown quantities
and
are to be determined, so we need an extra equation to be able to do so. The last row of the matrix-vector multiplication (cf. (33)) above comes from (32) since
If we substitute the expression (30) for
and (29) for
into the left hand side, then one obtains
Furthermore, by expanding the first term on the left hand side, using the properties of
, then
Consequently,
Observe that because
is real,
is nonzero. Accordingly, after dividing both sides by
(35)
We combine the two equations (34) and (35) below
(36)
to compute
and
simultaneously. The matrix on the left hand side is always nonsingular except at the root (in which case all entries are zero), this is because, its determinant is
. Equation (35) can now be applied to simplify
(37)
Notice that we have used the property
to arrive at the second step above and the definition (29) for
.
Next, we want to establish the orthogonality of
and
in the next key result. Before we do that, notice from Theorem (2.1) that
, at the root, is a scalar multiple of
and by the definition of
in (18), we can also write
, with
. This result holds when
is at the root or not, but because the analysis used to establish the orthogonality is based on the assumption that
is nonsinsingular when not at the root. As a result of this, after presenting Algorithm (1) (to follow shortly), we then show that the same result holds when
is singular at the root.
Theorem 3.1 Let
be an approximation to the exact nullvector
of
. Then,
is orthogonal to
.
Proof: To proof this, recall that
,
and
. This implies
(38)
showing that
. In arriving at the last step above, we have used the properties of
and a special case of (37) where
.
We present Algorithm (1), which involves the solution of two linear systems. The first is the 2n by 2n linear system of equations in (28), while the second is the 2 by 2 linear system (36).
Stop Algorithm 1 as soon as
where
. The above analysis shows that
when
is not at the root. Next, we want to show that the same result holds at the root.
In a manner analogous to the proof of Lemma (2.1), we postmultiply both sides of the above system of equation by
where
is the 2n by 2 real matrix defined by (14), consisting of the left nullvectors of
. If
and
are as defined respectively in (10) and (11), then, this is the same as
But by the definition of
, the first term on the left hand side of the equation above is zero, since
. It can be recalled from the proof of Lemma 2.1 that the 2 by 2 real matrix
is nonsingular at the root. This implies
Accordingly,
because of the nonsingularity of
. Therefore,
(39)
From the property of
at the root, it has two nonzero nullvectors and hence singular. The implication of this fact and the above is that,
is in the nullspace of
. But, we have already established that the nullspace of
consists of
and
. Hence, we can write
Now, from the last equation in (24),
Consequently,
But at the root,
. Therefore,
,
,
and
. This will now be used to deduce the following corollary at the root.
Corollary 3.1 Let
. Let
. Then,
is orthogonal to
at the root.
Proof: This follows from the second to the last line of the proof of Theorem 3.1 (cf., Equation (38)) where
. Hence,
4. Inexact Inverse Iteration with Preconditioning for Solving (28)
In Section 2, we found two nonzero nullvectors for
at the root. As a result of this property of
at the root, in this section, we will describe an inexact inverse iteration technique for solving the large sparse system of Equations (28) in step 2 of Algorithm 1 and present Algorithm 2 and Algorithm 3. Result of a numerical experiment is given which supports the theory in Section 5.
We give the following version of inexact inverse iteration in Algorithm 2. We will use a fixed tolerance. Note that because of the special nature of
at the root, the choice of a preconditioner is crucial for convergence to the desired accuracy to be achieved. The inexact linear solver that we use is the preconditioned GMRES [9] where we define the following block tridiagonal preconditioner,
(40)
Next, we present Algorithm 3, which is the inexact inverse iteration equivalence of Algorithm 1. The stopping criterion for the outer iteration in Algorithm 3 depends on the norm of the eigenvalue residuals, that is
and
and
5. Numerical Experiments
As mentioned earlier, the sparse linear system of equations in step 2 of Algorithm 1 is solved with an LU type factorization of
which is expensive and besides the L and U factors may have more nonzero elements than
. In this section, our main goal is to use preconditioned GMRES with the block triangular preconditioner
in (40) to solve the system
inexactly. We do this by considering a single numerical experiment with a fixed and a decreasing tolerance. Results are presented by means of tables and figures.
Example 5.1
Consider the 200 by 200 matrix
bwm200.mtx from the matrix market library [12] . It is the discretised Jacobian of the Brusselator wave model for a chemical reaction. The resulting eigenvalue problem with
was also studied in [13] and we are interested in finding the rightmost eigenvalue of
which is closest to the imaginary axis and its corresponding eigenvector.
In this example, we take
in line with [13] and took
and
, where
is the vector of all ones.
In Table 1 and Table 2, we present results for the computation of the eigenpairs
for the matrix in Example 5.1, and stop once the norms of
and the eigenvalue residuals
and
are smaller than the outer tolerance
. f represents the number of inner iterations used by preconditioned GMRES to satisfy the fixed inner tolerance
in Table 1, while d in Table 2 represents number of inner iterations used by preconditioned GMRES to satisfy the decreasing inner tolerance
. Quadratic convergence to
is easily observed from the second up to the seventh iterates in columns five and six of Table 2. However, this quadratic convergence is lost in the last iterate. This could be due to the fact that we are solving a nearly singular system as we approach the root. As shown in columns five and six of Table 1, as we approach the root, more number of inner iterations
Table 1. Table showing convergence to the eigenvalue
with a fixed inner tolerance for Example 5.1. The last column shows the number of inner iterations it took to satisfy the fixed inner tolerance
.
Table 2. Table showing quadratic convergence to the eigenvalue
with a decreasing tolerance for Example 5.1. The last column shows the number of inner iterations it took to satisfy the decreasing inner tolerance
.
Figure 1. Convergence history of the eigenvalue residuals on Example 5.1 with a fixed tolerance
.
were needed to satisfy the decreasing inner tolerance as against the fixed tolerance. From Table 1, we observed superlinear convergence in columns five and six.
Figure 1 and Figure 2 shows a plot of the norm of the eigenvalue residuals against the outer iterations with fixed and decreasing inner tolerances respectively. While the norm of the eigenvalue residuals decayed almost
Figure 2. Convergence history of the eigenvalue residuals on Example 5.1 with a decreasing tolerance
.
superlinearly in Figure 1, we observed a quadratic decrease in the norm of the eigenvalue residuals in Figure 2. It is quite surprising to see that Algorithm 3 works because
is singular at the root, which means we solved a singular system at the root.
6. Conclusions
While Ruhe ( [2] Section 3) used the normalisation
and solved the resulting
by
nonlinear system of equations to obtain
. We have been able to show that, with the addition of the non differentiable normalisation
, it is still possible to convert the resulting system of under-determined nonlinear equations into a square one.
Nevertheless, in this work, we obtained Algorithm 1 which consists of a combination of a 2n -by-2n system of equations (which is the same as those in [13] ) and a 2-by-2 system. A numerical example show that using an LU-solver on the one hand and preconditioned GMRES as an inexact solver on the other hand to solve the large sparse 2n-by-2n system of Equations (28), followed by the 2 by 2 system in each case, give similar results in the limit. By and large, the algorithms presented in this paper relies on good initial guesses to the desired eigenpair of interest.
Acknowledgements
The first author acknowledges funds provided by the University of Bath, United Kingdom during his PhD as well as an anonymous referee for useful comments.