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

$\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$

2 Answers

$\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$ 14 $\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$ 13

Your Answer

Sign up or log in

Sign up using Google Sign up using Facebook Sign up using Email and Password Submit

Post as a guest

Post Your Answer Discard

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