# A New Quotient and Iterative Detection Method in an Affine Krylov Subspace for Solving Eigenvalue Problems

Authors: Chein-Shan Liu; Jiang-Ren Chang; Jian-Hung Shen; Yung-Wei Chen
Journal: Journal of Mathematics
 (2023)

---

## Abstract

The paper solves the eigenvalues of a symmetric matrix by using three novel algorithms developed in them-dimensional affine Krylov subspace. The n-dimensional eigenvector is superposed by a constant shifting vector and an m-vector. In the first algorithm, the m-vector is derived as a function of eigenvalue by maximizing the Rayleigh quotient to generate the first characteristic equation, which, however, is not easy to determine the eigenvalues since its roots are not of simple ones, exhibiting turning points, spikes, and even no intersecting point to the zero line. To overcome that difficulty by the first algorithm, we propose the second characteristic equation through a new quotient with the inner product of the shifting vector to the eigen-equation. The Newton method and the fictitious time integration method are convergent very fast due to the simple roots of the second characteristic equation. For both symmetric and nonsymmetric eigenvalue problems solved by the third algorithm, we develop a simple iterative detection method to maximize the Euclidean norm of the eigenvector in terms of the eigen-parameter, of which the peaks of the response curve correspond to the eigenvalues. Through a few finer tunings to the smaller intervals sequentially, a very accurate eigenvalue and eigenvector can be obtained. The efficiency and accuracy of the proposed iterative algorithms are verified and compared to the Lanczos algorithm and the Rayleigh quotient iteration method.


---

## Body

## 1. Introduction

LetA=aij be an n×n real symmetric matrix. For a nonzero vector x=x1,…,xnT,(1)Rx≔=xTAxxTx,is the Rayleigh quotient [1, 2]. If x is an eigenvector of A, then the corresponding Rayleigh quotient is an eigenvalue λ, satisfying the vector form eigen-equation:(2)Ax=λx,where λ is an eigen-parameter and nontrivial solutions with x≠0 exist only for λ being an eigenvalue; otherwise, x=0.In the free vibrations of a multidegree mass-damping-spring structural system, the second-order ordinary differential equations system is usually transformed to the Rayleigh quotient to find the natural frequencies [3]. The Rayleigh quotient has important roles in physical applications, and the eigenvalue problems are ubiquitous to be the main mathematical tools to explain many physical phenomena [4, 5], such as vibration frequencies, stability of dynamical systems, and energy levels in quantum mechanics, and also in the Mathieu’s equations for describing the wave propagation, the electromagnetic, and the elastic membrane as well as the heat conduction [6].The Rayleigh quotient (1) can be used to find λ of a symmetric matrix eigenvalue problem (2). The related systems which lead to eigenvalue problems were discussed by Tisseur and Meerbergen [7], and different numerical methods to solve the eigenvalue problems were depicted in [8–10].The basic ingredient of the Rayleigh–Ritz method is the construction of an orthogonaln×m matrix V, and set x=Vy; such that the original n-dimensional eigenvalue problem is reduced to an eigenvalue problem in an m-dimensional subspace with Bmy=μy and Bm≔=VTAV, where μ,y is a Ritz pair to approximate the true eigenvalue and eigenvector. The computational cost in the m-dimensional subspace can be significantly reduced if m≪n.The Lanczos method is further reducedBm to a symmetric tridiagonal matrix. A lot of works were concerned with the Lanczos algorithm [11–14]. However, it is restricted to finding at most m dominant eigenvalues of A. Bai [15] gave an error analysis of the Lanczos algorithm for the nonsymmetric eigenvalue problem.The numerical computations in [16, 17] revealed that the methods based on the Krylov subspace could be very effective in the nonsymmetric eigenvalue problems by using the Lanczos biorthogonalization algorithm and Arnoldi’s algorithm. The Arnoldi and nonsymmetric Lanczos methods are both of the Krylov subspace methods. Among the many algorithms to solve the matrix eigenvalue problems, the Arnoldi method [17–20], the nonsymmetric Lanczos algorithm [21], and the subspace iteration method [22] are well known.The present paper will introduce some novel methods to find all eigenvalues ofA in an affine Krylov subspace. Two new characteristic equations are derived to solve λ, of which both the eigenvector x and the eigenvalue λ can be solved simultaneously. Moreover, we will develop a very powerful iterative detection method to solve nonsymmetric eigenvalue problems in the affine Krylov subspace using the Galerkin condition to derive the governing equation of the m-vector. The affine Krylov subspace method was first developed by Liu [23] to solve a linear equations system, which is, however, not yet used to solve the matrix eigenvalue problems. There are many applications of the eigenvalues in PageRank [24, 25] and in free vibration [26].Subsequently, we express the Rayleigh quotient in an affine Krylov subspace in Section2, where a main result is derived to present the n-dimensional eigenvector in terms of a constant shifting vector and an m-vector governed by a nonhomogeneous linear equations system. The numerical algorithm for a symmetric matrix is developed in Section 3, where we derive the first characteristic equation. To improve the property of the first characteristic equation, we propose a new quotient in Section 4, where the second characteristic equation is derived for the symmetric eigenvalue problem. For the sake of comparison, we describe the Lanczos algorithm and the Rayleigh quotient iteration method in Section 5. The numerical tests of symmetric eigenvalue problems are carried out in Section 6. In Section 7, we discuss the properties of these two derived characteristic equations. In Section 8, we propose an iterative detection method for the nonsymmetric eigenvalue problems by maximizing the Euclidean norm of the derived eigenvector in the affine Krylov subspace, which is very powerful. Finally, we highlight the major results in Section 9.


## 2. The Rayleigh Quotient in an Affine Krylov Subspace

### 2.1. Affine Krylov Subspace

Given a nonzero constant shifting vectorx0∈n, the m-dimensional Krylov subspace is generated by(3)Km≔SpanAx0,…,Amx0.We employ Arnoldi’s process (a modification of the Gram–Schmidt process) to orthogonalize and normalize the Krylov vectorsAjx0, j=1,…,m, whose resultant vectors ui, i=1,…,m satisfy ui⋅uj=δij, i,j=1,…,m, where δij is the Kronecker delta symbol. After giving x0 and the dimension m, the Arnoldi procedure is written as follows [27]:(4)u1=Ax0Ax0,Doj=1:m−1,wj=Auj,Doi=1:j,hij=ui⋅wj,wj=wj−hijui,Enddo ofi,hj+1,j=wj.Ifhj+1,j=0stop,uj+1=wjhj+1,j,Enddo ofj.LetU denote the Krylov matrix with dimension n×m(5)U≔u1,…,um,which possesses the orthogonal property(6)UTU=Im.


### 2.2. A Main Theorem

