Is NaN a sign of error in MATLAB's Jacobi method?

[evinda]

Hi :) I have written a code in MATLAB to implement the Jacobi method.I calculated the error \left \| x_{k}-D \right \| (where D is the real solution) of the last iteration for different values of n .For n\geq 250 ,I get this as result: NaN.
Does this mean that I have done something wrong??

 

[Ackbach]

Yes, there's probably something wrong if that's supposed to be an error. NaN means "Not a Number". MATLAB responds with that if the result is infinite, or if there's a category error like you ask it to do something with strings that you normally only do with numbers. Can you post your code?

 

[evinda]

Yes, there's probably something wrong if that's supposed to be an error. NaN means "Not a Number". MATLAB responds with that if the result is infinite, or if there's a category error like you ask it to do something with strings that you normally only do with numbers. Can you post your code?

That is my code:

it=1;
x2=-inv(F)*W*x1+inv(F)*b
while ((it<=maxit)&&(norm(x2-x1)<e)&&(norm(b-A*x2))<e))
        error=norm(x-D);
        x1=x2;
        x2=-inv(F)*W*x1+inv(F)*b
        it=it+1;
end

I get Nan as result when I apply the Jacobi method at the Hilbert matrix,for n\geq 250.Why does this happen?

 

[Ackbach]

There are a lot of variables here that I don't know what they are. Could you please provide a complete cheat sheet for all your variables? Ideally, you could have all their definitions. I especially want to know the dimensions of all variables (whether matrix or vector). Thanks!

 

[evinda]

There are a lot of variables here that I don't know what they are. Could you please provide a complete cheat sheet for all your variables? Ideally, you could have all their definitions. I especially want to know the dimensions of all variables (whether matrix or vector). Thanks!

it=1;
x2=-inv(F)*W*x1+inv(F)*b
while ((it<=maxit)&&(norm(x2-x1)<e)&&(norm(b-A*x2))<e))
        error=norm(x-D);
        x1=x2;
        x2=-inv(F)*W*x1+inv(F)*b
        it=it+1;
end

F:the diagonal matrix of the Hilbert matrix
W:the Hilbert matrix - F
b:Ax=b
it:number of iteration
x2:approximation of solution
x1: previous approximation of solution
maxit:maximum number of iterations
e:a small number given from the user and
error:the error we are looking for

 

[Ackbach]

Ok, I have a few comments/thoughts:

1. You could probably speed this code up just a bit by setting Finv = inv(F), and referencing Finv instead of inv(F). As you are not changing F in the loop, it would be faster simply to recall the contents of memory somewhere than to evaluate the inverse of F multiple times in the loop.

2. Refering to your OP, is n the number of iterations? If so, could you please give me a plot of several errors on the y-axis versus the corresponding number of iterations on the x axis? I'm wondering if you might not be getting an underflow. That is, 250 iterations of the Jacobi method might make your error so small that you run up against the machine epsilon, and the computer cannot represent an error that small.

3. There might be a logical error in your program: you have the line
error=norm(x-D);
Should you be saying, instead, that
error=norm(x1-D)? Or error=norm(x2-D)? (Probably the latter, due to first-time-through issues.)
I mean, does the computer know what x is? It looks to me like you only know what x1 and x2 are (x isn't in your list). That indicates to me that that it's probably a typo in your code.

4. The while statement has unbalanced parentheses. Try enclosing the last test in its own parentheses.

See if one of these things, especially 2 or 3, doesn't help.

 

[evinda]

Ok, I have a few comments/thoughts:

1. You could probably speed this code up just a bit by setting Finv = inv(F), and referencing Finv instead of inv(F). As you are not changing F in the loop, it would be faster simply to recall the contents of memory somewhere than to evaluate the inverse of F multiple times in the loop.

2. Refering to your OP, is n the number of iterations? If so, could you please give me a plot of several errors on the y-axis versus the corresponding number of iterations on the x axis? I'm wondering if you might not be getting an underflow. That is, 250 iterations of the Jacobi method might make your error so small that you run up against the machine epsilon, and the computer cannot represent an error that small.

3. There might be a logical error in your program: you have the line
error=norm(x-D);
Should you be saying, instead, that
error=norm(x1-D)? Or error=norm(x2-D)? (Probably the latter, due to first-time-through issues.)
I mean, does the computer know what x is? It looks to me like you only know what x1 and x2 are (x isn't in your list). That indicates to me that that it's probably a typo in your code.

4. The while statement has unbalanced parentheses. Try enclosing the last test in its own parentheses.

See if one of these things, especially 2 or 3, doesn't help.

1)Ok..And...could I also write it like that: x2=-F\W*x1+F\b ??

2)I wrote this code in a function and with n I mean the maxit I give each time when I call the function...The underflow,has it to do with the Hilbert matrix??Because,when I ran my code,using other matrices,I didn't get Nan as result...

