Matlab Code-Include Iteration to QR Algorithm Gram-Schmidt - The Iterations of A will converge to Eigenvalues

Question

Still need to add the iteration to the Matlab Code of the QR Algorithm using Gram-Schmidt to iterate until convergence as follows:

I am having trouble completing the code to be able to iterate the algorithm that displays $Q$. The eigenvalues will become clear in the diagonal after so many iterations of the formula noted below which will give the next $A$. So far I have that the code will calculate $Q_0$ from $A_0=A$, which will be used in the following:

The iteration is such that $$A_{m+1}=Q_{m}^T*A_{m}*Q_{m}$$

Then the first iteration will give me $A_1$, which I will now use the algorithm to obtain the next $Q_1$, but getting stuck on the iteration part of the code. I need it to take each new Q, calculate in the above formula and display the next A. But, for now my code only calculate $Q$ from my initial $A$.

I need for the code to display $A_1, A_2,...$ until the eigenvalues become apparant on the diagonal. The eigenvalues will display nicely with this process. (Cannot us the eig function)

Matlab code so far:

function [Q R]=QR_GramSchmidt(A) % QR decomposition by Gram-Schmidt method

A=[3,1,0;1,4,2;0,2,1]; % I am testing

n=3; % Testing

[n n]=size(A);

Q=zeros(n);

R=zeros(n);

R(1,1)=norm(A(:,1));

Q(:,1)=A(:,1)/R(1,1);

for j=2:n

    Q(:,j)=A(:,j);

    for i=1:j-1

        Q(:,j)=Q(:,j)-A(:,j)'*Q(:,i)*Q(:,i);

        R(i,j)=A(:,j)'*Q(:,i);

    end

    R(j,j)=norm(Q(:,j));

    Q(:,j)=Q(:,j)/norm(Q(:,j));

end   

end

See below picture of $A_8$ producing $Q_8$, then $Q_8$ is used in the iteration formula to produce $A_9$ which has the eigenvalues in the diagonal (You could see it converging previous to $A_9$).

[Result of Manually entering each new Q to obtain the next A-Eigenvalues are in the diagonal as it converged]1

Answer 1

Your code is correct. Run it. Then do the following multiplication $QR$, it will turn out to be $A$. You will also get that $Q$ is the orthogonal matrix $Q^TQ = I$ and $R$ is upper triangular. You can print $A_i$ (values of $A$ at iteration $i$) as follows

function [Q R]=QR_GramSchmidt(A)
[n n]=size(A);
Q=zeros(n);
    R=zeros(n);
    R(1,1)=norm(A(:,1));
    Q(:,1)=A(:,1)/R(1,1);

for j=2:n
    Q(:,j)=A(:,j);
    for i=1:j-1
    Q(:,j)=Q(:,j)-A(:,j)'*Q(:,i)*Q(:,i);

    R(i,j)=A(:,j)'*Q(:,i);

    fprintf(['At iteration i = ',num2str(i),' QR is equal to ']) % print Ai
    Q*R
end
R(j,j)=norm(Q(:,j));
Q(:,j)=Q(:,j)/norm(Q(:,j));
fprintf(['At iteration j = ',num2str(j),' QR is equal to '])
    Q*R
end   
end

Output will be as follows

>> A=[3,1,0;1,4,2;0,2,3];
>> [Q,R] = QR_GramSchmidt(A)

At iteration i = 1 QR is equal to ans =

3.0000    2.1000         0
1.0000    0.7000         0
     0         0         0

At iteration j = 2 QR is equal to ans =

 3     1     0
 1     4     0
 0     2     0

At iteration i = 1 QR is equal to ans =

3.0000    1.0000    0.6000
1.0000    4.0000    0.2000
     0    2.0000         0

At iteration i = 2 QR is equal to ans =

3.0000    1.0000   -0.2609
1.0000    4.0000    2.7826
     0    2.0000    1.5652

At iteration j = 3 QR is equal to ans =

 3     1     0
 1     4     2
 0     2     3

Q =

0.9487   -0.2741    0.1576
0.3162    0.8224   -0.4729
     0    0.4984    0.8669

R =

3.1623    2.2136    0.6325
     0    4.0125    3.1402
     0         0    1.6550

Answer 2

Just implementing some code here..

function [T,Q,R] =basicqr(A,niter)
%
%
%
T = A;
for i=1:niter

    [Q,R] =qr(T);
    T=R*Q;
end

end


A=[3,1,0;1,4,2;0,2,1];
n=1000;

[T,Q,R] = basicqr(A,n);
[V,D] = eigs(A)
T

V =

    0.1518    0.9201   -0.3610
   -0.4658   -0.2556   -0.8472
    0.8718   -0.2968   -0.3898


D =

   -0.0687         0         0
         0    2.7222         0
         0         0    5.3465


T =

    5.3465    0.0000   -0.0000
   -0.0000    2.7222    0.0000
         0         0   -0.0687

the eigenvalues are right...they're just switched...you see.

You need to use a deflation technique..

function d = qralgH(B)
% QRALGH QR algorithm with shifts and deflation  
%        after applying Hessenberg reduction.

% Iterate over eigenvalues

A = hess1(B);

for n = length(A):-1:2
  % QR iteration
  while sum( abs(A(n,1:n-1)) ) > eps
    s = A(n,n);
    [Q,R] = qr(A-s*eye(n));
    A = R*Q + s*eye(n);
  end 
  % Deflation
  d(n) = A(n,n);
  A = A(1:n-1,1:n-1);
end
% Last remaining eigenvalue
d(1) = A(1,1);

% d = flipud(sort(d));

function d = qralg(A)
% QR algorithm with shifts and deflation.

% Iterate over eigenvalues
for n = length(A):-1:2
  % QR iteration
  while sum( abs(A(n,1:n-1)) ) > eps
    s = A(n,n);
    [Q,R] = qr(A-s*eye(n));
    A = R*Q + s*eye(n);
  end 
  % Deflation
  d(n) = A(n,n);
  A = A(1:n-1,1:n-1);
end
% Last remaining eigenvalue
d(1) = A(1,1);
end
end

function H = hess1(A)
% HESS1  Reduce a matrix to Hessenberg form.
% (demo only--not efficient)
% Input:
%   A      n-by-n matrix
% Output:
%   H     m-by-n upper Hessenberg form

[n,n] = size(A);
H = A;
for k = 1:n-1
  Q = eye(n);
  z = H(k+1:n,k);
  v = [ -sign(z(1))*norm(z)-z(1); -z(2:end) ];
  P = eye(n-k) - 2*(v*v')/(v'*v);   % HH reflection
  Q(k+1:n,k+1:n) = P;
  H = Q*H*Q';
end
H = triu(H,-1);

end

A=[3,1,0;1,4,2;0,2,1];
n=1000;

[T,Q,R] = basicqr(A,n);
[V,D] = eigs(A)
T
d = qralgH(A)

T =

    5.3465    0.0000   -0.0000
   -0.0000    2.7222    0.0000
         0         0   -0.0687


d =

    5.3465    2.7222   -0.0687