Some numerical aspects in the approximation of eigenpairs of matrices by the Newton method

SOME NUMERICAL ASPECTS IN THE APPROXIMATION OF EIGENPAIRS OF MATRICES BY THE NEWTON METHOD

E. CǍTINAŞ, I. PǍVǍLOIU

Abstract. The eigenpairs of matrices may be approximated using the Newton method. The nonlinear mapping in such a case may be chosen in different ways, depending on the norming function. We investigate the difficulties of solving the linear systems from such iterations by evaluating the condition number at the solution for each of the resulted Jacobians in the case of two test matrices.

1. The eigenpairs of matrices as zeroes of nonlinear mappings

Denote $V={𝕂}^n$ and let $A=\left(a_{ij}\right)\in {𝕂}^{n\times n}$, where $𝕂=ℝ$ or $ℂ$. We recall that the scalar $\lambda \in 𝕂$ is an eigenvalue of $A$ if there exists $v\in V$, $v\ne 0$ such that

The vector $v$ is called an eigenvector corresponding to the eigenvalue $\lambda$. Since for an eigenvalue $\lambda$ the eigenpair $(v,\lambda )$ is not uniquely determined, it is necessary to impose a supplementary condition. As in [9] we shall consider a ”norming” function $G:V\to 𝕂,$ $G\left(0\right)\ne 1$ and besides (1) the equation

The function $G$ may be chosen in different ways:

where $\parallel \cdot {\parallel }_2$ stands for the euclidean norm.

In the paper [2] we proposed the choice

Considering ${𝕂}^{n+1}=V\times 𝕂$ endowed with the $\max$ norm $\parallel \cdot {\parallel }_{\mathrm{\infty }}$, we are led to the nonlinear mapping $F:{ℝ}^{n+1}\to {ℝ}^{n+1},$

Its first $n$ components are

The last component is determined by the choice of $G$. We get for each corresponding case:

It can be seen now that the eigenpairs of $A$ may be regarded as the solutions of the nonlinear system

For the solving of such a problem we may use the Newton method, which requires the derivative of $F$. It is known that the first and the second order derivatives of $F$ are given by

for all $x_0=\left(\genfrac{}{}{0.0pt}{}{v_0}{{\lambda }_0}\right)$, $h=\left(\genfrac{}{}{0.0pt}{}{u}{\alpha }\right)$, $k=\left(\genfrac{}{}{0.0pt}{}{w}{\beta }\right)\in {𝕂}^{n+1}$, where $e_i$ is the $i$-th vector from the canonical basis.

Some direct computations show that the norms of the second derivatives of $F$ are given by

regardless of $x\in {ℝ}^{n+1}$, since $F$ is a polynomial of degree $2$.

Different Newton-type methods for solving such nonlinear systems were studied in several papers [9], [6], [12], [14], [1], [2] and [3]. It is interesting to note that the Rayleigh quotient method is equivalent with a certain Newton method (”the scaled Newton method” [7]).

2. THE ROLE OF CERTAIN QUANTITIES IN THE CONVERGENCE OF THE NEWTON METHODS

We are interested how the different quantities specific to a nonlinear mapping influence the convergence of the Newton iterates.

The following result was obtained by Deuflhard and Potra

Theorem 1: [5] Let $F:D\subseteq {𝕂}^{n+1}\to {𝕂}^{n+1}$ be a nonlinear operator defined on the convex domain $D$. Suppose that $F$ is Fréchet differentiable on $D$, that $F^{\prime }\left(x\right)$ invertible for all $x\in D$, and that

$\forall x,y\in D$, $t\in [0,1]$. Moreover, assume that $x^{\ast }\in D$ is a solution of the nonlinear system $F\left(x\right)=0$. Let $x_0\in D$ be such that ${\overline{B}}_{\Vert x_0-x^{\ast }\Vert }\left(x^{\ast }\right)=\{x\in {𝕂}^{n+1}|\Vert x-x^{\ast }\Vert \le \Vert x_0-x^{\ast }\Vert \}\subset D$ and $x_0\in B_{2/\omega }\left(x^{\ast }\right)$. Then the Newton iterates

remain in the open ball $B_{\Vert x_0-x^{\ast }\Vert }\left(x^{\ast }\right)$, converge to $x^{\ast }$ and satisfy

Similar but somehow weaker results were obtained by Rheinboldt [10], Dennis and Schnabel [4] and by Ortega and Rheinboldt [8].

Remark. If $\Vert F^{\prime }{\left(x\right)}^{-1}\Vert \le \beta$ and $\Vert F^{\prime \prime }\left(x\right)\Vert \le K$ for all $x\in D$, one can take $\omega =\beta K$.

As explicitly noted by Dennis and Schnabel [4], the smaller are $\beta$ and $K$, the faster is the convergence of the iterates. On the other hand, as $\beta$ and $K$ decrease, the radius of the attraction ball becomes larger.

We see that in the case of large matrices, the choice III yields a much smaller constant $K$ compared to choice II (the choice I also yields small $K$, but its use may theoretically not be safe [2]). Of course that taking the norming function III, the constant $\beta$ may grow compared to the corresponding one from II but we do not analyse here this aspect.

There is another aspect important in the practical use of the Newton method. When the iterates converge, the matrices tend to $F^{\prime }\left(x^{\ast }\right)$. The condition number of a nonsingular matrix $A$, defined as $\kappa \left(A\right)=\Vert A\Vert \Vert A^{-1}\Vert$, is a measure of the sensitivity of a linear system $Ax=b$ to perturbations. The larger is $\kappa \left(A\right)$, the more inaccurate results one may obtain when using direct methods.

