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

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$).

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


• I am getting stuck on how to code the iteration so that I may run it so many times, thus showing A1, A2, .. – user565684 Jul 19 '18 at 17:13
• OK. I have edited. – Ahmad Bazzi Jul 19 '18 at 17:25
• Something is not quite right. For instance the matrix I am using will need to be iterated I believe it was 7 times when I manually tested it to obtain the eigenvalues. Where A_7=(Q_6)^T * A_6 * Q_6 resulting in eigenvalues of 5.346, 2.722, -0.068. So, I am in need of obtaining a new Q to obtain the next A, not the breakdown of Q and R. So QR method will calculate each new Q to obtain the next A and then display each A, ie: A1, A2,...., A7 in this case. Sorry - I hope that I explained well :) – user565684 Jul 19 '18 at 17:45
• But the eigenvalues of A are 1.2087, 3 , 5.7913 .. In the end, you get QR = A, where Q is orthogonal and R is upper triangular. Hence, it is correct. – Ahmad Bazzi Jul 19 '18 at 17:48
• So sorry my test matrix last digit is 1 not 3 like I noted – user565684 Jul 19 '18 at 17:51

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

• Can you determine based on my code what I need to add to have the algorithm take the Q it produces and plug it into the formula I have noted to obtain the next A value, continue iterating this way until convergence of the eigenvalues? I have run this manually, and have found convergence by A_7 – user565684 Jul 19 '18 at 19:59
• Im currently on the road and accidently put it in there. Typically you use hessenberg reduction. It isnt gram schmidt exactly. Halves the computations. – user3417 Jul 19 '18 at 20:16
• Youre missing multiples things. Once im home i can post more. But this is classic gram schmidt. It typically goes to hessenberg first. Then there is deflation. Sometimes you use block techniques – user3417 Jul 19 '18 at 20:18
• Thank you. It will converge with adding the basic iteration formula to the code, taking the current Q Matlab produces, obtaining the next A. Do this 7 times and it works :) – user565684 Jul 19 '18 at 20:28
• Right but theyre in the wrong place. – user3417 Jul 19 '18 at 20:32