Theorem 1.
Forx∈x0+Km, x0≠0, the eigenvector x for(7)maxxR=xTAxxTx,is given by(8)x=x0+UB−λIm−1UTλx0−Ax0,where B≔UTAU and the eigenvalue λ is to be solved from a nonlinear characteristic equation.Proof.
Because ofx∈x0+Km, we can expand x by(9)x=x0+Uα,where the coefficient α∈m is to be determined. By using equations (6) and (9), we can derive(10)μ1≔xTAx=x0TAx0+2x0TAUα+αTBα,(11)v1≔∇αμ1=2UTAx0+2Bα,(12)μ2≔x2=x02+2x0TUα+αTα,(13)v2≔∇αμ2=2UTx0+2α,where(14)B≔UTAU,is a constant m×m matrix.
The maximality in equation (7) renders.(15)∇αxTAxxTx=0⇒x2∇αxTAx−xTAx∇αx2=0,which by equations (10)–(13) leads to(16)x2v1−xTAxv2=0.
It follows from equation (16) that v1 is proportional to v2, and as R defined by equation (1), we can write.(17)v1=λv2,where(18)λ=xTAxx2=μ1μ2,is an eigenvalue to be determined. By equations (11), (13), and (17), α is solved by(19)B−λImα=λUTx0−UTAx0,hence, α is a nonlinear m-vector function of λ. Then, inserting α into equation (9), we can achieve equation (8).
From equation (18), it follows that(20)Fλ≔λμ2μ1−1=0.
Upon inserting equation (19) for α into equations (10) and (12) and substituting μ1 and μ2 into the above equation, we can obtain a nonlinear characteristic equation to solve λ. This ends the proof.
Theorem1 indicates that the eigenvector x can be expressed in terms of x0 and the m-vectorα, which is a vector function of λ as shown by equation (8). Therefore, the current computational cost for solving the Rayleigh quotient and symmetric eigenvalue problems is very cheap by merely solving an m-dimensional linear system (19) to determine α and a nonlinear equation (20) to determine the eigenvalue λ, and the eigenvector x is computed by equation (9). We will give examples to show that the zero points of Fλ coincide to the roots obtained from the determinant-form characteristic equation as follows:(21)detA−λIn=0,however, it is not an efficient method to solve the eigenvalues when n is large.
Equation (20) is the first characteristic equation derived in the m-dimensional affine Krylov subspace, which can overcome the inefficiency of the determinant-form characteristic equation. However, the Rayleigh quotient used in equation (20) may cause a weak property of the resulting characteristic equation as to be demonstrated by examples given in Section 6.Remark 1.
We may consider a larger subspace with(22)K^m+1≔Spanx0,Ax0,…,Amx0,such that x0+KmK^m+1, and set(23)x=Vy,where V is constructed from K^m+1 by the Arnold procedure. As that carried out in the proof of Theorem 1, we can derive(24)μ1=xTAx=yTBm+1y,v1=∇yμ1=2Bm+1y,μ2=x2=y2,v2=∇yμ2=2y,where Bm+1≔=VTAV.
Then, by the maximality in equation (7), v1 is proportional to v2 with v1=μ1/μ2v2, and thus we come to an eigenvalue problem for Bm+1:(25)Bm+1y=λy.
Unlike equation (19) for α, we cannot solve y explicitly in terms of λ from equation (25). Therefore, we can conclude that x∈x0+Km instead of xK^m+1 is necessary, and expanding x by equation (9), rather than by equation (23), is a key way to succeed the presented new method, which is named the Rayleigh quotient method (RQM).


## 2.1. Affine Krylov Subspace

Given a nonzero constant shifting vectorx0∈n, the m-dimensional Krylov subspace is generated by(3)Km≔SpanAx0,…,Amx0.We employ Arnoldi’s process (a modification of the Gram–Schmidt process) to orthogonalize and normalize the Krylov vectorsAjx0, j=1,…,m, whose resultant vectors ui, i=1,…,m satisfy ui⋅uj=δij, i,j=1,…,m, where δij is the Kronecker delta symbol. After giving x0 and the dimension m, the Arnoldi procedure is written as follows [27]:(4)u1=Ax0Ax0,Doj=1:m−1,wj=Auj,Doi=1:j,hij=ui⋅wj,wj=wj−hijui,Enddo ofi,hj+1,j=wj.Ifhj+1,j=0stop,uj+1=wjhj+1,j,Enddo ofj.LetU denote the Krylov matrix with dimension n×m(5)U≔u1,…,um,which possesses the orthogonal property(6)UTU=Im.


## 2.2. A Main Theorem

Theorem 1.
Forx∈x0+Km, x0≠0, the eigenvector x for(7)maxxR=xTAxxTx,is given by(8)x=x0+UB−λIm−1UTλx0−Ax0,where B≔UTAU and the eigenvalue λ is to be solved from a nonlinear characteristic equation.Proof.
Because ofx∈x0+Km, we can expand x by(9)x=x0+Uα,where the coefficient α∈m is to be determined. By using equations (6) and (9), we can derive(10)μ1≔xTAx=x0TAx0+2x0TAUα+αTBα,(11)v1≔∇αμ1=2UTAx0+2Bα,(12)μ2≔x2=x02+2x0TUα+αTα,(13)v2≔∇αμ2=2UTx0+2α,where(14)B≔UTAU,is a constant m×m matrix.
The maximality in equation (7) renders.(15)∇αxTAxxTx=0⇒x2∇αxTAx−xTAx∇αx2=0,which by equations (10)–(13) leads to(16)x2v1−xTAxv2=0.
It follows from equation (16) that v1 is proportional to v2, and as R defined by equation (1), we can write.(17)v1=λv2,where(18)λ=xTAxx2=μ1μ2,is an eigenvalue to be determined. By equations (11), (13), and (17), α is solved by(19)B−λImα=λUTx0−UTAx0,hence, α is a nonlinear m-vector function of λ. Then, inserting α into equation (9), we can achieve equation (8).
From equation (18), it follows that(20)Fλ≔λμ2μ1−1=0.
Upon inserting equation (19) for α into equations (10) and (12) and substituting μ1 and μ2 into the above equation, we can obtain a nonlinear characteristic equation to solve λ. This ends the proof.
Theorem1 indicates that the eigenvector x can be expressed in terms of x0 and the m-vectorα, which is a vector function of λ as shown by equation (8). Therefore, the current computational cost for solving the Rayleigh quotient and symmetric eigenvalue problems is very cheap by merely solving an m-dimensional linear system (19) to determine α and a nonlinear equation (20) to determine the eigenvalue λ, and the eigenvector x is computed by equation (9). We will give examples to show that the zero points of Fλ coincide to the roots obtained from the determinant-form characteristic equation as follows:(21)detA−λIn=0,however, it is not an efficient method to solve the eigenvalues when n is large.
Equation (20) is the first characteristic equation derived in the m-dimensional affine Krylov subspace, which can overcome the inefficiency of the determinant-form characteristic equation. However, the Rayleigh quotient used in equation (20) may cause a weak property of the resulting characteristic equation as to be demonstrated by examples given in Section 6.Remark 1.
We may consider a larger subspace with(22)K^m+1≔Spanx0,Ax0,…,Amx0,such that x0+KmK^m+1, and set(23)x=Vy,where V is constructed from K^m+1 by the Arnold procedure. As that carried out in the proof of Theorem 1, we can derive(24)μ1=xTAx=yTBm+1y,v1=∇yμ1=2Bm+1y,μ2=x2=y2,v2=∇yμ2=2y,where Bm+1≔=VTAV.
Then, by the maximality in equation (7), v1 is proportional to v2 with v1=μ1/μ2v2, and thus we come to an eigenvalue problem for Bm+1:(25)Bm+1y=λy.
Unlike equation (19) for α, we cannot solve y explicitly in terms of λ from equation (25). Therefore, we can conclude that x∈x0+Km instead of xK^m+1 is necessary, and expanding x by equation (9), rather than by equation (23), is a key way to succeed the presented new method, which is named the Rayleigh quotient method (RQM).


## 3. Numerical Algorithm

For equation (20), Liu and Atluri [28] considered the variables transformation by(26)ηt=1+tλ,where t is a fictitious time variable. Equation (20) is equivalent to(27)0=−νFλ,ν≠0.Addingη′=dη/dt=λ. by equation (26) on both sides of equation (27)(28)η′=λ−νFλ,and using equation (26) for λ, yields(29)η′=η1+t−νFη1+t.Multiplying equation (29) by 1/1+t, obtains(30)dη/1+tdt=−ν1+tFη1+t.Usingη/1+t further, leads to(31)λ′t=−ν1+tFλt,which is a first-order ODE for λξ.The numerical procedures of the proposed algorithm RQM based on Theorem1 are summarized as follows:(i)
GiveA, select m, x0 and λ0, give v, Δt and ε, and construct U and B=UTAU(ii)
Fork=0,1,2,…, execute the following steps by using the fictitious time integration method (FTIM) [28] of equation (31)(32)λk+1=λk−vΔtFλk1+kΔt.For a quick convergence, we can also apply the Newton method to solve equation (20) as follows:(33)λk+1=λk−FλkF′λk,where(34)Fλk=λkμ2kμ1k−1,F′λk=μ2kμ1k+λkμ12kμ1kμ2′k−μ2kμ1′k,in which(35)αk=B−λkIm−1λkUTx0−UTAx0,μ1k=x0TAx0+2x0TAUαk+αkTBαk,μ2k=x02+2x0TUαk+αkTαk,v1k=2UTAx0+2UTAUαk,v2k=2UTx0+2αk,μ1′k=v1kTB−λkIm−1αk+UTx0,μ2′k=v2kTB−λkIm−1αk+UTx0.The scriptk denotes the kth step value. If Fλk satisfies Fλk<ε, then stop, and the eigenvalue λ and the eigenvector x=x0+Uα are obtained; otherwise, go to step (ii).When the Krylov matrixU is constructed, the computational cost of RQM is quite saving with the computations of four scalars μ1, μ2, μ1′, and μ2′, and three m-dimensional vectors α, v1, and v2 at each iteration step.


