The QR Method for Determining All Eigenvalues of Real Square Matrices

Eigenvalues are special sets of scalars associated with a given matrix. In other words for a given matrix A, if there exist a non-zero vector V such that, AV= λV for some scalar λ, then λ is called the eigenvalue of matrix A with corresponding eigenvector V. The set of all nxm matrices over a field F is denoted by Mnm (F). If m = n, then the matrices are square, and which is denoted by Mn (F). We omit the field F = C and in this case we simply write Mnm or Mn as appropriate. Each square matrix AεMnm has a value in R associated with it and it is called its determinant which is use full for solving a system of linear equation and it is denoted by det (A). Consider a square matrix AεMn with eigenvalues λ, and then by definition the eigenvectors of A satisfy the equation, AV = λV, where v={v1, v2, v3............vn}. That is, AV=λV is equivalent to the homogeneous system of linear equation (A-λI) v=0. This homogeneous system can be written compactly as (A-λI) V = 0 and from Cramer’s rule, we know that a linear system of equation has a non-trivial solution if and only if its determinant is zero, so the solution λ is given by det (A-λI) =0. This is called the characteristic equation of matrix A and the left hand side of the characteristic equation is the characteristic polynomial whose roots are equals to λ.


Introduction
Each nxn square matrix A has a value associated with its determinant denoted by det (A), which is useful to solve a system of linear equation and moreover, if there exist a nonzero vector V and scalar λ such that AV = λ V, then λ is called eigenvalue of matrix A with corresponding eigenvector V. It is bare fact that, AV = λ V is equivalent to = where I is the identity matrix then, by Cramer's rule we know that, linear system of equation has non-trivial solution if det 0 A ≠ , so that λ can be computed from this equation.
Applications of eigenvalue problems play a great role in our real word. One class of applications which has recently gained considerable ground is that related to eigenvalues problems of a matrix. However, the most commonly solved eigenvalue problems today are those issues associated with the vibration analysis of large structures. The vibrational frequencies are therefore determined by the eigenvalues of a symmetric 3x3 matrix. This is an instance of simple eigenvalue problem AV V λ = that is common in practice.
However, Since computing the eigenvalue λ is equivalent to finding the roots of matrix's characteristic polynomial, we can see that task is quickly too difficult for larger dimensional matrices/especially for the case of matrices which have more than three dimensions/ even if we know characteristic polynomial and algorithms such as Newton's method for finding zeros cannot be depended upon produce all the zeros with reasonable speed and accuracy. Fortunately, numerical analysts have found an entirely different ways to calculate eigenvalues of a given square matrix. Among those methods QR method is the most widely used, important, accurate and speedy one. The primary reason that modern implementations of this method are efficient and reliable is that a QR factorization can be used to create each new matrix in the sequence and each QR factorization can be calculated quickly and accurately; it yields easily a new matrix orthogonally similar to the original matrix; and orthogonal similarities tend to minimize the effect of round off error on the eigenvalues.

