Calculating the eigen vectors without using eig function but using a data set, the transpose of the data set and a unit vector

Hello,

I am trying to calculate the eigen vectors using a looped equation without using the eig function but using a dataset, the transpose, and a random unit vector with the magnitude of 1.

here is what I have so far:

clc

clear

A = importdata('mariana_depth.csv');

B=transpose (A);

C=B*A;

u=rand(1440,1); u=u/norm(u);

for i=1:10

u_n=u;

u=C*u

v(:,i)=u;

u_n(:,i) =u.*v(:,i)/norm(u.*v(:,i));

diffu = norm(u_n-u);

i

diffu

end

the problem I am having is that it is counting up but not down so I dont believe it it calculating the proper values of my eigen values.

[Zachary David]

this is the equation i am trying to use

u= the unit vector with a magnitude of 1

u_n+1= transposeA * A* u / || transposeA * A * u||

repeating ||u_n+1 - u|| until it is very small

Answer 1

[Sam Chak]

Hi @Zachary David

Can consider an alternative approach by finding the characteristic polynomial of a matrix using poly()

https://www.mathworks.com/help/matlab/ref/poly.html

and then use some root-finding algorithms or roots() to find them.

https://www.mathworks.com/help/matlab/ref/roots.html

For example:

A = [0 1; -2 -3]
p = poly(A)
r = roots(p)

A =

0 1

-2 -3

p =

1 3 2

r =

-2

-1

[Zachary David]

Hello Sam,

unfortunately I am required to use the equation.

[Zachary David]

If i try it a different way it stops after one iteration. This is just a test code

A=[1 3 4; 5 6 7; 2 9 8];

u=rand(3,1); u= u/norm(u);

B=transpose(A);

C=B*A;

B;

C;

i=0

diff=0;

while diff<0.0001

i=i+1;

u_n=u;

u_n=C*u/norm(C*u)

diff= norm(u_n-u)

i

diff

end

[Zachary David]

Sorry i know that was probably confusing disregard my last post. Its late

[Sam Chak]

No worries, @Zachary David. Take a good rest tonight.

There must be something stopping the iteration. Need to study your proposed algorithm and to compare with other algorithms:

• Power iteration
• Inverse iteration
• Rayleigh quotient iteration
• Arnoldi iteration — based on Krylov subspaces
• Lanczos algorithm — Arnoldi, specialized for positive-definite matrices
• Block Lanczos algorithm — for when matrix is over a finite field
• QR algorithm
• Jacobi eigenvalue algorithm — select a small submatrix which can be diagonalized exactly, and repeat
• Jacobi rotation — the building block, almost a Givens rotation
• Jacobi method for complex Hermitian matrices
• Divide-and-conquer eigenvalue algorithm
• Folded spectrum method
• LOBPCG — Locally Optimal Block Preconditioned Conjugate Gradient Method
• Eigenvalue perturbation — stability of eigenvalues under perturbations of the matrix

https://en.wikipedia.org/wiki/Eigenvalue_algorithm#Iterative_algorithms

[Zachary David]

essentially the original code will iterate through but the eigen values are becoming larger not smaller

[Zachary David]

ok thanks ill try and look into it more.

[Sam Chak]

Think your method is related to the Power method. But the script below can only calculate the real largest eigenvalue.

function Demo_Eigen

close all

clc

A = [1 3 4; 5 6 7; 2 9 8]

% A = [0 1; -1 -1]

% A = inv(A) % Inverse Iteration Method to find smallest eigenvalue

n = size(A, 1);

x1 = ones(n, 1);

epsilon = 1e-6;

kmax = 100;

x(:, 1) = x1./norm(x1, 2); % initial eigenvector

x(:, 2) = A*x(:, 1);

lambda(2) = x(:, 1).'*x(:, 2); % initial eigenvalue

x(:, 2) = x(:, 2)./norm(x(:, 2), 2);

for k = 2:kmax

x(:, k+1) = A*x(:, k);

lambda(k + 1) = x(:, k).'*x(:, k + 1);

x(:, k + 1) = x(:, k + 1)./norm(x(:, k + 1), 2);

if abs(lambda(k + 1) - lambda(k)) < epsilon

break

end

end

Eigen_val = lambda(end)

Eigen_vec = x(:,end)

end

Answer 2

[Bruno Luong]

Your method give you the largest singular vector not eigen vector:

A=rand(5);

C=A'*A;
u=randn(size(A,2),1);
u=u/norm(u);
for i=1:100
    u_n=C*u;
    u_n=u_n/norm(u_n);
    diffu = norm(u_n-u);
    if diffu < 1e-10
        fprintf('converge afeter %d iterations\n', i)
        break
    end
    u=u_n;
end

converge afeter 9 iterations

u

u = 5×1

-0.5840 -0.5132 -0.2270 -0.2861 -0.5121

[U,S,V]=svd(A);
V(:,1)

ans = 5×1

-0.5840 -0.5132 -0.2270 -0.2861 -0.5121

[Bruno Luong]

If you want a smallest singular vector you need to use backslash instead of time

A = rand(5)

A = 5×5

0.7672 0.0288 0.8776 0.3047 0.9679 0.8588 0.8712 0.5263 0.7778 0.4880 0.5808 0.5344 0.5108 0.4657 0.6664 0.8760 0.8231 0.0007 0.5201 0.3654 0.2871 0.6496 0.1216 0.0663 0.6528

C=A'*A;
u=randn(size(A,2),1);
u=u/norm(u);
for i=1:100
    u_n=C\u;
    u_n=u_n/norm(u_n);
    diffu = norm(u_n-u);
    if diffu < 1e-10
        fprintf('converge afeter %d iterations\n', i)
        break
    end
    u=u_n;
end

converge afeter 7 iterations

u

u = 5×1

0.2794 0.3037 0.4550 -0.6649 -0.4249

[U,S,V]=svd(A);
V(:,end)

ans = 5×1

0.2794 0.3037 0.4550 -0.6649 -0.4249