## 4. A New Quotient Method

In addition to equation (18), by taking the inner product of x0 to the eigen-equation (2), we have a simpler quotient to express the eigenvalue by(36)λ=μ3μ4,where(37)μ3≔x0TAx,μ4≔x0Tx.Under the conditionμ4=x0Tx≠0, equation (36) is well-defined. Upon comparing to μ1 and μ2 in equations (10) and (12) as being quadratic nonlinear functions of x, μ3 and μ4 in equation (37) are linear functions of x. Therefore, from equation (36), a simpler nonlinear equation is available as follows:(38)Gλ≔λμ4μ3−1=0,which is the second characteristic equation derived in the m-dimensional affine Krylov subspace.The numerical procedures based on the new quotient method are depicted as follows:(i)
GiveA, select m, x0 and λ0, give v, Δt and ε, and construct U and B=UTAU(ii)
Fork=0,1,2,…, carry out the following steps by using the FTIM(39)λk+1=λk−vΔtGλk1+kΔt,or by the Newton method:(40)λk+1=λk−GλkG′λk,whereGλk=λkμ4k/μ3k−1, G′λk=μ4k/μ3k+λk/μ32kμ3kμ4′k−μ4kμ3′k, in which(41)αk=B−λkIm−1λkUTx0−UTAx0,μ3k=x0TAx0+x0TAUαk,μ4k=x0Tx0+x0TUαk,μ3′k=x0TAUB−λkIm−1αk+UTx0,μ4′k=x0TUB−λkIm−1αk+UTx0.IfGλk<ε, then stop; otherwise, go to step (ii).It can be seen that the presented iterative algorithm resorted on the new quotient is simpler than that in Section3. We name this technique a new quotient method (NQM). When the Krylov matrix U is constructed, the computational cost of the NQM is quite saving with the computations of four scalars μ3, μ4, μ3′, and μ4′, and an m-dimensional vector α at each iteration step.


## 5. Lanczos Algorithm and Rayleigh Quotient Iteration Method

As pointed in [8], the Lanczos algorithm can produce a tridiagonal matrix T from the symmetric matrix A by(42)selectq1q1=1,α1=q1TAq1,β0q0=0,Dok=1:m−1,rk=A−αkInqk−βk−1qk−1,βk=rk,Ifβk=0stop,qk+1=rk/βk,αk+1=qk+1TAqk+1,Enddo ofk.LetQ denote the Krylov matrix with dimension n×m(43)Q≔q1,…,qm,which satisfies QTQ=Im. Then, T=QTAQ is an m×m tridiagonal matrix, given by(44)T=α1β10…β1α2⋱0…0⋱⋱βm−10…βm−1αm.Letpkλ be the characteristic function of a k×k submatrix of T. Starting from p0=1 and p1=α1−λ, we can apply the recurrence formula for k=2,…,m(45)pkλ=αk−λpk−1λ−βk−12pk−2λ,pk′λ=αk−λpk−1′λ−pk−1λ−βk−12pk−2′λ,where p0′=0 and p1′=−1, and then the Newton method to find the eigenvalue of λTλA is given by(46)λj+1=λj−pmλjpm′λj,j=0,…,which is starting from an initial guess λ0 and ending until λj+1−λj<ε.For the purpose of comparison, we list the Rayleigh quotient iteration method (RQIM) as follows [8]:(47)Selectx0x0=1,μ0=Rx0,Dok=0,…,SolveA−μkInzk+1=xk,xk+1=zk+1zk+1,μk+1=Rxk+1,Ifμk+1−μk<εstop.Remark 2.
Instead of solving anm-dimensional linear system (19) by the presented two new methods, the RQIM solves an n-dimensional linear system A−μkInzk+1=xk at each iteration. In the RQM and NQM, we give an initial guess of λ0 to seek an exact eigenvalue λ; but in the Rayleigh quotient iteration method, the initial guess is x0, which is more difficult than that given an initial guess λ0 to seek certain eigenvalue λ.Remark 3.
Although the methods in Sections3 and 4 and the Lanczos algorithm in Section 5 are performed in the m-dimensional (affine) Krylov subspace, for an n-dimensional symmetric eigenvalue problem, the first two methods can find n eigenvalues with a suitable value of m<n, while the Lanczos algorithm can find at most m eigenvalues, since the Lanczos algorithm is reduced the symmetric eigenvalue problem to find the eigenvalues of an m×m tridiagonal symmetric matrix.


## 6. Symmetric Eigenvalue Problems

