Journal of Applied Mathematics and Physics
Vol.06 No.02(2018), Article ID:82705,17 pages
10.4236/jamp.2018.62040
A Block-Preconditioned Inexact Linear Solver for Computing the Complex Eigenpairs of a Large Sparse Matrix
Richard Olatokunbo Akinola1*, Stephen Yakubu Kutchin1, Ayodeji Sunday Ayodele1, Kingsley Obiajulu Muka2
1Department of Mathematics, Faculty of Natural Sciences, University of Jos, Jos, Nigeria
2Department of Mathematics, University of Benin, Benin, Nigeria
Copyright © 2018 by authors and Scientific Research Publishing Inc.
This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).
http://creativecommons.org/licenses/by/4.0/
Received: December 21, 2017; Accepted: February 25, 2018; Published: February 28, 2018
ABSTRACT
In computing eigenpairs of the matrix pencil, one obtains a linear system of equations. In this work, we show how a block triadiagonal preconditioner in GMRES can be used to solve a large sparse unsymmetric system of equations inexactly using a fixed and decreasing tolerance. While the fixed tolerance solver converged superlinearly to the eigenvalue of interest, the decreasing one converged quadratically. This surpasses an earlier result which converged harmonically.
Keywords:
Preconditioner, Superlinear Convergence, Quadratic Convergence, Eigenvalues
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.
Cite this paper
Akinola, R.O., Kutchin, S.Y., Ayodele, A.S. and Muka, K.O. (2018) A Block-Preconditioned Inexact Linear Solver for Computing the Complex Eigenpairs of a Large Sparse Matrix. Journal of Applied Mathematics and Physics, 6, 429-445. https://doi.org/10.4236/jamp.2018.62040
References
- 1. Stewart, G.W. (2001) Matrix Algorithms. Vol. II, Eigensystems, SIAM.
- 2. Ruhe, A. (1973) Algorithms for the Nonlinear Eigenvalue Problem. SIAM Journal on Matrix Analysis and Applications, 10, 674-689. https://doi.org/10.1137/0710059
- 3. Keller, H.B. (1977) Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. In: Rabinowitz, P., Ed., Applications of Bifurcation Theory, Academic Press, New York, 359-384.
- 4. Deuflhard, P. (2004) Newton Methods for Nonlinear Problems. Springer, 174-175.
- 5. Akinola, R.O. (2014) Theoretical Expression for the Nullvector of the Jacobian: Inverse Iteration with a Complex Shift. International Journal of Innovation in Science and Mathematics, 2, 367-371.
- 6. Akinola, R.O. and Spence, A. (2014) Two-Norm Normalization for the Matrix Pencil: Inverse Iteration with a Complex Shift. International Journal of Innovation in Science and Mathematics, 2, 435-439.
- 7. Akinola, R.O. and Spence, A. (2015) Numerical Computation of the Complex Eigenvalues of a Matrix by solving a Square System of Equations. Journal of Natural Sciences Research, 5, 144-157.
- 8. Akinola, R.O. (2015) Computing the Complex Eigenpair of a Large Sparse Matrix in Complex Arithmetic. International Journal of Pure & Engineering Mathematics (IJPEM), 3, 137-158.
- 9. Saad, Y. and Schultz, M.H. (1986) GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7, 856-869. https://doi.org/10.1137/0907058
- 10. Boyd, S. (2008/2009) Lecture 8: Least Norm Solutions of Under-Determined Equations, EE263 Autumn.
- 11. Freitag, M.A. and Spence, A. (2009) The Calculation of the Distance to Instability by the Computation of a Jordan Block, Linear and Nonlinear Eigen Problems for PDEs, 274-275.
- 12. Biosvert, B., Pozo, R., Remington, K., Miller, B. and Lipman, R. Matrix Market. http://math.nist.gov/MatrixMarket/
- 13. Parlett, B.N. and Saad, Y. (1987) Complex Shift and Invert Strategies for Real Matrices. Linear Algebra and Its Applications, 575-595.