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

gram-schmidtMATLABmatricesmatrix decompositionnumerical linear algebra

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

Best Answer

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

Related Question