Example 1.
We first demonstrate the usefulness of Theorem1 by(48)A=62−2250−207,which possesses the exact eigenvalues λA=3,6,9.
We use this simple example to compare the characteristic equations (20) and (38) to the determinant-form characteristic equation as follows:(49)detA−λI3=λ3−18λ2+99λ−162=0.
The parameters used in the RQM aren=3, m=2, and x0=1,2,3T. In Figure 1(a), Fλ is plotted with respect to the eigen-parameter λ∈2,10 by a solid line, where three remarked points correspond to the eigenvalues λ=3,6,9. Fλ traces across the zero line at the first point which is a simple root with Fλ=0 and F′λ≠0, but at the last two points λ=6 and λ=9, Fλ locates at the local minimal points where F′λ=0, which are not simple roots with turning points. The Newton method can be applied to find the first eigenvalue λ=3, but for the last two eigenvalues, it would be failed due to F′λ=0 at the root λ. When the Newton method is not applicable to find the last two eigenvalues λ=6 and λ=9, however, we can apply the FTIM to quickly find all eigenvalues.
In Figure1(b), det A−λIn is plotted with respect to λ∈2,10, where three intersecting points to the zero line correspond to the eigenvalues λ=3,6,9. The pattern of the characteristic curve obtained from det A−λIn is different from Fλ in Figure 1(a) by the solid line. The characteristic curve obtained from det A−λIn traces across those three points, which are all simple roots of det A−λIn=0.
Whenn is a large number, to realize the explicit form of the conventional characteristic equation det A−λIn=0 may be a tedious work and it is hard to employ the Newton method to solve the determinant-form characteristic equation to find all the eigenvalues.
But in the m-dimensional affine Krylov subspace withm≪n, the first characteristic equation Fλ=0 has the explicit form in equation (20), and its roots are easily found by using the FTIM. With v=300, Δt=0.01, ε=10−12, and λ0=0, the eigenvalue λ=3 is obtained by the FTIM through two steps, whose error to compare 3 is zero and the error of Ax−λx=1.79×10−14 is very small. With v=122, Δt=0.01, and λ0=7, the eigenvalue λ=6.00000000000404 is obtained by the FTIM through 92 steps, whose error to compare 6 is 4.04×10−12 and Ax−λx=6.89×10−11 is obtained. With v=300, Δt=0.01, and λ0=9.5, the eigenvalue λ=8.99999999999449 is obtained by the FTIM through 62 steps, whose error to compare 9 is 5.51×10−12 and Ax−λx=3.36×10−11 is obtained.
When we apply the Newton method to find the eigenvalue under a convergence criterionε=10−12 and with an initial guess λ0=1, we obtain λ=3.0000000000000001 through four steps, whose error to compare 3 is 4.44×10−16 and the error of Ax−λx=1.63×10−14. By starting from λ0=10, we obtain λ=9 through eight steps, whose error to compare 9 is zero and the error of Ax−λx=2.27×10−13. Obviously, when the Newton method is applicable, it converges faster and more accurate than the FTIM.
However, the Newton method with the Rayleigh quotient is failed to find the eigenvalueλ=6. In order to remedy this drawback, we employ the new quotient in equation (38). As shown in Figure 1(a) by the dashed line, at each root, the characteristic curve Gλ is crossing the zero line, which possesses the good property as that in equation (49) for the polynomial characteristic curve in Figure 1(b). It means that the roots of Gλ=0 are all simple roots with G′λ≠0. The property of the new quotient is better than the Rayleigh quotient used in equation (20). Then, we apply the NQM and the Newton method to solve the characteristic equation (38). By starting from λ0=0, through seven steps under ε=10−12, the error to compare 3 is 1.16×10−12 and Ax−λx=1.48×10−13 is obtained. By starting from λ0=7 and through five steps, the error to compare 6 is 8.08×10−14 and Ax−λx=5.72×10−13 is obtained. By starting from λ0=9.5 and through six steps, the error to compare 9 is 1.81×10−13 and Ax−λx=4.49×10−13 is obtained.Figure 1
For Example1, showing three roots of the nonlinear equations obtained from (a) the Rayleigh quotient and new quotient, and (b) characteristic equation corresponding to eigenvalues 3, 6, and 9.
(a)(b)Example 2.
We apply the new algorithms to solve equation (2) with(50)A=−1.64071.08141.20141.15391.08144.15737.4035−1.04631.20147.40352.7890−1.57371.1539−1.0463−1.57378.6944,where λA=−4,−2,8,12.
By takingn=4, m=3, and x0=14≔1,1,1,1T, we plot Fλ with respect to λ∈−4,12 in Figure 2(a), where four points correspond to the four eigenvalues λ=−4,−2,8,and12. It can be seen that there are several spikes and the quality of the first characteristic curve is not good. In contrast, we plot Gλ with respect to λ∈−4,12 in Figure 2(b) by a solid line, where four points obviously correspond to the four eigenvalues λ=−4,−2,8,and12. The improvement of the property from the spikes of Fλ to the crossing points with simple roots of Gλ is significant, which renders the work to find the eigenvalues quite easily.
We emphasize thatm must be large enough to find all eigenvalues. If we reduce m=3 to m=2, as shown in Figure 2(b) by a blue dashed line, the characteristic curve intersects to the zero line at only three points λ=−2,8,12. In this occasion, we cannot find the eigenvalue λ=−4 by using the proposed NQM with m=2. If one employs the Lanczos algorithm with m=2, only two eigenvalues can be found.
By using the Newton method to solveGλ=0 in equation (38) under a convergence criterion ε=10−12 and giving an initial guess λ0=13, the eigenvalue obtained is λ=12.0000500689830, which converges within ten steps. The error to compare 12 is 5.01×10−5 and the error of Ax−λx=1.08×10−10. Starting from λ0=7, the eigenvalue λ=7.99995888523901 is obtained with four steps. The error to compare 8 is 4.11×10−5 and the error of Ax−λx=5.8×10−10. Starting from λ0=−1, the eigenvalue λ=1.99993569548861 is obtained with six steps. The error to compare −2 is 6.43×10−5 and the error of Ax−λx=1.15×10−12. Starting from λ0=−5, the eigenvalue λ=−4.00007337793926 is obtained with six steps. The error to compare −4 is 7.34×10−5 and the error of Ax−λx=1.88×10−13.
On the other hand, we can solve this problem under a convergence criterionε=10−12, by applying the FTIM to solve Gλ=0 in equation (38). The parameters m=3, Δt=0.01, and ε=10−12 are taken. With v=−120 and λ0=13, the eigenvalue we obtain is λ=12.0000500689837 through 31 steps, whose error to compare 12 is 5.01×10−5. The error of Ax−λx=2.99×10−10. With v=150 and λ0=7, the eigenvalue obtained is λ=7.99995888524186 through 20 steps, whose error to compare 8 is 4.11×10−5 and the error of Ax−λx=4.21×10−10. With v=−250 and λ0=−1, the eigenvalue obtained is λ=−1.99993569548899 through 13 steps, whose error to compare −2 is 6.43×10−5 and the error of Ax−λx=3.81×10−12. With v=200 and λ0=−5, the eigenvalue obtained is λ=−4.00007337793938 through nine steps, whose error to compare 4 is 7.34×10−5. The error of Ax−λx=1.8×10−13. Overall, both the Newton method and the FTIM are convergent very fast and highly accurate in the solution of the eigenvalue problem with the new quotient in Gλ of equation (38) derived in the affine Krylov subspace.
Then, we apply the Lanczos algorithm together with the Newton method to determine the eigenvalues in Table1, where the parameters m and λ0 are listed, and the errors of eigenvalues are compared. When error 1 denotes the errors of eigenvalues obtained by the Lanczos algorithm, error 2 denotes the errors obtained by the present NQM. For the eigenvalues −4, 8, and 12, the present method is more accurate than the Lanczos algorithm. To obtain the eigenvalue λ=−2 by the Lanczos algorithm, m=4 was used.Figure 2
For Example2, showing four roots of the nonlinear equation obtained from (a) the Rayleigh quotient and (b) the new quotient, corresponding to eigenvalues −4, −2, 12, and 8 with m = 3, and −2, 12, and 8 with m = 2.
(a)(b)Table 1
For Example2, comparing the errors of eigenvalues solved by the Lanczos algorithm (error 1) and the NQM (error 2).
λ−4−2812λ0−5−2.511.512m3433Error 17.31×10−36.43×10−59.26×10−45.80×10−3Error 27.34×10−56.43×10−54.11×10−55.01×10−5Example 3.
LetA=aij, i,j=1,…,n. In equation (2), we take(51)aij=1i+j−1,i,j=1,…,n.
A is a Hilbert matrix with very large condition number.
We employ the Newton method together with the NQM in Section4 to find the eigenvalues. The parameters are n=5, m=4, and x0=15. In Figure 3(a), Gλ is plotted with respect to λ∈1.5,1.6, where the intersecting point to the zero line corresponds to an eigenvalue as a simple root of Gλ=0.
Starting fromλ0=2, the eigenvalue λ=1.56705069109793 is obtained with five steps under ε=10−12. The above eigenvalue is very close to 1.567050691098231 given in [29]. In contrast, when the Rayleigh quotient is employed in Fλ of equation (20), it gives no root in the interval λ∈1.5,1.6 as shown in Figure 3(b), which is incorrect.
Forn=6 and with m=4, 1.61889985892416 is obtained through five steps which is very close to 1.61889985892434 given in [29].
The smallest eigenvalue of the Hilbert matrix is difficult to be found. However, by using the NQM, we can do it easily. Forn=5 with m=4, starting from λ0=0, the eigenvalue λ=3.28792129597493×10−6 is obtained with 22 steps under ε=10−12, and the error of Ax−λx=3.47×10−14, which is very accurate. For n=6 with m=5, starting from λ0=0, the eigenvalue λ=1.0614278798386×10−7 is obtained with 26 steps under ε=10−10, and the error is Ax−λx=1.99×10−12.
In Table2, we list the largest eigenvalues for n=7,…,10, obtained by the Lanczos method and the NQM, which are convergent very fast with four or five steps by starting from the initial guess λ0=1.8. It can be seen that the presented results are very close to that given in [29] by using the cyclic Jacobi method [30]. Although, we solve the eigenvalues in a lower-dimensional subspace with dimension m=4, the results are very accurate.
To display the high efficiency of the new algorithm NQM, underε=10−10, we plot the number of iterations and the computed largest eigenvalues of the Hilbert matrices with n∈5,25 in Figure 4. Through 4 to 6 iterations the iterative algorithm is convergence as shown in Figure 4(a), and the largest eigenvalues increase with respect to n and tend to 1.96 as shown in Figure 4(b).Figure 3
For Example3, (a) showing one root of the nonlinear equation obtained from the new quotient, corresponding to the eigenvalue 1.56705069109793, and (b) obtained from the Rayleigh quotient giving no root in the interval which is incorrect.
(a)(b)Table 2
The largest eigenvalues of the Hilbert matrices withn = 7, 8, 9, and 10.
(n, m)(7, 4)(8, 4)(9, 4)(10, 4)Present1.660885338926581.695938996922561.725882660901951.75191967026518[29]1.660885338926931.695938996921951.725882660901851.75191967026518Lanczos1.660885338951891.695938996922331.725882660901841.75191967026516Figure 4
For Example3, (a) showing the number of iterations spent in the NQM, and (b) the eigenvalues for different n.
(a)(b)Example 4.
As a practical application of the presented method NQM in Section4, we consider Mathieu’s eigenvalue problem [31, 32](52)−y″x+qxyx=λyx,0<x<π,(53)y0=yπ=0,where qx=2e0cos2x with 2e0 a constant amplitude.
To convert equations (52) and (53) into a symmetric matrix eigenvalue problem, let(54)yx=∑k=1nxkϕkx,be the eigenfunction, where xk are coefficients to be determined, and ϕkx are the simplest orthonormal bases(55)ϕkx=2πsinkx.
Multiplying equation (52) by ϕkx, integrating it by parts from x=0 to x=π and using the corresponding boundary conditions, we can derive(56)∫0πy′xϕj′x+qxyxϕjxdx=λ∫0πyxϕjxdx.
Inserting equation (54) for yx into equation (56), taking j=1,…,n, and using the orthonormality of ϕjx, renders an n-dimensional eigenvalue problem (2), where the components aij of the symmetric coefficient matrix A are given by [6](57)aij=1−e0Ifi=j=1,i2Ifi=j≥2,e0Ifi−j=±2.
We employ the NQM in Section4 to find the eigenvalues with n=10, m=9, and x0=110. In Table 3, we list the first ten eigenvalues computed from the present method and compared to that obtained by Liu [6]. By the NQM method, they are convergence within five steps under ε=10−10, and the errors of Ax−λx are between 10−12 to 10−11. Then, we apply the Lanczos algorithm in Section 5 to find the eigenvalues; however, under ε=10−10 the Lanczos algorithm does not converge within 1000 iterations. The main reason is that the characteristic equation generated from the Lanczos algorithm is highly ill-posed; for example, we have p9100=1.04×1014 and p9101=−5.52×1013. The function undergoes a huge jump near the root of the characteristic equation. To compute the first four eigenvalues by the Lanczos algorithm, we must take m=10; otherwise, the eigenvalues obtained are incorrect.
The asymptotic eigenvalues are given by [33].(58)limn⟶∞λn=n2.
In Table4, we list the asymptotic eigenvalues from n=50 to n=300, where the parameters used in the NQM are x0=1n, m=40 for n≥50 and m=80 for n≥100. They are convergence within five steps under ε=10−10. It can be seen that the asymptotic eigenvalues calculated by the present method are coincident to that in equation (58) when n is large.Table 3
For Example4 with e0 = 5, for the first ten eigenvalues comparing the present results (third row) to those obtained by Liu [6] (second row) and the Lanczos algorithm (last row).
λ1λ2λ3λ4λ5λ6λ7λ8λ9λ10−5.79012.09959.236316.648225.510936.358949.269664.204781.7732100.6887−5.790082.099469.2363316.6482225.5109036.3589149.2695864.2047181.77327100.6887−5.790082.099469.2363316.6482222.6388936.2304049.0484564.1877081.77295100.6887Table 4
For Example4 with e0 = 5, the asymptotic eigenvalues for different n.
n50100150200250300λ2500.1310000.0622500.0440000.0362500.0390000.02