The Newton methods for eigenproblems were also considered in [13], in the context of preconditioning techniques.

The following section will be devoted to the numerical study of the condition numbers of two test matrices.

NUMERICAL ASPECTS

We shall consider two test matrices (these matrices are available from MatrixMarket at the following address: http://math.nist.gov/MatrixMarket/) in order to study the condition numbers. The programs were written in Matlab 5.0 and were run on a PC.

Pores1 matrix. This matrix arise from oil reservoir simulation. It is real, unsymmetric, of dimension 30 and has 20 real eigenvalues. We have plotted in fig. 1 the evaluations of the condition numbers of $F^{\prime }$ at the real solutions (the eigenpairs were computed by Matlab).

One can see that in this case the condition numbers are very close. At the nineth eigenvalue, the condition number of $F^{\prime }\left(x^{\ast }\right)$ is slightly better conditioned with choice III than with choice II.

Fidap002 matrix. This real symmetric matrix of dimension $n=441$ arise from finite element modeling. Its eigenvalues are all simple and range from $-7\cdot {10}^8$ to $3\cdot {10}^6$. Figure 2 contains the plot of the difference between the condition numbers computed with choices II and III.

One can see that in this case the differences that appear between the two condition numbers are on the whole small, but for some certain few eigenpairs the differences are large in favour of choice II. However, we cannot conclude that one of the choices has better overall condition numbers.

REFERENCES

1. E. Cǎtinaş and I. Pǎvǎloiu, On the Chebyshev Method for Approximating the Eigenvalues of Linear Operators, Rev. Anal. Numér. Théor. Approx., 25 (1996) no. 1-2, pp. 43-56.

2. E. Cǎtinaş and I. Pǎvǎloiu, On a Chebyshev-Type Method for Approximating the Solutions of Polynomial Operators Equations of Degree 2, Approximation and Optimization. Proceedings of the International Conference on Approximation and Optimization, Romania (ICAOR), Cluj-Napoca, July 29 - August 1, 1996, vol. I, pp. 219-226.

3. E. Cǎtinaş and I. Pǎvǎloiu, On Approximating the Eigenvalues and Eigenvectors of Linear Continuous Operators, Rev. Anal. Numér. Théor. Approx., 26 (1997) no. 1-2, pp. 19-27.

4. J.E. Dennis, Jr. and R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice-Hall Series in Computational Mathematics, Englewood Cliffs, NJ, 1983.

5. P. Deuflhard and F.A. Potra, Asymptotic Mesh Independence of Newton-Galerkin Methods via a Refined Mysovskii Theorem, SIAM J. Numer. Anal., 29 (1992) no.5, pp. 1395-1412.

6. J.J. Dongarra, C.B. Moler and J.H. Wilkinson, Improving the Accuracy of the Computed Eigenvalues and Eigenvectors, SIAM J. Numer. Anal., 20 (1983) no. 1, pp. 23-45.

7. S.M. Grzegórski, On the Scaled Newton Method for the Symmetric Eigenvalue Problem, Computing, 45 (1990), pp. 277-282.

8. J.M. Ortega and W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.

9. G. Peters and J.H. Wilkinson, Inverse Iteration, Ill-Conditioned Equations and Newton’s Method, SIAM Review, 21 (1979) no. 3, pp. 339-360.

10. W.C. Rheinboldt, An Adaptive Continuation Process for Solving Systems of Nonlinear Equations, Polish Academy of Science, Banach Ctr. Publ., 3 (1977), pp. 129-142.

11. M.C. Santos, A Note on the Newton Iteration for the Algebraic Eigenvalue Problem, SIAM J. Matrix Anal. Appl., 9 (1988) no. 4, pp. 561-569.

12. R.A. Tapia and L.D. Whitley, The Projected Newton Method has Order $1+\sqrt{2}$ for the Symmetric Eigenvalue Problem, SIAM J. Numer. Anal., 25 (1988) no. 6, pp. 1376-1382.

13. K. Wu, Y. Saad and A. Stathopoulos, Inexact Newton Preconditioning Techniques for Eigenvalue Problems, Lawrence Berkeley National Laboratory report number 41382 and Minnesota Supercomputing Institute report number UMSI 98-10, 1998.

14. T. Yamamoto, Error Bounds for Computed Eigenvalues and Eigenvectors, Numer. Math. 34 (1980), pp. 189-199.

Aspecte numerice în aproximarea valorilor şi vectorilor proprii ai matricilor prin metoda lui Newton

Valorile şi vectorii proprii ai matricilor pot fi aproximaţ prin metoda lui Newton. Aplicaţia neliniarǎ în acest caz poate fi aleasǎ în mai multe moduri, depinzând de funcţia de normare. In prezenta notǎ investigǎm dificultǎţile de rezolvare a sistemelor liniare din astfel de iteraţii prin evaluarea numǎrului de condiţionare în soluţie pentru fiecare din Jacobienii rezultaţi în cazul unor matrici test.

Cerc. pr. gr. III Emil Cǎtinaş, ”T. Popoviciu” Institute of Numerical Analysis, P.O. Box 68, 3400 Cluj-Napoca, Romania (ecatinas@ictp.math.ubbcluj.ro)

Cerc. pr. gr. I Ion Pǎvǎloiu, ”T. Popoviciu” Institute of Numerical Analysis, P.O. Box 68, 3400 Cluj-Napoca, Romania (pavaloiu@ictp.math.ubbcluj.ro)