3)Oh sorry,I accidentally wrote it like that. At the code I ran ,I wrote it like that: error=norm(x2-D) ..

 

[Ackbach]

I mean something more like this:

it=1;
Finv=inv(F);
x2=-Finv*W*x1+Finv*b
while ((it<=maxit)&&(norm(x2-x1)<e)&&(norm(b-A*x2)<e))
        error=norm(x2-D);
        x1=x2;
        x2=-Finv*W*x1+Finv*b
        it=it+1;
end

 

[evinda]

I mean something more like this:

it=1;
Finv=inv(F);
x2=-Finv*W*x1+Finv*b
while ((it<=maxit)&&(norm(x2-x1)<e)&&(norm(b-A*x2)<e))
        error=norm(x2-D);
        x1=x2;
        x2=-Finv*W*x1+Finv*b
        it=it+1;
end

I wrote it like that right now..but I get the same result again(NaN)..Have I done maybe an other error??

 

[Ackbach]

There's something rather puzzling about your code: you check whether x1 is close to x2, which you wouldn't necessarily expect early in your iterations. So I can easily imagine a scenario where you wouldn't even enter the loop. However, presumably the loop is where you want the iterations to get closer and closer to the final answer. I don't have the time right this second to think about the proper loop condition, but think carefully about what comes right after the 'while' keyword.

 

[evinda]

There's something rather puzzling about your code: you check whether x1 is close to x2, which you wouldn't necessarily expect early in your iterations. So I can easily imagine a scenario where you wouldn't even enter the loop. However, presumably the loop is where you want the iterations to get closer and closer to the final answer. I don't have the time right this second to think about the proper loop condition, but think carefully about what comes right after the 'while' keyword.

I checked if,giving the Hilbert matrix with dimension 250,maxit=250,my code enters the loop,and it does...it makes 250 iterations...Have you an idea what I could have done wrong??

 

[evinda]

Also,the method does not converge...Has it maybe to do something with that??

 

[evinda]

But if I want to find the error of the last iteration of the Gauss-Seidel method,using the Hilbert Matrix,I get a number and not Nan as result...And,in this case the method does also not converge...

Why does this happen?

 

[evinda]

I know where the error is now...Using matlab,the determinant of the matrix equals to zero..But,what could I do,so that I have more precision of digits,to have a determinant \neq 0 ?

 

[Ackbach]

I know where the error is now...Using matlab,the determinant of the matrix equals to zero..But,what could I do,so that I have more precision of digits,to have a determinant \neq 0 ?

Having more precision of digits isn't necessarily going to help with the determinant being zero.

So... if the determinant is zero, then you most likely have infinite solutions with that particular matrix. Are you sure you are using the correct matrix? Can you double-check all the entries (I know there are a lot of them!)? If you're going to solve $Ax=b$, then having the right $A$ is definitely of the essence.

 

[evinda]

Having more precision of digits isn't necessarily going to help with the determinant being zero.

So... if the determinant is zero, then you most likely have infinite solutions with that particular matrix. Are you sure you are using the correct matrix? Can you double-check all the entries (I know there are a lot of them!)? If you're going to solve $Ax=b$, then having the right $A$ is definitely of the essence.

How could I check this? I use the command hilb(n)..

 

[Ackbach]

Hmm. Hilbert matrices are technically invertible, and hence have non-zero determinant. I wonder if there is a MATLAB setting you can change to use more precision. We are fast approaching the limits of my knowledge.

 