## 7. Discussions

One may assert that the shortcoming of the characteristic curveFλ in equation (20) is attributed to the Rayleigh quotient computed in the m-dimensional affine Krylov subspace. Even with m=n for a full space, the shortcoming of Fλ still exists. Now, we equal m to n and solve the eigenvalue problem in the n-dimensional space. For Example 1, we plot the two characteristic curves Fλ and Gλ with respect to λ∈2,10 in Figure 5(a), where we take m=n=3 and x0=1,−1,1T.Figure 5
For Example1, (a) showing characteristic curves obtained in n-dimensional space, (b) the peaks of μ1 and μ2, and (c) the peaks of μ3 and μ4.
(a)(b)(c)Comparing Figures1(a)–5(a), we can see that the characteristic curve Fλ based on the Rayleigh quotient is still exhibiting two crossing points near to λ=3 and λ=6 and one turning point near to λ=9. Moreover, a fictitious root is happened near to 3.488, which is not the eigenvalue. The zigzags of the characteristic curve of Fλ occur because the values of μ1 and μ2 are very small in the order of 10−27 as shown in Figure 5(b), such that the ratio μ2/μ1 in Fλ renders the zigzags.On the other hand, theGλ curve crosses the zero line at λ=3,6,9 as expected, and for the same reason the characteristic curve of Gλ exhibits the smaller zigzags, because the values of μ3 and μ4 are quite small in the order of 10−13 as shown in Figure 5(c), where three peaks of μ3 and μ4 are happened at λ=3,6,9. Besides these three points at λ=3,6,9, the values of μ3 and μ4 are very small owing to the small value of x≈0.In order to reflect these observations, we can prove the following results.Theorem 2.
Under the conditionm=n, we have x=0, IfλλA, x≠0andα+UTx0∈NullB−λIn, Ifλ∈λA, where λA denotes the set of all eigenvalues of A.Proof.
Because ofm=n, U is an n×n matrix, whose inversion is equal to UT in view of equation (6):(59)U−1=UT.
We first prove that(60)α=−U−1x0,is a solution of equation (19) with m=n. From equations (59) and (60), it follows that(61)B−λInα=−B−λInU−1x0=λUTx0−UTAx0,which is just equation (19) with m=n. Thus, equation (60) is true. Next, inserting it into equation (9) yields(62)x=x0+Uα=x0−UU−1x0=0.
On the other hand, we can recast equation (61) to(63)B−λInα=−B−λInUTx0,due to equation (59). Removing the right side to the left side, yields(64)B−λInα+UTx0=0.
Whenλ is an eigenvalue of A, it is also an eigenvalue of UTAU, and meanwhile the null space of UTAU−λIn is not empty. Hence, it leads to(65)α+UTx0≠0.
Then, by equations (9) and (59), we have(66)x=x0+Uα=UUTx0+α≠0.
We complete the proof of this theorem.
Owing to the above theoretical results, we can explain why the values ofμ1 and μ2 in Figure 5(b), and μ3 and μ4 in Figure 5(c) are very small in the whole interval besides at the eigenvalue points, at which the picks appear.
When we apply the pivoting Gaussian elimination method to solve equation (61) whose right-hand side is not zero, the nonzero solution of x is obtained at the eigenvalue point. For this reason, we can explain why μ1 and μ2 in Figure 5(b) and μ3 and μ4 in Figure 5(c) exhibit pick values at the eigenvalue points. When m=n and the matrix A is given, we can check the values of x to detect the approximate positions of eigenvalues. For Example 2, we plot x versus the eigen-parameter λ in Figure 6, where we can observe four peaks happened at the eigenvalues λ=−4,−2,8,and12.Figure 6
For Example2, showing four picks of x, corresponding to eigenvalues −2, −4, 12, and 8.Example 5.
We apply the abovementioned detection method to [8].(67)T=1200234004560067.
We plot the response curve ofx versus the eigen-parameter λ in Figure 7, where we can observe four peaks happened at the eigenvalues λ=−2.4848,0.7046,4.9366,and12.831.
In order to detect the eigenvalue in a more precise location to obtain more accurate eigenvalue, we will develop an iterative detection algorithm in the next section.Figure 7
For Example5, showing four picks of x, corresponding to eigenvalues −2.4848, 0.7046, 4.9366, and 12.831.


## 8. Nonsymmetric Eigenvalue Problems

### 8.1. A New Detection Method of Eigenvalues