Householder Transformation
There is no simple way to calculate eigenvalues for a matrix larger than 2x2 dimensions. The method of calculating the characteristic polynomial and then finding its zeros is not good numerically and moreover finding the roots of characteristic polynomial involves taking a determinant which uses large amount of computing time. The QR method is one of the most important methods which used to find eigenvalues of real square matrix. Therefore, the main idea of this chapter is determining all eigenvalues of real square matrix by using QR factorization (where Q is orthogonal and R is upper triangular matrices). To this end, suppose a real square A is given. Let A= Q 0 R 0 be QR factorization of A (where Q is orthogonal and R is upper triangular matrices), and create A 1 = R 0 Q 0 . Let A 1 = Q 1 R 1 be QR factorization of A 1 and similarly create A 2 = R 1 Q 1 , continue this process in the same fashion for ≥ 1, ( ℎ = ). Once A m has been created such that, A m = Q m R m , and A m+1 = R m Q m. Thus, the sequence {A m } will usually converges to something from which the eigenvalues can be computed easily. Moreover, A 2 is similar to A and A m+1 is similar to A m and so on. Hence A m and A m+1 have the same eigenvalues as matrix A by definition of similarity of matrices. This process will be stopped when the entries below the main diagonal of current matrix A m are sufficiently small, or if it appears that convergence will not happen. This implies that, QR factorization sequence process can fail to converge or the convergence can be also extremely slow and expensive. If it converges to certain matrix, then the diagonal entries of this current matrix are tends to be eigenvalues of matrix A, if not, it can be modified in order to speed up convergence or to accelerate the rate of convergence of the given real square matrix dramatically we use methods like shifting of origin.
The method of orthogonal triangularization (orthogonal transformation) is analogous to Gaussian elimination. Since orthogonal transformation will not worsen the condition or stability of eigenvalues of a non-symmetric matrix we will attempt to decompose an arbitrary real square matrix A into a product QR, where Q is orthogonal and R is upper triangular matrices. Now, let us begin our discussion by looking at a special class of orthogonal matrices known as Householder matrices. Theorem 1: Let U be Householder matrix (U = I-2ww*). Then i). U is Hermitian (U*=U) ii). U is orthogonal (U*U=I). The Householder's matrices will be used to transform a nonzero vector in to a new vector containing mainly zeros. Let b 0 ≠ be given, b ∈ ℜ and suppose we want to produce U of form (1) such that Ub contains zeros in position r+1 through n, for some given 1 r ≥ . Choose w as in (2). Then the first r-1 columns of b and Ub are the same. To simplify the later work, let m = n-r+1 then Then our restriction on the form Ub requires the first r-1 components of Ub to be c, and (I − 2ww 5 )6 = (7, 0, … ,0) % , = 1 for some 7 (3) Since I − 2ww 5 is orthogonal the length of d is preserved, and thus |7| = 6 = : 7 = ±: = ±<6 + 6 + ⋯ + 6 > (4) Define: From (3) we have; (I − 2ww 5 )6 = 6 − 2ww 5 6 = (7, 0, … ,0) % (6) Multiply by % and using = 1 implies that, P =−7 Substituting this in to the first component of (5) give us; 6 + 27 = 7 which imply that, Choose the sign of 7 in (4) by Sign (7) = -sign (6 ). This choice maximizes J , and it avoids any possible loss of significance error in the calculation of . The sign for is irrelevant. Having it is possible to obtain p from (6). Return to equation (5), and using components 2 through m, The equation (4) and (6) to (8) completely define w and thus w and U. The operation count is 2m + 2 multiplication and division, and one square root. A sequence of such operation is used to systematically reduce the given matrix to simpler form. Now, given a real square matrix A, we will show there is an orthogonal matrix Q and an upper triangular matrix R for which, Let ? P = Q − 2 (P) (P)% , where r = 1, 2, 3…, With (P) as is defined in equation (2) with r-1 leading zeros. Writing A interims of its column * , * , … , * , then we have, ? = (? * , ? * , … , ? * ) % . Pick ? and using the preceding construction (4) to (8) with b= * . Then, ? contains zeros below the diagonal in its first column. Choose ? similarly, so that R R will contain zeros in its second column below the diagonal. Note that because contains a zero in position one, and R is zero in the first column below position one, thus the product R R and R contains the same element in row and column one. Now choose ? and as before in (4) to (8) with b equal to the second column of ? by carrying out with each column of A, we obtain an upper triangular matrix, If at step r of the construction, all elements below the diagonal of column r are zero, then just choose ? P = Q and go on to the next step. To complete the construction, define, This is orthogonal. Then A = QR as desired. It is useful to know to what extent the factorization A = QR is unique?. For non-singular matrix A, suppose Then, S and S must also be nonsingular, and U % U = S S T The inverse of an upper triangular matrix is upper triangular, and the product of two upper triangular matrices is upper triangular. Thus S S T is upper triangular. Also the product of two orthogonal matrices is orthogonal; and thus the product U % U is orthogonal. But it is not hard to show that the only upper triangular orthogonal matrices are the diagonal matrices for some diagonal matrix D.
Since we are only dealing real matrices, D has a diagonal element equals to +1 or -1 combining the results we have, Another practical matter is deciding how to evaluate the matrix R of (11). Let, P = R P PT = (Q − 2 (P) (P)% ) PT , where r = 1, 2, 3, … , n (15) With X = , T = S. If we calculate R P and then multiply it by T to form P , the number of multiplication will be; There is a much more efficient method for calculating P . By rewriting (15) as; First calculate
Multiplications, which shows (16) is a preferable way to evaluate each P and finally S = T .

