QR decomposition & QR algorithm

Gypsophila
  2025-03-28 22:18:40 2025-03-28 22:18:40 Created   2025-03-30 17:35:33 2025-03-30 17:35:33 Updated    

QR decomposition

We give a brief introduction to QR decomposition, which is a method for decomposing a matrix into an orthogonal matrix and an upper triangular matrix. The QR decomposition is useful in various applications, including solving linear systems, least squares problems, and eigenvalue problems.

The implementation of QR decomposition can be done using various methods, including Classical Gram-Schmidt orthogonalization, modified Gram-Schmidt, Householder transformation, and Givens rotation. Each method has its own advantages and disadvantages in terms of numerical stability and computational efficiency.

Classical Gram-Schmidt Orthogonalization

CGS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
% CGS. Classical Gram-Schmidt Process
% Input:    V = [v1, v2, ..., vm]
% Output:   Q = [q1, q2, ..., qm]
% s.t.  Q'*Q = I, span(Q) = span(V)
%       V = Q*R
%       R is an upper triangular m x m matrix.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.4, P229.

function [Q, R] = CGS(V)
[n, m] = size(V);                               % assume n >= m
if n < m
    fprintf("column number is lager than row number!")
    return
end

Q = zeros(n, m);
R = zeros(m, m);

%%%%%%%%%%%%%%%%%%%% Classical Cram-Schmidt Process %%%%%%%%%%%%%%%%%
%---------------------- specify the first column -------------------%
R(1, 1) = norm(V(:, 1));

% Check whether linear dependent to aovid 0 as denominator
if R(1, 1) < 2^(-50)
    fprintf("columns of V are linearly dependent!")
    return
end

Q(:, 1) = V(:, 1) / R(1, 1);                    % normalization

%--------------------- deal with following columns -----------------%
% CGS
for j = 2:m
    Q(:, j) = V(:, j);                          % Orthogonalization initializes
    for i = 1:j-1                               % Orthogonalization begins
        R(i, j) = V(:, j)' * Q(:, i);           % using Unorthogonalized V(:,j)
        Q(:, j) = Q(:, j) - R(i, j) * Q(:, i);  % delete projections on previous vectors
    end
    R(j, j) = norm(Q(:, j));                    % diagonal entry

    % Check whether linear dependent to aovid 0 as denominator
    if R(j, j) < 2^(-50)
        fprintf("columns of V are linearly dependent!")
        return
    end

    Q(:, j) = Q(:, j) / R(j, j);                % normalization
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Modified Gram-Schmidt Orthogonalization

MGS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
% MGS. Modified Gram-Schmidt Process
% Input:    V = [v1, v2, ..., vm]
% Output:   Q = [q1, q2, ..., qm]
% s.t.  Q'*Q = I, span(Q) = span(V)
%       V = Q*R
%       R is an upper triangular m x m matrix.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.4, P237.

function [Q, R] = MGS(V)
[n, m] = size(V);                               % assume n >= m
if n < m
    fprintf("column number is lager than row number!")
    return
end

Q = zeros(n, m);
R = zeros(m, m);

%%%%%%%%%%%%%%%%%%%% Modified Cram-Schmidt Process %%%%%%%%%%%%%%%%%
%---------------------- specify the first column -------------------%
R(1, 1) = norm(V(:, 1));

% Check whether linear dependent to aovid 0 as denominator
if R(1, 1) < 2^(-50)
    fprintf("columns of V are linearly dependent!")
    return
end

Q(:, 1) = V(:, 1) / R(1, 1);                    % normalization

%--------------------- deal with following columns -----------------%
% MGS
for j = 2:m
    Q(:, j) = V(:, j);                          % Orthogonalization initializes
    for i = 1:j-1                               % Orthogonalization begins
        %----------- Main Change compared with CGS ---------------%
        R(i, j) = Q(:, j)' * Q(:, i);           % using *partially* orthogonalized Q(:,j)
        %---------------------------------------------------------%
        Q(:, j) = Q(:, j) - R(i, j) * Q(:, i);  % delete projections on previous vectors
    end
    R(j, j) = norm(Q(:, j));                    % diagonal entry

    % Check whether linear dependent to aovid 0 as denominator
    if R(j, j) < 2^(-50)
        fprintf("columns of V are linearly dependent!")
        return
    end

    Q(:, j) = Q(:, j) / R(j, j);                % normalization
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Householder transformation