WhenA is not a symmetric matrix, the Rayleigh quotient is not applicable to determine the eigenvalue. However, we can derive equation (19) by projecting the vector eigen-equation (2) into the Krylov subspace Km by using the Galerkin condition(68)UTAx−λx=0m.Inserting equation (9) for x and using equation (6) yields(69)UTAx0+Uα−λUTx0+Uα=0m,which by equations (6) and (14) leads to equation (19) again. Obviously, two different approaches render the same equation (19), which is the basic equation for solving the matrix eigenvalue problems.Upon comparing to the original vector eigen-equationA−λInx=0n in equations (2) and (19) is different for the appearance of a nonhomogeneous term λUTx0−UTAx0 in the right side, which is the projection of equation (2) in the affine Krylov subspace, not in the Krylov subspace. To excite the nonzero response of x, we must give a nonzero exciting vector x0≠0.In the proof of Theorem2, we do not need the symmetric property of A. When A is a nonsymmetric matrix having real eigenvalues, α is solved from equation (19) with m=n:(70)B−λInα=UTλx0+Ax0,which is not necessary given by equation (62). Therefore, x is not a zero vector when λ is not an eigenvalue.Later, we give examples to show that the peaks ofxλ at the eigenvalues are happened in the response curve of xλ versus the eigen-parameter λ, which motivates us using a simple maximum method to determine the eigenvalue by collocating points inside an interval by(71)maxλ∈a,bxλ,where the size of a,b must be large enough to include at least one real eigenvalue λ∈a,b. Therefore, the numerical procedures for determining the eigenvalue of a given nonsymmetric matrix are summarized as follows: (i) select m, x0, a, and b. (ii) Construct U. (iii) For each collocating point λi∈a,b solving equation (19), setting x=x0+Uα and taking the optimal λi to meet equation (71).In equation (71), we carry out the computation of the peak point by collocating points at λi=a+ib−a/N0,i=0,…,N0 and pick up the maximal point denoted by λ. We can write a program by gradually decaying the size of the interval centered at the previous peak point by renewing the interval to a finer one. Give an initial interval ak,bk with k=0, and we fix the collocating points by λi=ak+ibk−ak/N0,i=0,…,N0 and pick up the maximal point denoted by λk. Then, a finer interval is given by ak+1=λk−bk−ak/N1 and bk+1=λk+bk−ak/N1, and in that interval, we pick up the new maximum point denoted by λk+1. Continuing this process until λk+1−λk<ε, we can obtain the eigenvalue with high precision. We fix N0=5000 and N1=1000. This algorithm is shortened as an iterative detection method (IDM). The technique of IDM is restricted to detect the real eigenvalues.In order to construct the response curve, we choose a large interval ofa,b to include all eigenvalues, such that the rough locations of the eigenvalues can be observed in the response curve as the peaks. Then, to precisely determine the individual eigenvalue, we choose a small initial interval a0,b0 to include that eigenvalue as internal point. A few iterations by the IDM can compute very accurate eigenvalue.When the eigenvalueλ is obtained, if one wants to compute the eigenvector, we can normalize the nth-component of x by xn=1. Let(72)cij=aij−λδij,dj=−cjn,where aij are the components of A and δij is the Kronecker delta symbol. Then, it follows from equation (2) an (n − 1)-dimensional linear system:(73)cijxj=dj,i,j=1,…,n−1.We can apply the Gaussian elimination method to solvex.


### 8.2. Examples

Example 6.
Now, we can apply the detection method to a nonsymmetric matrix [8].(74)A=−261209−49−530422−98−800631−144.
We plotx versus the eigen-parameter λ in Figure 8, where we can observe three peaks happened at the eigenvalues λ=3,4,and10.
For this example, withm=3, a=2.5, and b=3.5, we can obtain the eigenvalue λ=3 exactly by using equation (71), and with a,b=9,11, we can obtain the eigenvalue λ=10. They are exact eigenvalues.Figure 8
For Example6 of a nonsymmetric matrix, showing three picks of x, corresponding to eigenvalues 3, 4, and 10.Example 7.
We consider another nonsymmetric matrix:(75)A=3.548815.5938.5775−4.01232.359524.52614.596−5.81570.08995327.59921.483−5.84151.922755.66739.717−10.558,where λA=1,2,6,30. With m=4, we plot x versus the eigen-parameter λ in Figure 9(a), where we can observe four peaks happened at the eigenvalues λ=1,2,6,30, while with m=3, we plot x versus the eigen-parameter λ in Figure 9(b), where we can observe three peaks happened at the eigenvalues λ=2,6,and30.
For this example, withm=4, a=1.8, and b=2.3, we can obtain the eigenvalue λ=2 by using the IDM in equation (71), and with a=5.5 and b=6.3, we can obtain the eigenvalue λ=6. They are exact eigenvalues. But with m=3, a=5.5, and b=6.3, we can obtain the eigenvalue λ=6.0208, whose accuracy is reduced.Figure 9
For Example7 of a nonsymmetric matrix, (a) with m=4 showing four picks of x, corresponding to eigenvalues 1, 2, 6, and 30, and (b) with m=3 showing three picks of x, corresponding to eigenvalues 2, 6, and 30.
(a)(b)Example 8.
Consider the Frank nonsymmetric matrix [15, 34, 35](76)A=111…1122…2023…3⋮0⋱…⋮00…n−1n.
Forn=30, the maximum eigenvalue λ1 is given by(77)λ1=96.200622293285.
We takem=20 and plot xλ versus the eigen-parameter λ in Figure 10 by a solid line, where we can observe three peaks happened at the first three eigenvalues. We find that m can be reduced to m=8, whose response curve still has three picks near to the one with m=20 as shown in Figure 10 by a dashed line.
Starting froma0=95 and b0=97 and through six iterations under a convergence criterion 10−15, we can obtain λ1=96.200622293285620, which is very close to the one in equation (77). It is remarkable that by using equation (73), we can obtain Ax−λx=2.99×10−11. Starting from a0=77 and b0=78 and through six iterations, we can obtain λ2=77.34401431966266 and Ax−λx=9.42×10−11. Starting from a0=63 and b0=64 and through six iterations, we can obtain λ2=63.21561592564849 and Ax−λx=2.42×10−10.Figure 10
For Example8 of the Frank nonsymmetric matrix with n=30, showing three picks of x in the interval [60, 99], m=20 solid line, and m=8 dashed line.Example 9.
The abovementioned techniques are applicable to the eigenvalue problems of symmetric matrices. We apply the new algorithm IDM to solve equation (2) with [8].(78)A=111111123456136101521141020355615153570126162156126252.
We takem=3 and plot xλ versus λ in Figure 11(a) in an interval [0, 16], where we can observe two peaks happened at the first two eigenvalues. For m=4, a new peak for the third eigenvalue appears but the peak for the first eigenvalue seems disappeared. However, when m=4, the pick for the first eigenvalue is quite small as shown in Figure 11(b) in an interval [15, 16]. The values of xλ are greatly reduced when m is changed from m=3 to m=4.
By takingx0=16/6 and under ε=10−10, the RQIM through seven iterations to find λ=15.5534732738 with Ax−λx=3.85×10−8.
By takingx0=16 and m=5 in the IDM, starting from a0=15 and b0=16 and through seven iterations under a convergence criterion 10−15, we can obtain λ1=15.55347327375157, which is very close to λ=15.5534732738 obtained by the RQIM. By using equation (73), we can obtain the corresponding eigenvector x whose error Ax−λx to satisfy the eigen-equation (2) is 2.74×10−14, which is more accurate than 3.85×10−8 obtained by the RQIM.
Starting froma0=2 and b0=2, with m=5 and through seven iterations, we can obtain λ2=2.043573780089090 and Ax−λx=7.64×10−14. Starting from a0=0.4 and b0=1, with m=5 and through six iterations, we can obtain λ3=0.489338828743265 and Ax−λx=1.7×10−11.
Starting froma0=0 and b0=0.004, with m=6 and through six iterations, we can obtain a small eigenvalue λ4=0.003004389574740893 and Ax−λx=1.23×10−13.Figure 11
For Example9, (a) in the interval 0,16 showing two picks of x with m=3 and with m=4, and (b) a pick of x with m=4 is enlarged in the interval 15,16.
(a)(b)Example 10.
We revisit Example2 again. We take m=4 and plot xλ versus the eigen-parameter λ in Figure 12 in an interval −5,13, where we can observe four peaks happened at four eigenvalues. In Table 5, the IDM is more accurate than the Lanczos algorithm.Figure 12
For Example10, showing four picks of x, corresponding to eigenvalues −4, −2, 8, and 12.Table 5
For example 10 solved by the IDM and the Lanczos algorithm, comparing the errors of eigenvalues.
λ−4−2812IDM06.43×10−54.11×10−55.01×10−5Lanczos7.31×10−36.43×10−59.26×10−45.80×10−3Example 11.
We revisit the eigenvalue problem for the Hilbert matrix; however, we solve this eigenvalue problem by using the IDM in Section8.1. We take n=m=5 and plot xλ versus the eigen-parameter in Figure 13 in an interval [0, 2], where we can observe five peaks happened at five eigenvalues.
Due to highly ill-conditioned nature of the Hilbert matrix, it is a quite difficult eigenvalue problem. For this problem, we taken=100, m=30, and j0 = 50 to compute the largest eigenvalue, which is given as 2.182696097757424. The affine Krylov subspace method is convergent very fast with six iterations, and the error of eigen-equation is Ax−λx=2.005×10−13. The CPU time for the affine Krylov subspace method together with the IDM is 1.52 s.
However, finding the largest eigenvalue in the full space with dimensionn=100, it does not converge within 100 iterations. By using the MATLAB, we obtain 2.182696097757423 which is close to that obtained by the IDM. The QR method implemented in MATLAB computes all the eigenvalues in the diagonal matrix D and the eigenvectors V of A in 0.000812 s with an accuracy equal to AV−VD=1.1714×10−15, which in comparison to the IDM is much fast and more accurate than the IDM for the largest eigenvalue.
Notice that the smallest eigenvalue is very difficult to be computed, since it is very close to zero. However, we can obtain the smallest eigenvalue3.1×10−19 with one iteration and the error Ax−λx=9.412×10−17 is obtained. For the eigenvalue problem of the Hilbert matrix, the MATLAB leads to a wrong eigenvalue −9.693294591063901×10−17, which is negative and contradicts to the positive eigenvalues of the Hilbert matrix. The eigenfunction in MATLAB cannot guarantee to obtain a positive eigenvalue for the positive definite Hilbert matrix. The first 41 eigenvalues are all negative. The first 73 eigenvalues are less than 10−17. So, most of these eigenvalues computed by MATLAB should be spurious. The MATLAB is effective for general purpose eigenvalue problem with normal matrices, but for the highly ill-conditioned matrices, the effectiveness of the MATLAB might be lost. Due to the highly ill-posed nature of the smallest eigenvalue of the Hilbert matrix, we may consider the QR algorithm with shift proposed by Francis [36, 37] to resolve this difficulty.
Then, we apply the Lanczos algorithm together with the Newton method to determine the largest eigenvalue of the Hilbert matrix withn=100, where we take m=30. The Lanczos algorithm spent 26 iterations to satisfy the convergence criterion ε=10−15. The largest eigenvalue obtained is 2.182664504835603 with a poor accuracy of Ax−λx=2.313×10−3. For the IDM, six iterations and Ax−λx=2.005×10−13 are obtained, which is much more accurate and efficient than the Lanczos algorithm.Figure 13
For Example11, showing five picks of x, corresponding to the eigenvalues of Hilbert matrix with n = 5.


