Insignificant modification of code leads to NaN

I am an utter beginner at programming and I wrote a small script to inplement the Power Method for finding eigenvalues of a matrix. The iterates seem to converge fine, to the vector [0.7071, 0.7071]. However, iterate 511 (column 511 in the matrix R) suddenly becomes [NaN; 0], and then everything else is [NaN, NaN]. Here is the code:

%Starting point of the iteration
x = [1; 4];

%The matrix whose eigenvalue we approximate
A = [-1, 5;
    1, 3];
num_iterates = 1000;
for i = 1:num_iterates
    x = A*x;
    
    %Store the iterates into the matrix R
    R(:,i) = x/norm(x); %%%% Problem Line %%%%
end

%Show the last iterate
v = R(:, num_iterates)
A*v

Now, I know the error is with the line labeled "Problem Line". If I change the code to the following, then everything is fine. Could you please help me understand why the second code is good, while the first is not?

%Starting point of the iteration
x = [1; 4];

%The matrix whose eigenvalue we approximate
A = [-1, 5;
    1, 3];
num_iterates = 1000;
for i = 1:num_iterates
    x = A*x;
    x = x/norm(x)
    %Store the iterates into the matrix R
    R(:,i) = x;
end

%Show the last iterate
v = R(:, num_iterates)
A*v

Answer 1

[James Tursa]

Think about what is happening with this loop:

for i = 1:num_iterates
    x = A*x;
    R(:,i) = x/norm(x); %%%% Problem Line %%%%
end

Suppose x starts out as the unit eigenvector e with associated max eigenvalue lambda. Then

1st iteration: x = lambda*e

2nd iteration: x = lambda^2*e

3rd iteration: x = lambda^3*e

etc.

At each iteration x grows in magnitude by factor lambda. Eventually it will overflow to infinity or underflow to 0 (depending on value of lambda) and lead to the NaN issue.

In your other loop, you are storing unit vectors for x at each iteration, which avoids this problem.

In your particular problem, the max eigenvalue is 4. And if we look at the following:

>> norm([1;4])*4^510
ans =
  4.6325e+307
>> norm([1;4])*4^511
ans =
   Inf

You can see that we would expect an overflow at about iteration 511, and that's what you got.

[Jack]

Thank you very much! Yes, I knew we have to normalize x at each iteration. My brain just didn't register that we are not updating x correctly in the erroneous code.