HO.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
% HO. Orthogonalization by Householder Transform
% Input:    V = [v1, v2, ..., vm]  (n x m matrix)
% Output:   vectors tau and gamma, matrix W
% s.t.  V = Q*R, where Q is implictly saved as multiplication of some
% Householder reflectors, given by tau, gamma and uk's (stored in *strictly
% lower triangular part* of W: uk = [1; W(k+1:n, k)]
%       Qk = I - gamma(k) * uk * uk'
%       Q = Q1 * Q2 * ... * Q_min(m, n-1)
%       Q'*Q = I, span(Q) = span(V)
%       R is an upper triangular matrix, which is stored in *upper triangular
%       part* of W.
% Remark. Explicit construction of Q need additional 4m^3/3 flops, since we
% almost always do not need the matrix Q, for example, we only need to apply Qk
% to both sides of equation and then solve the final triangular system to
% solve a linear system, we only need keep the two vector tau and gamma for
% forming Qk's.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P205.
function [tau, gamma, W] = HO(V)
    [n, m] = size(V);                               % assume n >= m
    if n < m
        fprintf("column number is lager than row number!")
        return
    end
    
    % Initialization
    tau = zeros(m, 1);
    gamma = zeros(m, 1);

    %===================== Orthogonalization =====================%
    for k = 1:min(m, n-1)
        %------ construct Householder Transform reflector --------%
        [tau(k), gamma(k), V(k:end, k)] = reflector(V(k:end, k)); % orthogonalize kth column of V  
        %---------------------------------------------------------%

        %---------- Update k+1 to m columns of V -----------------%
        V(k:end, k+1:end) = V(k:end, k+1:end) ...                 % V(k:m, k+1:end) = Qk * V(k:m, k+1:end)
                          - gamma(k) * V(k:end, k) * V(k:end, k)' *V(k:end, k+1:end);
        %---------------------------------------------------------%

        V(k, k) = -tau(k);              % update the diagonal entry of V
    end
    %=============================================================%

    if m == n                           % check whether orthogonalized the last column
        gamma(n) = V(n, n);
    end

    W = V;                              % W is stored in V
end