[evinda]

Hmm. Hilbert matrices are technically invertible, and hence have non-zero determinant. I wonder if there is a MATLAB setting you can change to use more precision. We are fast approaching the limits of my knowledge.

Where could I find such a command??

 

[evinda]

And...something else...is this right that both of the methods,the Jacobi and Gauss-Seidel method,do not converge for the Hilbert matrix?Because,no matter with dimension and which small number ε I give,the methods do not converge...

 

[Opalg]

I know where the error is now...Using matlab,the determinant of the matrix equals to zero..But,what could I do,so that I have more precision of digits,to have a determinant \neq 0 ?

The determinant of the Hilbert matrix is not zero, but it tends to zero exponentially fast as the size of the matrix increases. I know nothing about matlab, but I doubt whether it could possibly have enough precision to distinguish $\det H$ from zero for a $250\times 250$ Hilbert matrix. As is mentioned here, the Hilbert matrix is notoriously ill-conditioned for numerical computation.

 

[evinda]

The determinant of the Hilbert matrix is not zero, but it tends to zero exponentially fast as the size of the matrix increases. I know nothing about matlab, but I doubt whether it could possibly have enough precision to distinguish $\det H$ from zero for a $250\times 250$ Hilbert matrix. As is mentioned here, the Hilbert matrix is notoriously ill-conditioned for numerical computation.

I wrote the command det(H) and the result is that it equals to zero

 

[evinda]

I have searched in Google and found this:
Hilbert matrices are notoriously ill-conditioned: the computation of the determinant is subject to severe cancellation effects. The following results, both with HardwareFloats as well as with SoftwareFloats, are marred by numerical roundoff:

A := linalg::hilbert(15):
float(numeric::det(A, Symbolic)),
numeric::det(A, HardwareFloats),
numeric::det(A, SoftwareFloats)

How could I apply this,to check the result??Because I get an error message,when I write it...

 

[evinda]

The determinant of the Hilbert matrix is not zero, but it tends to zero exponentially fast as the size of the matrix increases. I know nothing about matlab, but I doubt whether it could possibly have enough precision to distinguish $\det H$ from zero for a $250\times 250$ Hilbert matrix. As is mentioned here, the Hilbert matrix is notoriously ill-conditioned for numerical computation.

Hasn't it maybe to do with the determinant,but with the condition number of the matrix?

 

[evinda]

Anyway... :) I also calculated the spectral radius of the hilbert matrix with dimension 250.Is this right that it is 217.3320,using the Jacobi method?

 

[Opalg]

Anyway... :) I also calculated the spectral radius of the hilbert matrix with dimension 250. Is this right that it is 217.3320,using the Jacobi method?

The Hilbert matrices are positive, so the spectral radius is equal to the norm, which is $\leqslant \pi$. Asymptotically, as the size of the matrix increases, the norm converges to $\pi.$ See Man-Duen Choi's paper "Tricks or Treats with the Hilbert Matrix".

 

[evinda]

The Hilbert matrices are positive, so the spectral radius is equal to the norm, which is $\leqslant \pi$. Asymptotically, as the size of the matrix increases, the norm converges to $\pi.$ See Man-Duen Choi's paper "Tricks or Treats with the Hilbert Matrix".

I mean the spectral radius of the iteration matrix..Is that what you told me valid for the Hilbert matrix or for the iteration matrix??

 

[Opalg]

I mean the spectral radius of the iteration matrix..Is that what you told me valid for the Hilbert matrix or for the iteration matrix??

I meant the Hilbert matrix. I know nothing about the iteration matrix.

 

[I like Serena]

If this were not clear yet, your algorithm is diverging.
As a result the values become bigger and bigger.
At some point they become so big that they are considered infinite (Inf). The same happens if you divide by zero somewhere.
Then, sooner or later, Inf is subtracted from Inf or something like that, resulting in the special value Not-a-Number (NaN).
Afterward, those values will usually remain NaN, although there are some unexpected cases where they may turn into a normal numbers again.

It's good practice to stop the algorithm when NaN occurs before the result becomes faulty.