Symmetric QR Algorithm with Permutations

In this paper, we present the QR Algorithm with Permutations that shows an improved convergence rate compared to the classical QR algorithm. We determine a bound for performance based on best instantaneous convergence, and develop low complexity methods for computing the permutation matrices at every iteration. We use simulations to verify the improvement, and to compare the performance of proposed algorithms to the classical QR algorithm.


I. INTRODUCTION
Let A be an NxN symmetric invertible matrix, spectral theorem states that the matrix can be diagonalised as follows: where Λ is a matrix with eigenvalues along the diagonals and V = [v 1 , v 2 , . . . , v N ] is a matrix with corresponding eigenvectors v 1 , v 2 , . . . , v N . QR algorithm is an iterative algorithm based on QR decomposition to find the eigenvalues and eigenvectors simultaneously and is as follows: Further, A ∞ = Λ, Q ∞ = I, and V = Q 0 Q 1 ...Q ∞ . For ease of analysis, we require in equation (3) that uppertriangular matrix R of QR factorization to contain positive values along the diagonal. This can be easily arranged by multiplying the row containing the negative value by '-1' and compensating the matrix Q. Such upper-triangular matrix R can be shown to be the upper-triangular Cholesky factor of the matrix A 2 .
For an introduction to the QR algorithm and a proof for convergence, and for modified algorithms: QR algorithm based on Hessenberg form and QR algorithm with shifts, refer [1] or [3].

A. Cholesky Iteration
A closely related algorithm to QR algorithm is Cholesky iterations based on Cholesky decomposition, given as follows: Aravindh Krishnamoorthy is currently with Ericsson Modem Nuremberg GmbH working in the area of Wireless Communication. Contact e-mail: aravindh.krishnamoorthy@ericsson.com, aravindh.k@ieee.org.
The sets {A 0 , A 1 , A 2 , ...} and {B 0 , B 1 , B 2 , ...} are equivalent sets, and A k is the principal square root of the matrix B k . This relation can be used to study the convergence properties of QR algorithm through Cholesky iteration which is more tractable for analysis; for e.g. see [4].

B. Permutations
A permutation matrix P is a matrix containing a single '1' in every row and column, and zeroes elsewhere. Let P be the set of such matrices of order N, then the cardinality of set P is |P| = N ! A symmetric permutation P AP H on matrix A defined above, maintains the symmetry of the matrix. It can be easily verified that the permuted matrix has the same eigenvalues as the original matrix, and the eigenvectors are P V .

C. Terminology
The following special functions are used in this paper: • dsort : R N → R N takes a real vector of order N as input, and returns a vector with the elements of input sorted in descending order. • diag : C N xN → R N takes a symmetric matrix of order N as input, and returns a vector containing the diagonal elements of the matrix in the order along the main diagonal. • diagmat : R N → R N xN takes a real vector of order N as input, and returns a real matrix containing the elements of the input vector along the main diagonal, and zeroes elsewhere.
In the following section, we present the QR algorithm with permutations along with a theoretical bound on performance. In section III we discuss computation of the permutation matrices, in section IV we describe the simulation setup used to verify the performance of the algorithms. In section V we discuss the performance results for simulation runs involving positive-definite and symmetric matrices, followed by concluding remarks in section VI.

II. QR ALGORITHM WITH PERMUTATIONS
Applying a permutation P on equations (2, 3, 4) yields an iteration of QR algorithm as follows: showing that a permutation on the matrix A has an additional factor P H in the applied orthonormal matrix.

arXiv:1402.5086v1 [cs.NA] 20 Feb 2014
The effect of applying a permutation matrix on convergence can be intuitively understood when the QR algorithm is viewed as power iteration algorithm with a complete set of orthonormal basis-vectors instead of a single vector.
Let n-th column of the matrix Q k be q n = a n,1 v 1 +a n,2 v 2 + ... + a n,N v N , the convergence rate of the n-th eigenvalue λ n depends upon a n,n /a n,n−1 λ n /λ n−1 . With some luck, permutation changes the values a n,1 , a n,2 , ..., a n,N , such that it offers a better convergence compared to the original matrix.
Assume that P is a permutation matrix that can ensure an optimum arrangement of the basis-vectors of Q, then the permutation causes the maximum convergence in that QR algorithm's iteration. This optimal arrangement can also aid the convergence of repeated eigenvalues or eigenvalues close to each other. However, this convergence is only speeded up by a linear factor! Let A be an NxN symmetric invertible matrix as defined above, the QR algorithm with permutations with a permutation matrix P k at each iteration is given as follows: Further, A ∞ = Λ, Q ∞ = I, and In general, in the k-th iteration, the estimates of eigenvalue matrix Λ k , eigenvector matrix V k can be given as follows: The error in eigenvalue estimation E k for any algorithm can be given as follows: QR algorithm with permutations requires the following additional computation compared to the classical QR algorithm: 1) computation of the permutation matrix, 2) symmetric permutation of matrix A k before QR decomposition, 3) permutation of the matrix Q k before multiplication, if eigenvectors are desired.