% REFLECTOR. construct Householder Transform reflector
% Input:    vector x with length n
% Output:   parameters tau, gamma and vector u
% s.t.  Q = (I - gamma * u*u') is a reflector for which
%       Qx = [-tau, 0, ..., 0]'
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P201.
function [tau, gamma, u] = reflector(x)
    beta = norm(x, "inf");              % check whether x=0
    if beta < 2^(-52)                   % if x=0, no need to reflect x
        gamma = 0;                      % gamma is set to be zero, giving Q=I
        tau = 0;
    else
        % ------------ x!=0 ------------%
        x = x / beta;                   % normalization

        tau = norm(x);                  % Q is unitary, norm(Qx) = tau = norm(x)

        if x(1) < 0                     % check if sign of x1 is same with tau
            tau = - tau;                % reverse sign of tau to aovid cancellation
        end
        x(1) = tau + x(1);              % same sign -> without cancellation

        gamma = x(1) / tau;             % recall gamma = (tau + x1) / tau

        x(2:end) = x(2:end) / x(1);     % [u2, ..., un] = [x2, ..., xn] / (tau + x1)
        x(1) = 1;                       % u1 = 1

        tau = tau * beta;               % re-scaling
    end
    u = x;                              % u is stored in x
end
HOqr.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
% HOQR. QR decompostion by Householder Transform
% Input:    V = [v1, v2, ..., vm]  (n x m matrix)
% Output:   Q = [q1, q2, ..., qm], R
% s.t.  Q'*Q = I, span(Q) = span(V)
%       V = Q*R
%       R is an upper triangular n x m matrix.
% Remark. Explicit Q is not necessary for most applications. See HO.m
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P205.
function [Q, R] = HOqr(V)
    [n, m] = size(V);

    [~, gamma, W] = HO(V);                  % using HO.m

    R = triu(W);                            % R stored in upper part of W
    
    Q = eye(n);                             % Initialization
    for k = 1:min(m, n-1)
        uk = [1; W(k+1:end, k)];            % Reflect vector

        % Q = Q1 * Q2 * ... * Q_min(m, n-1)
        Q(:, k:end) = Q(:, k:end) - gamma(k)* (Q(:, k:end) * uk) * uk';
    end

end

Givens rotation

GR.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
% GR. Orthogonalization by Givens Rotation
% Input:    V = [v1, v2, ..., vm]  (n x m matrix)
% Output:   W
% s.t.  V = Q*R, where Q is implictly saved as multiplication of some Givens 
% rotations, given by c(i,j)'s (stored in *strictly lower triangular part* of 
% W: c(i,j) = W(i,j), i>j
%       Qij is an identity matrix, except for an 2x2 submatrix to be
%       [c(i,j), -s(i,j); s(i,j), c(i,j)], where s(i,j) = sqrt(1-c(i,j)^2)
%       by applying Qij' to V(:, j), we introduce zero entry at V(i,j).
%       R = Qn1' * Q(n-1)1' * ... * Q21' * Qn2' * ... * Qnm' * ... *
%       Q(m+1)m' * V, therefore Q = Q(m+1)m * ... * Qn1, Q'*Q = I
%       R is an upper triangular matrix, which is stored in *upper triangular
%       part* of W.
% Remark. Explicit construction of Q need additional 4m^3/3 flops, since we
% almost always do not need matrix Q, for example, we only need to apply Qij
% to both sides of equation and then solve the final triangular system to
% solve a linear system, we only need keep the cij for forming Qij's.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P190.
function W = GR(V)
    [n, m] = size(V);                               % assume n >= m
    if n < m
        fprintf("column number is lager than row number!")
        return
    end

    %===================== Orthogonalization =====================%
    for j = 1:min(m, n-1)
        for i = n:-1:j+1
            %------------ construct Givens Rotator ---------------% 
            % using (i-1,j) and (i,j) elements       -------------> rotate in {x(i-1)-xi} plane
            [c, s, V(i-1, j)] = rotator([V(i-1, j); V(i,j)]);     % Introduce 0 at (i,j)
            Q = [c, -s; s, c];
            %--- Update j+1 to m columns of V([i-1,i], :) --------%
            V([i-1, i], j+1:end) = Q'*V([i-1, i], j+1:end);
            V(i, j) = c;                                          % restrore c
        end
    end
    %=============================================================%
    W = V;                                      % W is stored in V
end

% ROTATOR. construct Givens Rotator
% Input:    vector x = [x1; x2]
% Output:   parameters c, s, nu
% s.t.  Q = [c, -s; s, c]
%       Q'x = [nu; 0], nu = norm(x, 2)
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P195.
function [c, s, nu] = rotator(x)
    beta = norm(x, "inf");              % check whether x=0
    if beta < 2^(-52)                   % if x=0, no need to rotate x
        c = 1; s = 0;                   % Q=I
    else
        % ------------ x!=0 ------------%
        x(1) = x(1) / beta;             % normalization
        x(2) = x(2) / beta;             % scaling to aovid overflow or harmful underflow

        nu = norm(x);                   % Q is unitary, norm(Qx) = tau = norm(x)

        c = x(1) / nu;
        s = x(2) / nu;
    end

    nu = beta * nu;                     % re-scraling
end
GRqr.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
% GRQR. Orthogonalization by Givens Rotation
% Input:    V = [v1, v2, ..., vm]  (n x m matrix)
% Output:   Q = [q1, q2, ..., qm], R
% s.t.  Q'*Q = I, span(Q) = span(V)
%       V = Q*R
%       R is an upper triangular n x m matrix.
% Remark. Explicit Q is not necessary for most applications. See GR.m
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P190.
function [Q, R] = GRqr(V)
    [n, m] = size(V);

    W = GR(V);                  % using GR.m
    R = triu(W);                % R stored in upper part of W

    C = tril(W, -1);            % c in rotators
    S = sqrt((1-C).*(1+C));     % s in rotators

    Q = eye(n);                 % Intialization
    for j = 1:min(m, n-1)
        for i = n:-1:j+1
            c = C(i,j);
            s = S(i,j);
            Qij = [c, -s; s, c];% rotator
            % Q = Qn1*Q(n-1)1*...*Q21 * Qn2*Q(n-1)2*...*Q32 * ...
            Q(:, [i-1, i]) = Q(:, [i-1, i]) * Qij;
        end
    end

end

QR algorithm

upperHenssenberg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
% UPPERHENSSENBERG. Reduce a given matrix to upper Henssenberg form
% Input:    matrix A
% Output:   matrix H, vector gamma
% s.t.  triu(H, -1) is the upper Henssenberg matrix reduced from A
%       uk = [1; H(k+1:n, k)] is the k-th reflect vector
%       Qk = diag(I(k-1), (I - gamma(k) * uk * uk')) is the k-th reflector
%            applied to A
%       triu(H, -1) = Q(n-2)'*...*Q2'*Q1' *A* Q1*Q2*...*Q(n-2) = Q'*A*Q
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.5, P354.
function [H, gamma] = upperHenssenberg(A)
    [n, ~] = size(A);

    if n <= 2                       % there is no need to do anything if n <= 2
        H = A;                      % already upper Henssenberg
        return
    end
    
    gamma = zeros(n-2, 1);                              % initalization
    %================= Reduction by Householder Reflector =================%
    for k = 1:n-2
        [tau, gamma(k), u] = reflector(A(k+1:end, k));  % Set up k-th reflector
        A(k+1:end, k) = u;          % save reflect vector in lower triangular part
        
        %----------- Left multiplication by reflector Qk'=Qk---------------%
        A(k+1:end, k+1:end) = A(k+1:end, k+1:end) - gamma(k)* u * (u' * A(k+1:end, k+1:end));
        %------------- Right multiplication by reflector Qk ---------------%
        A(:, k+1:end) = A(:, k+1:end) - gamma(k)* (A(:, k+1:end) * u) * u';
        %------------------------------------------------------------------%

        % replace first 1 in uk by (k+1, k) entry of upper Henssenberg matrix
        A(k+1, k) = -tau;
    end
    %=======================================================================%
    % tau = -A(n, n-1);

    H = A;                      % H is stored in A
end

% REFLECTOR. construct Householder Transform reflector
% Input:    vector x with length n
% Output:   parameters tau, gamma and vector u
% s.t.  Q = (I - gamma * u*u') is a reflector for which
%       Qx = [-tau, 0, ..., 0]'
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P201.
function [tau, gamma, u] = reflector(x)
    beta = norm(x, "inf");              % check whether x=0
    if beta < 2^(-52)                   % if x=0, no need to reflect x
        gamma = 0;                      % gamma is set to be zero, giving Q=I
        tau = 0;
    else
        % ------------ x!=0 ------------%
        x = x / beta;                   % normalization

        tau = norm(x);                  % Q is unitary, norm(Qx) = tau = norm(x)

        if x(1) < 0                     % check if sign of x1 is same with tau
            tau = - tau;                % reverse sign of tau to aovid cancellation
        end
        x(1) = tau + x(1);              % same sign -> without cancellation

        gamma = x(1) / tau;             % recall gamma = (tau + x1) / tau

        x(2:end) = x(2:end) / x(1);     % [u2, ..., un] = [x2, ..., xn] / (tau + x1)
        x(1) = 1;                       % u1 = 1

        tau = tau * beta;               % re-scaling
    end
    u = x;                              % u is stored in x
end

QR Iteration

eqr.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% EQR. Explicit QR Iteration Algorithm
% Input:    A    - matrix concerned
%           maxN - max iteration number
%           tol  - tolance
% Output:   M    - approximate Schur form of A
% -------------------------------------------------------------------------
% Theorem. Let A <- Q'*A*Q = R*Q, where [Q, R] = qr(A), then A converges to
% Schur form of A if all the eigenvalues have different absolute values.
function M = eqr(A, maxN, tol)
    M = A;
    n = 1;
    e = 1;
    while n <= maxN && e > tol
        [Q, R] = qr(M);         % Explicitly compute the Q matrix
        M = R*Q;                % M <- Q'*M*Q = R*Q
        n = n+1;
        e = norm(tril(M));
    end
end
eqr1S.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
% EQR1S. Explicit QR Iteration Algorithm with a Shift
% Input:    A    - matrix concerned
%           maxN - max iteration number
%           tol  - tolance
% Output:   M    - approximate Schur form of A
% -------------------------------------------------------------------------
% Theorem. Let A <- Q'*A*Q = R*Q, where [Q, R] = qr(A), then A converges to
% Schur form of A if all the eigenvalues have different absolute values.
function M = eqr1S(A, maxN, tol)
    M = A;
    [n, ~] = size(M);
    t = 1;                              % counter
    e = 1;                              % error
    while t <= maxN && e > tol
        shift = wilkinson(M);           % Wilkinson shift
        % can also using Rayleigh quotient shift
        % shift = M(n, n);

        [Q, R] = qr(M - shift*eye(n));  % Explicitly compute the Q matrix
        M = R*Q + shift*eye(n);         % M <- Q'*M*Q = R*Q
        t = t+1;
        e = norm(tril(M));
    end
end

function shift = wilkinson(A)
    [n, ~] = size(A);
    htr = .5*(A(n-1, n-1) + A(n, n));   % half a 2x2 trace
    dscr = sqrt((.5*(A(n-1,n-1)-A(n,n)))^2 + A(n,n-1)^2);
    % discriminant
    if htr < 0
        dscr = -dscr;                   % to avoid cancellation
    end 
    root1 = htr + dscr;                 % quadratic formula
    if root1 == 0                       % almost never happens
        root2 = 0;
    else                                % almost always happens
        det = A(n-1,n-1)*A(n,n) - A(n,n-1)^2;
        % 2x2 determinant = product of roots
        root2 = det/root1;
    end
    if abs(A(n,n)-root1) < abs(A(n,n)-root2)
        shift = root1;
    else
        shift = root2;
    end
end

Francis’s Implicitly Shifted QR Algorithm

Single shift for Symmetric Matrix

iqr1S4S.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
% IQR1S4S. Implicitly QR Iteration Algorithm with single shift for Symmetric Matrix
% Input:    symmetric matrix A
% Output:   matrix M
% s.t.  By using Francis's Implicit Iteration, M(n,n) will become extremely  
% close to an eigenvalue of A, and M(n,n-1) is nearly 0 as iteration continues.
% This approach is applied repeatly to compute all eigenvalues of A.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P359.
% Test example: n=6;row=[2,1,zeros(1,n-2)];A=toeplitz(row);
function M = iqr1S4S(A)
    [n, ~] = size(A);

    %--------- Reduction to Upper Henssenberg form -----------------%
    A = upperHenssenberg(A);    % call upperHenssenberg.m
    A = triu(A, -1);            % reduce A to upper Henssenberg form
    %---------------------------------------------------------------%

    M = A;          % initalization
    maxN = 100;     % max iteration number

    % Reduce (-1)-subdiagonal entries from (n, n-1) to (2, 1)
    for m = n:-1:2
        % using Implicitly QR Iteration to reduce M(m, m-1) to 0
        % in next loop only need to consider a smaller (n-1)x(n-1) martix
        if abs(M(m, m-1)) < 1e-16 * (abs(M(m, m)) + abs(M(m-1, m)))
            M(m, m-1) = 0;      % check whether *proper* upper Henssenberg
        else
            M(1:m, 1:m) = deflation(M(1:m, 1:m), maxN);
        end
    end
end

% DEFLATION. Eliminate A(n, n-1) to make A(n,n) approach an
% eigenvalue of A by using Francis's Implicit Iteration.
% Input:    symmetric matrix A, max iteration number maxN
% Output:   matrix M
function M = deflation(A, maxN)
    [n, ~] = size(A);
    M = A;
    t = 1;                                      % counter
    while t <= maxN && abs(M(n, n-1)) > 1e-16   % stop if M(n,n-1) is small enough
        M = SymmetricFrancisIteration(M);       % Francis's Iteration
        t = t+1;
    end
end

% SYMMETRICFRANCISITERATION. Perform a symmetric Francis's Iteration.
% Input:    symmetric matrix A
% Output:   matrix M
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P362.
function A = SymmetricFrancisIteration(A)
    [n, ~] = size(A);

    % assign the single shift to be Wilkinsion shift
    shift = wilkinson(A);
    % can also using Rayleigh quotient shift
    % shift = M(n, n);

    %----------- Construct Givens rotator -----------%
    cs = A(1, 1) - shift;
    sn = A(2, 1);
    % normalization
    r = norm([cs, sn]);
    cs = cs / r;
    sn = sn / r;
    q0 = [cs, -sn; sn, cs];     % Givens rotator
    %------------------------------------------------%
    
    %----------- Introduce a bulge ------------------%
    A(1:2, :) = q0'*A(1:2, :);  % left multiplication
    A(:, 1:2) = A(:, 1:2)*q0;   % right multiplication
    %------------------------------------------------%

    %================= Bulge-Chasing ================%
    for i = 1:n-2
        %--------- Construct Givens rotator ---------%
        cs = A(i+1, i);
        sn = A(i+2, i);
        % normalization
        r = norm([cs, sn]);
        cs = cs / r;
        sn = sn / r;
        qi = [cs, -sn; sn, cs];% Givens rotator
        %--------------------------------------------%

        %---------- eliminate the old bulge----------%
        A(i+1, i) = r;
        A(i+2, i) = 0;
        %--------------------------------------------%

        %--introduce a new bulge along (-1)-subdiag--% Similarity Transform
        A(i+1:i+2, i+1:n) = qi' * A(i+1:i+2, i+1:n); % Left multiplication
        A(:, i+1:i+2) = A(:, i+1:i+2) * qi;          % Right multiplication
        %--------------------------------------------%
    end
    %================================================%
end

% WILKINSON. Compute Wilkinson shift with respect to matrix A
% Input:    matrix A
% Output:   shift
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P361.
function shift = wilkinson(A)
    [n, ~] = size(A);
    htr = .5*(A(n-1, n-1) + A(n, n));   % half a 2x2 trace
    dscr = sqrt((.5*(A(n-1,n-1)-A(n,n)))^2 + A(n,n-1)^2);
    % discriminant
    if htr < 0
        dscr = -dscr;                   % to avoid cancellation
    end 
    root1 = htr + dscr;                 % quadratic formula
    if root1 == 0                       % almost never happens
        root2 = 0;
    else                                % almost always happens
        det = A(n-1,n-1)*A(n,n) - A(n,n-1)^2;
        % 2x2 determinant = product of roots
        root2 = det/root1;
    end
    if abs(A(n,n)-root1) < abs(A(n,n)-root2)
        shift = root1;
    else
        shift = root2;
    end
end

The above code is for symmetric matrix. For general real matrix, the following code may be useful. However, if there are paired conjugate imaginary eigenvalues, then some (m, m-1) entry will not be eliminated by using Francis’s Implicit Iteration with only one shift, therefore, the following code can be unconvergent in this case. To deal with general real matrix, we can use the double shift strategy in next subsection.

iqr1S.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
% IQR1S. Implicitly QR Iteration Algorithm with single shift
% Input:    matrix A
% Output:   matrix M
% s.t.  By using Francis's Implicit Iteration, M(n,n) will become extremely  
% close to an eigenvalue of A, and M(n,n-1) is nearly 0 as iteration continues.
% This approach is applied repeatly to compute all eigenvalues of A.
% Remark. Only useful for matrix with real eigenvalues. If there are paired
% conjugate imaginary eigenvalues, then some (m, m-1) entry will not be 
% eliminated by using Francis's Implicit Iteration with only one shift.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P359.
% Test example: V=[2,2,3;2,3,4;3,4,5];S=diag([1,2,3]);A=V*S*inv(V);
function M = iqr1S(A)
    [n, ~] = size(A);

    %--------- Reduction to Upper Henssenberg form -----------------%
    A = upperHenssenberg(A);    % call upperHenssenberg.m
    A = triu(A, -1);            % reduce A to upper Henssenberg form
    %---------------------------------------------------------------%

    M = A;          % initalization
    maxN = 100;     % max iteration number

    % Reduce (-1)-subdiagonal entries from (n, n-1) to (2, 1)
    for m = n:-1:2
        % using Implicitly QR Iteration to reduce M(m, m-1) to 0
        % in next loop only need to consider a smaller (n-1)x(n-1) martix
        if abs(M(m, m-1)) < 1e-16 * (abs(M(m, m)) + abs(M(m-1, m)))
            M(m, m-1) = 0;      % check whether *proper* upper Henssenberg
        else
            M(1:m, 1:m) = deflation(M(1:m, 1:m), maxN);
        end
    end
end

% DEFLATION. Eliminate A(n, n-1) to make A(n,n) approach an
% eigenvalue of A by using Francis's Implicit Iteration.
% Input:    symmetric matrix A, max iteration number maxN
% Output:   matrix M
function M = deflation(A, maxN)
    [n, ~] = size(A);
    M = A;
    t = 1;                                      % counter
    while t <= maxN && abs(M(n, n-1)) > 1e-16   % stop if M(n,n-1) is small enough
        M = FrancisIteration(M);       % Francis's Iteration
        t = t+1;
    end
end

% FRANCISITERATION. Perform a Francis's Iteration.
% Input:    symmetric matrix A
% Output:   matrix M
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P362.
function A = FrancisIteration(A)
    [n, ~] = size(A);

    % assign the single shift to be Wilkinsion shift
    shift = wilkinson(A);
    % can also using Rayleigh quotient shift
    % shift = M(n, n);

    %----------- Construct Givens rotator -----------%
    cs = A(1, 1) - shift;
    sn = A(2, 1);
    % normalization
    r = norm([cs, sn]);
    cs = cs / r;
    sn = sn / r;
    q0 = [cs, -sn; sn, cs];     % Givens rotator
    %------------------------------------------------%
    
    %----------- Introduce a bulge ------------------%
    A(1:2, :) = q0'*A(1:2, :);  % left multiplication
    A(:, 1:2) = A(:, 1:2)*q0;   % right multiplication
    %------------------------------------------------%

    %================= Bulge-Chasing ================%
    for i = 1:n-2
        %--------- Construct Givens rotator ---------%
        cs = A(i+1, i);
        sn = A(i+2, i);
        % normalization
        r = norm([cs, sn]);
        cs = cs / r;
        sn = sn / r;
        qi = [cs, -sn; sn, cs];% Givens rotator
        %--------------------------------------------%

        %---------- eliminate the old bulge----------%
        A(i+1, i) = r;
        A(i+2, i) = 0;
        %--------------------------------------------%

        %--introduce a new bulge along (-1)-subdiag--% Similarity Transform
        A(i+1:i+2, i+1:n) = qi' * A(i+1:i+2, i+1:n); % Left multiplication
        A(:, i+1:i+2) = A(:, i+1:i+2) * qi;          % Right multiplication
        %--------------------------------------------%
    end
    %================================================%
end

% WILKINSON. Compute Wilkinson shift with respect to matrix A
% Input:    matrix A
% Output:   shift
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P361.
function shift = wilkinson(A)
    [n, ~] = size(A);
    htr = .5*(A(n-1, n-1) + A(n, n));   % half a 2x2 trace
    dscr = sqrt((.5*(A(n-1,n-1)-A(n,n)))^2 + A(n,n-1)^2);
    % discriminant
    if htr < 0
        dscr = -dscr;                   % to avoid cancellation
    end 
    root1 = htr + dscr;                 % quadratic formula
    if root1 == 0                       % almost never happens
        root2 = 0;
    else                                % almost always happens
        det = A(n-1,n-1)*A(n,n) - A(n,n-1)^2;
        % 2x2 determinant = product of roots
        root2 = det/root1;
    end
    if abs(A(n,n)-root1) < abs(A(n,n)-root2)
        shift = root1;
    else
        shift = root2;
    end
end

Double shift for General Real Matrix

iqr2S.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
% IQR2S. Implicitly QR Iteration Algorithm with double shift
% Input:    matrix A
% Output:   matrix M
% s.t.  The eigenvalues of A are approxed by the eigenvalue(s) of the 2x2
% or 1x1 submatrices in diag of M.
% Remark. To be improved: this code might break down if the matrix deduces
% to an *improper* upper Henssenberg matrix. 
% Unforturnately the convergence is not guaranteed, which is rare forturnately
% and can usually be avoided by an exceptional shift strategy: whenever the
% iterations seem not to be converging, take one step with random shifts to
% shake things up.
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P366.
% Test example:  A = [0 0 0 1; 1 0 0 0; 0 1 0 0; 0 0 1 0]
function M = iqr2S(A)
    [n, ~] = size(A);

    %--------- Reduction to Upper Henssenberg form -----------------%
    A = upperHenssenberg(A);    % call upperHenssenberg.m
    A = triu(A, -1);            % reduce A to upper Henssenberg form
    %---------------------------------------------------------------%

    M = A;          % initalization
    maxN = 100;     % max iteration number
    flag = 0;       % normal case

    m = n;
    % Reduce (-1)-subdiagonal entries from (n, n-1) to (2, 1)
    while m >= 3
        % using Implicitly QR Iteration to reduce M(m, m-1) to 0
        % in next loop only need to consider a smaller (n-1)x(n-1) martix
        if abs(M(m, m-1)) < 1e-16 * (abs(M(m, m)) + abs(M(m-1, m)))
            M(m, m-1) = 0;      % check whether *proper* upper Henssenberg
            d = 1;
        else
            [M(1:m, 1:m), d, t] = deflation(M(1:m, 1:m), maxN, flag);
            if t == maxN + 1
                flag = 1;       % Convergence Failure, set an exception flag
                % try another special deflation with one time random shifts
                [M(1:m, 1:m), d, ~] = deflation(M(1:m, 1:m), maxN, flag);
            end
            flag = 0;
        end
        m = m - d;
    end
end

% DEFLATION. Eliminate A(n, n-1) to make A(n,n) approach an
% eigenvalue of A by using Francis's Implicit Iteration.
% Input:    symmetric matrix A, max iteration number maxN
% Output:   matrix M
function [M, d, t] = deflation(A, maxN, flag)
    [n, ~] = size(A);
    M = A;
    t = 1;                                      % counter
    
    d = 2;
    while t <= maxN && abs(M(n-1, n-2)) > 1e-16 % stop if M(n-1,n-2) is small enough

        if abs(M(n, n-1)) < 1e-16             
            d = 1;                              % stop if M(n,n-1) is small enough
            break
        end

        M = FrancisIteration(M, flag);          % Francis's Iteration

        if flag ~= 0
            flag = 0;
        end

        t = t+1;
    end
end

% FRANCISITERATION. Perform a Francis's Iteration.
% Input:    symmetric matrix A
% Output:   matrix M
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P362.
function A = FrancisIteration(A, flag)
    [n, ~] = size(A);

    if flag == 0        % Normal case
        % assign the single shift to be Wilkinsion shift
        rho = RayleighQuotient(A);
        % can also using Rayleigh quotient shift
        % shift = M(n, n);
    else                % Convergence Failure, try a random shift
        rho = randn(2,1);
    end

    % First column
    p = [(A(1,1)-rho(1))*(A(1,1)-rho(2))+A(1,2)*A(2,1);
        A(2,1)*((A(1,1)+A(2,2)) - rho(1) - rho(2));
        A(3,2) * A(2,1)];
    p = real(p);                                     % p is real if A is a real matrix

    %----------- Householder Reflector --------------%
    [~, gamma, u] = reflector(p);
    % Q = (I - gamma * u*u') is a reflector for which
    % Qx = [-tau, 0, ..., 0]'
    %------------------------------------------------%
    
    %----------- Introduce a bulge ------------------%
    A(1:3, :) = A(1:3, :) - gamma*u* (u'*A(1:3, :));
    A(:, 1:3) = A(:, 1:3) - gamma* (A(:, 1:3)*u)*u';
    %------------------------------------------------%

    %================= Bulge-Chasing ================%
    for i = 1:n-3
        %--------- Householder Reflector ---------%
        [tau, gamma, u] = reflector(A(i+1:i+3, i));
        %--------------------------------------------%

        %---------- eliminate the old bulge----------%
        A(i+1, i) = -tau;
        A(i+2, i) = 0;
        A(i+3, i) = 0;
        %--------------------------------------------%

        %--introduce a new bulge along (-1)-subdiag--% Similarity Transform
        A(i+1:i+3, i+1:n) = A(i+1:i+3, i+1:n) - gamma*u* (u'*A(i+1:i+3, i+1:n));% Left multiplication
        A(:, i+1:i+3) = A(:, i+1:i+3) - gamma* (A(:, i+1:i+3)*u)*u';            % Right multiplication
        %--------------------------------------------%
    end

    % Last reflector only need to deal with the last *two* rows and columns
    [tau, gamma, u] = reflector(A(n-1:n, n-2));
    A(n-1, n-2) = -tau;
    A(n, n-2) = 0;
    A(n-1:n, n-1:n) = A(n-1:n, n-1:n) - gamma*u* (u'*A(n-1:n, n-1:n));          % Left multiplication
    A(:,n-1:n) = A(:, n-1:n) - gamma* (A(:, n-1:n)*u)*u';                       % Right multiplication
    %================================================%
end

% RAYLEIGHQUOTIENT. Compute generalized Rayleigh Quotient shift
% Input:    matrix A
% Output:   vector shifts = [rho1; rho2]
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 5.6, P368.
function shifts = RayleighQuotient(A)
    [n, ~] = size(A);
    htr = .5*(A(n-1, n-1) + A(n, n));   % half a 2x2 trace
    dscr = sqrt((.5*(A(n-1,n-1)-A(n,n)))^2 + A(n,n-1)*A(n-1,n));
    % discriminant
    if htr < 0
        dscr = -dscr;                   % to avoid cancellation
    end 
    root1 = htr + dscr;                 % quadratic formula
    if root1 == 0                       % almost never happens
        root2 = 0;
    else                                % almost always happens
        det = A(n-1,n-1)*A(n,n) - A(n,n-1)*A(n-1,n);
        % 2x2 determinant = product of roots
        root2 = det/root1;
    end
    shifts = [root1; root2];
end

% REFLECTOR. construct Householder Transform reflector
% Input:    vector x with length n
% Output:   parameters tau, gamma and vector u
% s.t.  Q = (I - gamma * u*u') is a reflector for which
%       Qx = [-tau, 0, ..., 0]'
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P201.
function [tau, gamma, u] = reflector(x)
    beta = norm(x, "inf");              % check whether x=0
    if beta < 2^(-52)                   % if x=0, no need to reflect x
        gamma = 0;                      % gamma is set to be zero, giving Q=I
        tau = 0;
    else
        % ------------ x!=0 ------------%
        x = x / beta;                   % normalization

        tau = norm(x);                  % Q is unitary, norm(Qx) = tau = norm(x)

        if x(1) < 0                     % check if sign of x1 is same with tau
            tau = - tau;                % reverse sign of tau to aovid cancellation
        end
        x(1) = tau + x(1);              % same sign -> without cancellation

        gamma = x(1) / tau;             % recall gamma = (tau + x1) / tau

        x(2:end) = x(2:end) / x(1);     % [u2, ..., un] = [x2, ..., xn] / (tau + x1)
        x(1) = 1;                       % u1 = 1

        tau = tau * beta;               % re-scaling
    end
    u = x;                              % u is stored in x
end
  • Title: QR decomposition & QR algorithm
  • Author: Gypsophila
  • Created at : 2025-03-28 22:18:40
  • Updated at : 2025-03-30 17:35:33
  • Link: https://chenx.space/2025/03/28/QR/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments
On this page
QR decomposition & QR algorithm
  1. QR decomposition
    1. Classical Gram-Schmidt Orthogonalization
    2. Modified Gram-Schmidt Orthogonalization
    3. Householder transformation
    4. Givens rotation
  2. QR algorithm
    1. QR Iteration
    2. Francis’s Implicitly Shifted QR Algorithm
      1. Single shift for Symmetric Matrix
      2. Double shift for General Real Matrix