## 8.1. A New Detection Method of Eigenvalues

WhenA is not a symmetric matrix, the Rayleigh quotient is not applicable to determine the eigenvalue. However, we can derive equation (19) by projecting the vector eigen-equation (2) into the Krylov subspace Km by using the Galerkin condition(68)UTAx−λx=0m.Inserting equation (9) for x and using equation (6) yields(69)UTAx0+Uα−λUTx0+Uα=0m,which by equations (6) and (14) leads to equation (19) again. Obviously, two different approaches render the same equation (19), which is the basic equation for solving the matrix eigenvalue problems.Upon comparing to the original vector eigen-equationA−λInx=0n in equations (2) and (19) is different for the appearance of a nonhomogeneous term λUTx0−UTAx0 in the right side, which is the projection of equation (2) in the affine Krylov subspace, not in the Krylov subspace. To excite the nonzero response of x, we must give a nonzero exciting vector x0≠0.In the proof of Theorem2, we do not need the symmetric property of A. When A is a nonsymmetric matrix having real eigenvalues, α is solved from equation (19) with m=n:(70)B−λInα=UTλx0+Ax0,which is not necessary given by equation (62). Therefore, x is not a zero vector when λ is not an eigenvalue.Later, we give examples to show that the peaks ofxλ at the eigenvalues are happened in the response curve of xλ versus the eigen-parameter λ, which motivates us using a simple maximum method to determine the eigenvalue by collocating points inside an interval by(71)maxλ∈a,bxλ,where the size of a,b must be large enough to include at least one real eigenvalue λ∈a,b. Therefore, the numerical procedures for determining the eigenvalue of a given nonsymmetric matrix are summarized as follows: (i) select m, x0, a, and b. (ii) Construct U. (iii) For each collocating point λi∈a,b solving equation (19), setting x=x0+Uα and taking the optimal λi to meet equation (71).In equation (71), we carry out the computation of the peak point by collocating points at λi=a+ib−a/N0,i=0,…,N0 and pick up the maximal point denoted by λ. We can write a program by gradually decaying the size of the interval centered at the previous peak point by renewing the interval to a finer one. Give an initial interval ak,bk with k=0, and we fix the collocating points by λi=ak+ibk−ak/N0,i=0,…,N0 and pick up the maximal point denoted by λk. Then, a finer interval is given by ak+1=λk−bk−ak/N1 and bk+1=λk+bk−ak/N1, and in that interval, we pick up the new maximum point denoted by λk+1. Continuing this process until λk+1−λk<ε, we can obtain the eigenvalue with high precision. We fix N0=5000 and N1=1000. This algorithm is shortened as an iterative detection method (IDM). The technique of IDM is restricted to detect the real eigenvalues.In order to construct the response curve, we choose a large interval ofa,b to include all eigenvalues, such that the rough locations of the eigenvalues can be observed in the response curve as the peaks. Then, to precisely determine the individual eigenvalue, we choose a small initial interval a0,b0 to include that eigenvalue as internal point. A few iterations by the IDM can compute very accurate eigenvalue.When the eigenvalueλ is obtained, if one wants to compute the eigenvector, we can normalize the nth-component of x by xn=1. Let(72)cij=aij−λδij,dj=−cjn,where aij are the components of A and δij is the Kronecker delta symbol. Then, it follows from equation (2) an (n − 1)-dimensional linear system:(73)cijxj=dj,i,j=1,…,n−1.We can apply the Gaussian elimination method to solvex.


## 8.2. Examples