A. Performance Bound
The performance of the algorithm depends on the selection of the permutation matrices, i.e. E k = f (A, P 0 , P 1 , ..., P k ).
A simple bound for E k can be given by requiring the best instantaneous convergence at every iteration. Therefore, the selection criterion for selecting the permutation matrix in the k-th iteration becomes: That is, the permutation matrix P k of the available permutation matrices of order N, which when used in the k-th iteration, causes the least error in eigenvalue estimation E k .

III. COMPUTING THE PERMUTATION MATRICES
Computing the permutation matrix P k is non-trivial. In this section we present two simple methods for computing the permutation matrices which compute them without iterating through the set P. Other methods for computing P k may be established easily, and their choice is a trade-off between convergence rate and computational complexity.

A. Diagonal Ordering
In this method, we choose the permutation matrix P k which orders the diagonal elements matrix A k in descending order of magnitude (absolute value). The MATLAB compatible pseudo-code to determine this permutation matrix is as follows:

B. Column Ordering
In this method, we choose the permutation matrix P k which orders the diagonal elements matrix B k = A 2 k in descending order. Such ordering also orders the columns of the matrix A k in descending order of their norm value.
This permutation matrix can be computed with a computational cost of N 2 for computing the N norm values of the columns of A k . The MATLAB compatible pseudo-code to determine this permutation matrix is P = diord(A * A).
Both diagonal ordering and column ordering methods do not preserve the shape of the matrix in subsequent QR iterations. Therefore, computational complexity reduction due to conversion to Hessenberg form cannot be applied with these two methods.
IV. SIMULATION SETUP In order to compare the performance of various algorithms, we use MATLAB based simulations. MATLAB compatible source code of the simulation is available from [7].
The average eigenvalue estimation error E 2 k equation (16) for 4x4 symmetric invertible matrices, each with 50 iterations are plotted. The first figure (1a) uses a random subset of size 10,000 of the set of symmetric matrices which includes positive-definite, negative-definite, and indefinite matrices. The second figure (1b) uses a random subset of size 25,000 of the set of positive-definite symmetric matrices.
Label 'QR' is the error plot for the classical QR algorithm, 'QRH' first converts the given matrix into Hessenberg (tridiagonal) form and then performs the classical QR algorithm, 'QRS' uses the QRH algorithm (A is in Hessenberg form) with shifts as given below (A 0 = A): Here the notation A k−1 (N, N ) refers to the (N,N)-th element of the matrix A k−1 .
Label 'DO' refers to the QR algorithm with permutations which uses diagonal ordering method for computing the permutation matrix at each iteration. 'CO' uses the column ordering method.
Label 'BIC' refers to the best instantaneous convergence bound as described in equation (17).

V. SIMULATION RESULTS
From figure (1a), we observe that permutations increase the convergence rate on an average for symmetric and positive definite matrices.
The best instantaneous convergence (BIC) bound shows the highest rate of convergence, showing that permutation has an effect of speeding up the convergence of QR algorithm. However, meeting this bound using the definition is not feasible as it requires prior knowledge of eigenvalues.
QR algorithm with permutations with diagonal ordering (DO) shows good convergence rate for positive-definite matrices, but saturates after the first few iterations on an average, for symmetric matrices in general. DO performs better than classical QR algorithm and variations for positive-definite matrices.
QR algorithm with permutations with column ordering (CO) shows good convergence rate for positive-definite matrices and in general for symmetric matrices. CO performs better than classical QR algorithm and variations for positive-definite matrices and in general for symmetric matrices.
Therefore, for positive-definite matrices, DO provides the best convergence rate, and for symmetric matrices in general, CO provides the best convergence at low complexity for calculation of permutation matrices, both providing a 2x convergence speedup.

VI. CONCLUSION
In this paper, we presented the QR algorithm with permutations that shows a significantly improved convergence rate compared to classical QR algorithm and variations. We determined a bound for performance based on best instantaneous convergence, and developed two low complexity methods for computing the permutation matrices at every iteration. We used simulations to verify that the performance bound is significantly improved upon using permutations, and compared the performance of QR algorithm with permutations using the two methods for computing permutation matrices.
We find that a 2x convergence speedup is obtained by using QR algorithm with permutations. Further, we find that the diagonal ordering provides the best trade-off for positivedefinite matrices, and column ordering provides best trade-off for symmetric matrices in general.