12.4.6. Iterative QR Algorithm

The iterative QR algorithm is a well-known algorithm for finding the eigenvalues of a matrix. We are specifically referring to the explicitly shifted QR algorithm. The iterative algorithm relies on the QR matrix factorization (QR Decomposition) to move a matrix from upper Hessenberg form to upper triangular form.

We begin by describing the general concept of the algorithm before introducing refinements. First, convert the matrix to upper Hessenberg form by similarity transforms. Next, use QR factorization, \(\mathbf{A} = \mathbf{Q\,R}\), to find a new matrix \(\mathbf{A}_1 = \mathbf{R\,Q}\). The original \(\bf{A}\) matrix and the modified \(\mathbf{A}_1\) matrix are similar. Since \(\mathbf{A} = \mathbf{Q\,R}\), then \(\mathbf{R} = \mathbf{Q}^H\,\mathbf{A}\) and \(\mathbf{R\,Q} = \mathbf{Q}^H\,\mathbf{A\,Q}\) establishes similarity between \(\m{A}_i\) and \(\m{A_{i+1}} = \m{R}_i\, \m{Q}_i\). QR factorization is repeated until \(\mathbf{A}_k\) is upper triangular, and additional iterations will not change the result. As the \(\mathbf{A}_k\) matrix converges toward an upper triangular matrix with each iteration, the \(\bf{Q}\) matrix converges toward the identity matrix.

Two algorithm refinements lead to faster convergence. The first refinement reduces the diagonal by a scalar shift, \(\rho\), before each QR decomposition. As shown in equation (12.7), the shift is added back to the diagonal values after the decomposition. When the next eigenvalue to find is real, the value of the last diagonal element works for the shift. However, a Wilkinson shift is better for quicker convergence and for finding complex eigenvalues. Calculate the two eigenvalues of the \(2{\times}2\) lower-right submatrix. The eigenvalue closest to the last diagonal value is selected as the shift [WILKINSON65], [WILKINSON68shift]. Convergence often occurs with the next QR step after finding an eigenvalue from the \(2{\times}2\) submatrix. The algorithm finds eigenvalue \(\lambda_k\) at \(\m{A}(k, k)\) when the absolute value of the element at \(\m{A}(k, k-1)\) converges to near zero. The second algorithm refinement is deflation, which reduces the size of the matrix to exclude the last row after finding each eigenvalue.

(12.7)\[\begin{split}\begin{aligned} &\mathbf{A}_{j-1} - \rho\,\mathbf{I} = \mathbf{Q\,R} \\ &\mathbf{A}_j = \mathbf{R\,Q} + \rho\,\mathbf{I} \end{aligned}\end{split}\]

Note

The last row and column can be removed with deflation if we only want the eigenvalues. However, we can only remove the last row if we use the Schur vectors in the calculation of the eigenvectors. The multiplications performed during the remaining QR decompositions must affect every element above the last row. Thus, deflation causes us to perform the QR decomposition on under-determined matrices.

The eigenvectors may then be found after the Schur decomposition is finished from the Schur vectors and the eigenvectors of the upper triangular matrix, which are computed using the procedure in Eigenvectors Following the Schur Decomposition.

12.4.7. The eigQR Function

The eigQR function is a simple implementation of the iterative QR algorithm that first finds the Schur decomposition, and then finds the eigenvalues and the eigenvectors.

function [X, Lambda] = eigQR(A)
% EIGQR - Iterative QR eigen{value, vector} algorithm.
%         This is the explicitly-shifted QR algorithm.
%
% A - Square matrix
% EIGQR(A) or L = EIGQR(A) returns the eigenvalues of A.
% [X, L] = EIGQR(A) returns the eigenvectors of A as the
%     columns of the X matrix and the eigenvalues of A
%     are along the diagonal of matrix L.
%
% NOTE: This code is for educational purposes.
% For other purposes, use MATLAB's eig function.
%
% See also: EIG, EIG22, HESS, QR, EIGFRANCIS, UTEIGVECT

    [m, n] = size(A);
    if m ~= n
      error("Matrix must be square to find eigenvalues")
    end
    % Start with Hessenberg - zeros below first diagonal
    [Q,T] = hess(A); tol = 1e-14; eigs = zeros(n,1);
    k = n; I = eye(n);
    while k > 1
        [rho1, rho2] = eig22(T(k-1:k,k-1:k));
        if abs(T(k,k)-rho1) < abs(T(k,k)-rho2)
            shift = rho1*I; else, shift = rho2*I;
        end
        % MATLAB's qr is faster, but qrFactor allows for
        % breakpoints to understand the algorithm better.
        B = T(1:k,1:n) - shift; [Qi, R] = qrFactor(B); % qr(B);
        Qn = eye(n); Qn(1:k,1:k) = Qi;
        T(1:k,1:n) = R*Qn + shift; Q = Q*Qn;
        % Converge when kth subdiagonal near zero
        if abs(T(k,k-1)) < ...
                (tol*(abs(T(k,k)) + abs(T(k-1,k-1))))
            T(k,k-1) = 0; eigs(k) = T(k,k); % kth eigenvalue
            k = k - 1; I = eye(k, n);       % deflation
        end
    end
    eigs(1) = T(1,1);
    if nargout <= 1 % only want eigenvalues
        X = sort(eigs);
    else            % eigenvectors also
        Lambda = diag(eigs); V = UTeigVect(T); X = Q*V;
    end
end

12.4.8. The eig22 Function

The eigenvalues of a \(2{\times}2\) matrix are unique as the only eigenvalues found from a direct calculation–the quadratic equation factoring of its characteristic polynomial. The eig22 function is called many times by the iterative algorithms for diagonalizing larger matrices.

function [lambda1, lambda2] = eig22(A)
% EIG22 - Quick eigenvalue calculation of a 2-by-2 matrix
% [lambda1, lambda2] = EIG22(A) returns two eigenvalues

    m = (A(1,1) + A(2,2))/2;
    p = A(1,1)*A(2,2) - A(1,2)*A(2,1);
    s = sqrt(m^2 - p); lambda1 = m + s; lambda2 = m - s;
end