0
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ I am getting stuck on how to code the iteration so that I may run it so many times, thus showing A1, A2, .. $\endgroup$ – user565684 Jul 19 '18 at 17:13
  • $\begingroup$ OK. I have edited. $\endgroup$ – Ahmad Bazzi Jul 19 '18 at 17:25
  • $\begingroup$ 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 :) $\endgroup$ – user565684 Jul 19 '18 at 17:45
  • $\begingroup$ 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. $\endgroup$ – Ahmad Bazzi Jul 19 '18 at 17:48
  • $\begingroup$ So sorry my test matrix last digit is 1 not 3 like I noted $\endgroup$ – user565684 Jul 19 '18 at 17:51
0
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – user565684 Jul 19 '18 at 19:59
  • $\begingroup$ Im currently on the road and accidently put it in there. Typically you use hessenberg reduction. It isnt gram schmidt exactly. Halves the computations. $\endgroup$ – user3417 Jul 19 '18 at 20:16
  • $\begingroup$ 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 $\endgroup$ – user3417 Jul 19 '18 at 20:18
  • $\begingroup$ 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 :) $\endgroup$ – user565684 Jul 19 '18 at 20:28
  • $\begingroup$ Right but theyre in the wrong place. $\endgroup$ – user3417 Jul 19 '18 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.