Algorithms about eigenvalues

  1. Implement Matlab functions which can compute the corresponding eigenpair(s) via:

(a) Power iteration.

(b) Inverse iteration.

(c) Rayleign Quotient iteration.

(d) Simultaneous iteration.

(e) QR algorithm.

  1. Compare the differences and performances of these algorithms on a large-scaled problem.

Experimental purpose and requirement

Implement 5 funcions about eigonvalue and eigonvector.

Main equipment:

Hardware: MacBook Pro 2020, Catalina 11.15.7, 2 GHz Quad-Core Intel Core i5, 16 GB 3733 MHz LPDDR4X

MATLAB R2021b

Experimental procedure:

Write 5 function as required, use tic toc to record the time overhead of each calculating method.

Experimental result and analysis:

clc
sz = 50;
epsi = 1e-10;
maxx = 10000;
maxn = 5000;
A = randi([1,10],[sz,sz]);
A = A*A.';

% eig(A)
tic
for i = 1:maxn
PowerIteration(A,sz);
end
toc
tic
for i = 1:maxn
InverseIteration(A,sz);
end
toc
tic
for i = 1:maxn
RQIteration(A,sz);
end
toc
tic
for i = 1:maxn
SimulIteration(A);
end
toc
tic
for i = 1:maxn
QRALgorithm(A);
end
toc


% Power iteration
function [lamb] = PowerIteration(A,sz)
    epsi = 1e-10;
    maxx = 10000;
    v = randi([1,10],[sz,1]);
    v = v/norm(v);

    for i = 1:maxx
        w = A*v;
        v = w/norm(w);
        if mod(i,5)==0
            lamb = v.' * A * v; 
            if norm(A*v - lamb.*v) < epsi
%                 disp(['loop ',num2str(i),' time(s)']);
                break;
            end
        end
    end
end


% Inverse iteration
function [lamb] = InverseIteration(A,sz)

    epsi = 1e-10;
    maxx = 10000;
    v = randi([1,10],[sz,1]);
    v = v/norm(v);
    mu = 100;
    XX = (A-mu.*eye([sz,sz]))^(-1);
    for i = 1:maxx
        w = XX*v;
        v = w/norm(w);
        if mod(i,5)==0
            lamb = v.' * A * v; 
            if norm(A*v - lamb.*v) < epsi
%                 disp(['loop ',num2str(i),' time(s)']);
                break;
            end
        end
    end
end


% Rayleigh Quotient Iteration 
function [lamb] = RQIteration(A,sz)
    epsi = 1e-10;
    maxx = 10000;
    v = randi([1,10],[sz,1]);
    lamb = v.' * A * v;
    for i = 1:maxx
        w = (A - lamb.* eye([sz,sz]))\v;
        v = w/norm(w);
        if mod(i,5)==0
            lamb = v.' * A * v; 
            if norm(A*v - lamb.*v) < epsi
%                 disp(['loop ',num2str(i),' time(s)']);
                break;
            end
        end
    end
end


% Simultaneous Iteration
function [Q] = SimulIteration(A)
    maxx = 1000;
    epsi = 1e-10;
    [Q,R] = qr(A);
    for i = 1:maxx
        Z = A*Q;
        [Q,R] = qr(Z);
        if mod(i,5)==0
            v = Q(:,1);
            lamb = v.' * A * v; 
            if norm(A*v - lamb.*v) < epsi
%                 disp(['loop ',num2str(i),' time(s)']);             
                break;
            end
        end
    end
end

% QR Algorithm
function [A] = QRALgorithm(A)
    maxx = 1000;
    epsi = 1e-10;
    for i = 1:maxx
        [Q,R] = qr(A);
        A = R*Q;
        if mod(i,5)==0
            v = Q(:,1);
            lamb = v.' * A * v; 
            if norm(A*v - lamb.*v) < epsi
%                 disp(['loop ',num2str(i),' time(s)']);             
                break;
            end
        end
    end
end

Elapsed time is 0.051586 seconds.
Elapsed time is 0.314361 seconds.
Elapsed time is 3.314072 seconds.
Elapsed time is 3.584557 seconds.
Elapsed time is 3.077736 seconds.
>> 
posted @ 2022-03-07 13:15  miyasaka  阅读(48)  评论(0)    收藏  举报