Power method to determine largest eigenvalue and corresponding eigenvector had some wrong

I used MATLAB eig function to check answer, the answer should be 3.3876 for largest eigenvalue and [-0.371748 0.601501 -0.601501 0.371748] for the corresponding eigenvector. I have no idea why it is wrong.

Thanks for helping.

Code here:

% Matlab eig function
clc;clear
A = [1.7696        -1           0          0
    -1         1.7696      -1          0
    0         -1         1.7696      -1
    0          0           -1      1.7696 ];
[eigvec,eigval] = eig(A);
for i=1:4
    fprintf("Eigenvalue %f, corresponding eigenvector:[%f %f %f %f]\n",eigval(i,i),eigvec(:,i))
end

Eigenvalue 0.151566, corresponding eigenvector:[0.371748 0.601501 0.601501 0.371748]
Eigenvalue 1.151566, corresponding eigenvector:[-0.601501 -0.371748 0.371748 0.601501]
Eigenvalue 2.387634, corresponding eigenvector:[-0.601501 0.371748 0.371748 -0.601501]
Eigenvalue 3.387634, corresponding eigenvector:[-0.371748 0.601501 -0.601501 0.371748]

fprintf("\n")
%power method
[powval,powvec] = powereig(A,10^(-6));
fprintf("The largest eigenvalue and its corresponding eigenvector by power method:\n")

The largest eigenvalue and its corresponding eigenvector by power method:

fprintf("Eigenvalue:%f\n",powval)

Eigenvalue:2.387634

fprintf("Corresponding vector:[%f %f %f %f]\n",powvec)

Corresponding vector:[1.000000 -0.618034 -0.618034 1.000000]

function [eval, evect] = powereig(A,es)
n=length(A);
evect=ones(n,1);eval=1;iter=0;ea=100; %initialize
tol = es;
while(1)
  evalold=eval;           %save old eigenvalue value
  evectold = evect;       %save old eigenvalue vector
  evect=A*evect;          %determine eigenvector as [A]*{x)
  eval=max(abs(evect));   %determine new eigenvalue
  evect=evect./eval;      %normalize eigenvector to eigenvalue
  iter=iter+1;
  if eval~=0, ea = abs((eval-evalold)/eval)*100; end
  if ea<=es  || norm(evect - evectold,2) <= tol
      break
  end
end
end

Answer 1

[John D'Errico]

I would claim you have two problems here. One is a pretty common mistake, the other more subtle.

Do you understand that eig normalizes the eigenvectors, to have a unit 2-norm?

V = [1.000000 -0.618034 -0.618034 1.000000];
V./norm(V)

ans = 1×4

0.6015 -0.3717 -0.3717 0.6015

And indeed, that is the vector corresponding to the second largest eigenvalue. So the code worked. You just needed to normalize the vector.

Still, it looks like your power method did not find the largest eigenvalue however. But that is not immediately relevant, since your code did work, in a sense. However, why did your method not find the largest eigenvalue/eigenvector pair?

The answer is a subtle one. In fact, if you look at the actual eigenvector, it turns out to be orthogonal to the starting vector you used. You started with a vector of all ones. That is a really bad idea in this case, because made up examples like this will often cause problems.

Instead, you would have been far better off using a randomly chosen vector of the correct length. That makes the probability equal to zero of starting with a vector orthogonal to one of your eigenvectors.

[Yen Hung-lu]

Oh, I didn't notice that eig normalizes the eigenvectors.

I change evect=ones(n,1) to evect=rand(n,1), then answer match eig function.

Thanks for helping !

I learn a lot from this question.