Example 6.
Now, we can apply the detection method to a nonsymmetric matrix [8].(74)A=−261209−49−530422−98−800631−144.
We plotx versus the eigen-parameter λ in Figure 8, where we can observe three peaks happened at the eigenvalues λ=3,4,and10.
For this example, withm=3, a=2.5, and b=3.5, we can obtain the eigenvalue λ=3 exactly by using equation (71), and with a,b=9,11, we can obtain the eigenvalue λ=10. They are exact eigenvalues.Figure 8
For Example6 of a nonsymmetric matrix, showing three picks of x, corresponding to eigenvalues 3, 4, and 10.Example 7.
We consider another nonsymmetric matrix:(75)A=3.548815.5938.5775−4.01232.359524.52614.596−5.81570.08995327.59921.483−5.84151.922755.66739.717−10.558,where λA=1,2,6,30. With m=4, we plot x versus the eigen-parameter λ in Figure 9(a), where we can observe four peaks happened at the eigenvalues λ=1,2,6,30, while with m=3, we plot x versus the eigen-parameter λ in Figure 9(b), where we can observe three peaks happened at the eigenvalues λ=2,6,and30.
For this example, withm=4, a=1.8, and b=2.3, we can obtain the eigenvalue λ=2 by using the IDM in equation (71), and with a=5.5 and b=6.3, we can obtain the eigenvalue λ=6. They are exact eigenvalues. But with m=3, a=5.5, and b=6.3, we can obtain the eigenvalue λ=6.0208, whose accuracy is reduced.Figure 9
For Example7 of a nonsymmetric matrix, (a) with m=4 showing four picks of x, corresponding to eigenvalues 1, 2, 6, and 30, and (b) with m=3 showing three picks of x, corresponding to eigenvalues 2, 6, and 30.
(a)(b)Example 8.
Consider the Frank nonsymmetric matrix [15, 34, 35](76)A=111…1122…2023…3⋮0⋱…⋮00…n−1n.
Forn=30, the maximum eigenvalue λ1 is given by(77)λ1=96.200622293285.
We takem=20 and plot xλ versus the eigen-parameter λ in Figure 10 by a solid line, where we can observe three peaks happened at the first three eigenvalues. We find that m can be reduced to m=8, whose response curve still has three picks near to the one with m=20 as shown in Figure 10 by a dashed line.
Starting froma0=95 and b0=97 and through six iterations under a convergence criterion 10−15, we can obtain λ1=96.200622293285620, which is very close to the one in equation (77). It is remarkable that by using equation (73), we can obtain Ax−λx=2.99×10−11. Starting from a0=77 and b0=78 and through six iterations, we can obtain λ2=77.34401431966266 and Ax−λx=9.42×10−11. Starting from a0=63 and b0=64 and through six iterations, we can obtain λ2=63.21561592564849 and Ax−λx=2.42×10−10.Figure 10
For Example8 of the Frank nonsymmetric matrix with n=30, showing three picks of x in the interval [60, 99], m=20 solid line, and m=8 dashed line.Example 9.
The abovementioned techniques are applicable to the eigenvalue problems of symmetric matrices. We apply the new algorithm IDM to solve equation (2) with [8].(78)A=111111123456136101521141020355615153570126162156126252.
We takem=3 and plot xλ versus λ in Figure 11(a) in an interval [0, 16], where we can observe two peaks happened at the first two eigenvalues. For m=4, a new peak for the third eigenvalue appears but the peak for the first eigenvalue seems disappeared. However, when m=4, the pick for the first eigenvalue is quite small as shown in Figure 11(b) in an interval [15, 16]. The values of xλ are greatly reduced when m is changed from m=3 to m=4.
By takingx0=16/6 and under ε=10−10, the RQIM through seven iterations to find λ=15.5534732738 with Ax−λx=3.85×10−8.
By takingx0=16 and m=5 in the IDM, starting from a0=15 and b0=16 and through seven iterations under a convergence criterion 10−15, we can obtain λ1=15.55347327375157, which is very close to λ=15.5534732738 obtained by the RQIM. By using equation (73), we can obtain the corresponding eigenvector x whose error Ax−λx to satisfy the eigen-equation (2) is 2.74×10−14, which is more accurate than 3.85×10−8 obtained by the RQIM.
Starting froma0=2 and b0=2, with m=5 and through seven iterations, we can obtain λ2=2.043573780089090 and Ax−λx=7.64×10−14. Starting from a0=0.4 and b0=1, with m=5 and through six iterations, we can obtain λ3=0.489338828743265 and Ax−λx=1.7×10−11.
Starting froma0=0 and b0=0.004, with m=6 and through six iterations, we can obtain a small eigenvalue λ4=0.003004389574740893 and Ax−λx=1.23×10−13.Figure 11
For Example9, (a) in the interval 0,16 showing two picks of x with m=3 and with m=4, and (b) a pick of x with m=4 is enlarged in the interval 15,16.
(a)(b)Example 10.
We revisit Example2 again. We take m=4 and plot xλ versus the eigen-parameter λ in Figure 12 in an interval −5,13, where we can observe four peaks happened at four eigenvalues. In Table 5, the IDM is more accurate than the Lanczos algorithm.Figure 12
For Example10, showing four picks of x, corresponding to eigenvalues −4, −2, 8, and 12.Table 5
For example 10 solved by the IDM and the Lanczos algorithm, comparing the errors of eigenvalues.
λ−4−2812IDM06.43×10−54.11×10−55.01×10−5Lanczos7.31×10−36.43×10−59.26×10−45.80×10−3Example 11.
We revisit the eigenvalue problem for the Hilbert matrix; however, we solve this eigenvalue problem by using the IDM in Section8.1. We take n=m=5 and plot xλ versus the eigen-parameter in Figure 13 in an interval [0, 2], where we can observe five peaks happened at five eigenvalues.
Due to highly ill-conditioned nature of the Hilbert matrix, it is a quite difficult eigenvalue problem. For this problem, we taken=100, m=30, and j0 = 50 to compute the largest eigenvalue, which is given as 2.182696097757424. The affine Krylov subspace method is convergent very fast with six iterations, and the error of eigen-equation is Ax−λx=2.005×10−13. The CPU time for the affine Krylov subspace method together with the IDM is 1.52 s.
However, finding the largest eigenvalue in the full space with dimensionn=100, it does not converge within 100 iterations. By using the MATLAB, we obtain 2.182696097757423 which is close to that obtained by the IDM. The QR method implemented in MATLAB computes all the eigenvalues in the diagonal matrix D and the eigenvectors V of A in 0.000812 s with an accuracy equal to AV−VD=1.1714×10−15, which in comparison to the IDM is much fast and more accurate than the IDM for the largest eigenvalue.
Notice that the smallest eigenvalue is very difficult to be computed, since it is very close to zero. However, we can obtain the smallest eigenvalue3.1×10−19 with one iteration and the error Ax−λx=9.412×10−17 is obtained. For the eigenvalue problem of the Hilbert matrix, the MATLAB leads to a wrong eigenvalue −9.693294591063901×10−17, which is negative and contradicts to the positive eigenvalues of the Hilbert matrix. The eigenfunction in MATLAB cannot guarantee to obtain a positive eigenvalue for the positive definite Hilbert matrix. The first 41 eigenvalues are all negative. The first 73 eigenvalues are less than 10−17. So, most of these eigenvalues computed by MATLAB should be spurious. The MATLAB is effective for general purpose eigenvalue problem with normal matrices, but for the highly ill-conditioned matrices, the effectiveness of the MATLAB might be lost. Due to the highly ill-posed nature of the smallest eigenvalue of the Hilbert matrix, we may consider the QR algorithm with shift proposed by Francis [36, 37] to resolve this difficulty.
Then, we apply the Lanczos algorithm together with the Newton method to determine the largest eigenvalue of the Hilbert matrix withn=100, where we take m=30. The Lanczos algorithm spent 26 iterations to satisfy the convergence criterion ε=10−15. The largest eigenvalue obtained is 2.182664504835603 with a poor accuracy of Ax−λx=2.313×10−3. For the IDM, six iterations and Ax−λx=2.005×10−13 are obtained, which is much more accurate and efficient than the Lanczos algorithm.Figure 13
For Example11, showing five picks of x, corresponding to the eigenvalues of Hilbert matrix with n = 5.


## 9. Conclusions

In conclusion, we can claim that the weakness of the Rayleigh quotient used to generate the first characteristic equation is inherent, which is happened no matter whichm≤n is taken. Upon using the new quotient in the second characteristic equation, we can improve this property which can help us to quickly find the eigenvalue and the eigenvector. There already existed many excellent methods and algorithms to solve the symmetric eigenvalue problems, and the current paper proposed an alternative approach, which merely needs an initial guess of λ0. Therefore, we can guess different initial values λ0 to quickly find all eigenvalues. In many algorithms such as the Rayleigh quotient iteration method, they required the initial guess of the eigenvector, which is more difficult than that giving the initial guess of λ0 to achieve certain eigenvalue. The presented algorithms RQM and NQM can find the eigen-pairλ,x simultaneously. The computational cost is saving for solving merely a lower m-dimensional linear system and the derived characteristic equation, whose high efficiency and high accuracy were identified by the numerical tests. Through some tests, we found that the proposed new quotient method can find more eigenvalues than that obtained by the Lanczos algorithm, if the same value of m is used. If m is large enough and m<n, the presented method can find all eigenvalues very effectively. For the eigenvalue problems of nonsymmetric matrices, we have developed a simple and powerful technique to determine the eigenvalues as the peaks of the response curve of xλ with respect to the eigen-parameter λ in a selected interval. The vector x was expressed in terms of a nonzero exciting vector and an m-vector developed in the m-dimensional affine Krylov subspace, whose governing linear equation was derived by using the Galerkin condition. An iterative detection method was developed, where through a few finer tuning in a smaller interval, very accurate eigenvalue can be obtained. Also then, we reduced the eigen-equation with one dimension to a nonhomogeneous linear system to determine the eigenvector with high accuracy.



---

*Source: 9859889-2023-03-06.xml*