Reduction of a Matrix to Tridiagonal and Hessenberg Matrix Form
Let A be a real square matrix. To find the eigenvalues of A, it is usually first reduced to tridiagonal /if the given matrix is symmetric/or upper Hessenberg matrix form by orthogonal similarity transformations (since similar matrices have the same eigenvalues). Then, the eigenvalues of the reduced matrix can be calculated by using the QR algorithm more efficiently.
Suppose A is a real square matrix. For the orthogonal matrices, we will use the Householder matrices of the above. Let The matrix = ? % ? = ? ? is similar to . The element \ is unchanged, and will be symmetric. Produce ( ) and ? to obtain, ? * = (\ , \ ] , 0, … ,0) % for some \ ] . The vector * is the first column of A. use (4) to (8)  Since is symmetric, the construction on the first column of A will imply that has zeros in position 3 through n of both the first row and column. Continue this process, and letting, With = A. pick R P to introduce zeros in to position r+2 through n of column r. column one through r-1 will remain undisturbed in calculating R P PT , due to the special form of R P .
Pick the vector (P[ ) in analogy with the above description for ( ) . The final matrix a ≡ PT is tridiagonal and symmetric.
If A is non-symmetric, it is reduced to a similar Hessenbeg matrix form using the same algorithm. Reduction to symmetric tridiagonal form requires $ Y $ operations, where as the transformation to Hessenberg matrix form in the nonsymmetric case requires l $ Y $ operations.
This will be much more convenient form for the calculation of the eigenvalues of A; and the eigenvectors of A can easily be obtained from those of T.
Matrix A is related to T b; a = U % U, Where As before with the QR factorization, we seldom produce Q explicitly, preferring to work with the individual matrices R P in analogy with (16). For an eigenvector x of A, say Ax, we have, If we produce an orthogonal set of eigenvectors om p q for T, then oUm p q will be an orthogonal set of eigenvectors for A, since Q preserve length and angles.
Example. 2: Use the Householder's method to reduce the following given matrix A into the tridiagonal form. For an error analysis of this reduction to tridiagonal form, we give some result below. Let the computer arithmetic be binary floating-point with rounding, with t binary digits in the mantissa. Furthermore, assume that all inner products ∑ \ p ƒ p > p& that occur in the calculations are accumulated in the double precision, with rounding to single precision at this completion of the summation. These inner products occur in the variety of places in the computation of T from A. Let a " denote the actual symmetric tridiagonal matrix that is computed from A using the above computer arithmetic. Let R " P denote the actual matrix produced in converting PT … P , let R P be theoretical exact version of this matrix if no rounding errors occurred, and let U = ? . ? … ? PT be the exact product of these R P ′:, an orthogonal matrix.
Let A be a real symmetric matrix of order n. Let a " be the real symmetric tridiagonal matrices resulting from Householder similarity transformation method (17) to A, as in (18).
Assume the floating point arithmetic used has the characteristic described in the preceding paragraph. Let For small and moderate values of n, • ≡ 25(Y − 1), the result above shows that the reduction to tridiagonal form is extremely stable operation, with little new errors introduced for the eigenvalues.

The QR Algorithm
Nowadays, QR method is the most efficient and widely used general method for the calculation of all eigenvalues of a matrix. The QR method is quite complex in both its theory and application, and hence this paper able to give only an introduction to the theory and method.
Given real square matrix A, we know from our previous discussion that, there is a factorization, A = QR With R upper triangular and Q orthogonal matrices, with A real, Q and R can be chosen real; and their construction is as we have discussed before. We assume that A is real throughout this discussion. The basic QR algorithm is given as follows.
QR algorithm: Suppose a real square matrix A is given. Let A= Q 0 R 0 be QR factorization of A, and create A 1 = R 0 Q 0 . Let A 1 = Q 1 R 1 be QR factorization of A and create A 2 = R 1 Q 1. , Continue this process again. For ≥ 1, ( ℎ = ). Once A m has been created such that, A m = Q m R m , and A m+1 = R m Q m (22)  The diagonal entries of the current matrix are approximations of eigenvalues of the given matrix.
The elements in the above matrices row-1, column-2 position decrease geometrically with ratio of about 0.64 per iterate, and those in the row-2, column-3 position decrease with a ratio of about 0.42 per iterate. The value in the row-3, column-3 position of l will be 1.2679, which is correct to four decimal places. The QR algorithm (22) can be relatively expensive because the QR factorization is time consuming when repeated many times. To decrease the expense the matrix is prepared for the QR method by reducing it a simpler form, one for which the QR factorization is much less expensive. If A is symmetric, it is reduced to a similar tridiagonal matrix. If A is non-symmetric, it is reduced to a similar Hessenberg matrix exactly as described in discussion (2). It is upper triangular except for a single non-zero sub diagonal. The matrix A is reduced to Hessenberg matrix form by using the same algorithm as we have done for reducing symmetric matrices to tridiagonal form. With A diagonal or Hessenberg, the Householder matrices take a simple form when calculating the QR factorization. Having produced A 1 = Q 1 R 1 and A 2 = R 1 Q 1, we have to know that the form of A is the same as that of A 1 in order to continue using the less expensive form of QR factorization.
Suppose A 1 is Hessenberg form, then from discussion (1) the factorization has following value for Q 1 : With each • ž a Householder matrix (2. 13) of discussion (2.1), • ž Q " 2 ž ž % , 1 Ž Ÿ Ž Y " 1 (25) Because the matrix A 1 is of the Heisenberg form, the vector ž can be shown to have the special form; This can be shown from the equation for the components of ž , and in particular (18) of equation (26). From equation (25), the matrix • ž will differ from the identity in only the four elements in position (k, k), (k, k+1), (k+1, k) and (k+1, k+1). And from this it is a fairly straight forward computation to show that U must be Hessenberg in form.
Another necessary lemma is that the product of a triangular matrix and a Hessenberg matrix is again Hessenberg. Just multiply the two forms of matrices, observing the perspective patterns of zeros, in order to prove this lemma. Combining these results, observing that S is upper triangular, we have that S U must be in Heisenberg form. If A 1 is symmetric and tri-diagonal, then it is trivially Heisenberg. From the preceding result must also be Hessenberg. But is symmetric since, % U % U % U % % U U % U A and since any symmetric Heisenberg matrix is tridiagonal, we have shown that A tridiagonal. Theorem 3: Let A be a real matrix of order n, and let its eigenvalues satisfying that

The QR Algorithm with Shift/Improvement of Convergence
The basic QR algorithm is too slow for it to be competitive, especially if some of the eigenvalues are close together, but two refinements will accelerate the convergence dramatically. First, we perform a preliminary reduction of the matrix into tridiagonal or upper Heisenberg form, which decreases the computations cost of the future QR iterations. Second, we use multiple translational shift of origin of the matrix to reduce the total number of iterations before convergence. In other word, the QR algorithm is generally applied with a shift of origin for the eigenvalues in order to increase the speed of convergence.
During QR iterations on the matrix m A , where m is the number of the iteration current step, the diagonal entry 7 (>) ∈ m A will approach the eigenvalues n λ , assuming we meet the convergence criteria. Using our current methods, the rate of convergence will be linear with the ratio Similarly, Since the off-diagonal elements of t are approaching to zero or small number, the eigenvalues of A are approaching to diagonal elements of t . Now, let us take the diagonal entries of iteration t as approximated eigenvalues in both cases (which is computed with and without shift/improvement or acceleration) and let us compare those with exact eigenvalues of the same given square matrix A. As it is shown in the table (1), the error with shift is smaller as compare to without shift. It is possible to improve the eigenvalues more than this by applying some more shifting methods and hence the round-off error will decrease with shift as compare to without shift. This implies that, shifting technique is very important in QR method to reduce the cost per iteration but QR method barley by itself is too slow, too expensive and because of this reason it is impractical without improvement technique.

The Double QR Algorithm
Now, let us consider the important special case of real square matrix with complex eigenvalues. In order to realize the faster rate of convergence of general QR method (using shifts) on such matrices complex shifts must be employed at some stage. After such stage the matrices will contain complex elements and so will all the extra work involved in treating complex matrices (double storage, at least four times the work). There is thus storage incentive for seeking an algorithm which produces a real matrix from which the complex eigenvalues can be calculated easily as conjugate root of 2x2 real principal sub matrices.
Consider two steps of the generalized QR transformation of general real matrix .
We have, Thus, setting U X = U U and S X = S S , the formulas above says that $ = U X % U X where, by uniqueness, U X S X is the QR factorization of the non-singular matrix G. It is now clear what should be done. Suppose that • is complex shift and that • is forced to be •¤ , then G will be real, U X will be its (real) orthogonal factor, $ will be real, and the computation of complex will be circumvented. The solution is therefore to use the generalized QR transformation with the providing every shift with positive imaginary part is followed by its conjugate, the intermediate complex matrix being avoided as indicated above. Ignoring the storage of the auxiliary matrix G it turns out the formation of G, its decomposition, and the formation of $ (one double steps) require at least as much as work as two real QR steps. Note that apart from the decided advantage of keeping the arithmetic real there no extra gain.

Conclusion
We have seen that calculating eigenvalues a matrix is an important problem in mathematics and the science but the native approach of solving characteristic polynomial is inefficient for large dimensional matrices. Rather, using the elegant QR methods provides an effective answer. The QR method is one of the most wildly used methods which used to find all eigenvalues of square matrix. Therefore, we are intended to determine all eigenvalues of real square matrix by using QR factorization method. Suppose a square matrix A is given. Now, let A = Q 0 R 0 be QR factorization of A, and then create A 1 = R 0 Q 0 . Let A 1 = Q 1 R 1 be QR factorization of A and create A 2 = R 1 Q 1 . Continue this process again in the same fashion. For ≥ 1, ( ℎ = ). Once A m has been created such that, A m = Q m R m , and A m+1 = R m Q m . Thus, the sequence {A m } will usually converges to something from which the eigenvalues can be computed easily. Moreover, A 2 is similar to A and A m+1 is similar to A m and so on. Hence A m and A m+1 have the same eigenvalues by definition of similarity of matrices. This process will be stopped when the entries below the main diagonal of current matrix A m are sufficiently small, or if it appears that convergence will not happen. This implies that, QR factorization sequence process can fail to converge or the convergence can be also extremely slow and hence expensive. If it converges to certain matrix, then the diagonal entries of this current matrix are tends to eigenvalues of matrix A, if not, it can be modified in order to improve convergence/ to accelerate the rate of convergence dramatically by using shifting methods.