Warning: Matrix is close to singular or badly scaled.

Hello I am playing around with a code for that has been uploaded by Biswa Datta https://uk.mathworks.com/matlabcentral/fileexchange/2155-numerical-linear-algebra-and-applications/content/datta/CHAP8/rayqot.m. If I use

B=[16 2 3 13;5 11 10 8;9 7 6 12;4 14 15 1];
ep=0.0000005;
nmax=1000;
x0=ones(n,1);

I get this message: Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.691450e-17.

How would one go about this problem?

[John BG]

JR

the link you have supplied in your question returns

would it be possible for you to supply code instead?

regards

John BG <mailto:[e-mail] [e-mail]>

[Jan]

@John BG: I've removed the parenthesis from the link and in is working now.

@JR: What is n? Do you call the function like this:

[sigma,x,iter] = rayqot(B, x0, ep, nmax)

[John D'Errico]

But the link you show and the small amount of code you supply does not include how you are trying to call the function. Yes, I think I know how you are calling it. But be helpful to others, making it easy for them to help you.

[JR]

My apologies, I will put more effort when I post questions in the future!

@Jan Simon: many thanks for fixing it!

[Jan]

@JR: You are welcome. Questions for clarifications belong to the standard process of solving problems.

Answer 1

[John D'Errico]

How would you go about it?

If you call this function as it is supposed to be called, you will get this output:

[sigma,x,iter] = rayqot(B,x0,ep,nmax)
> In rayqot (line 24) 
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.691450e-17. 
sigma =
    34
x =
                         1
                         1
                         1
                         1
iter =
     0

So, what happened? Lets compute the eigenvalues/vectors of B.

[V,D] = eig(B)
V =
         -0.5     -0.82361      0.37639     -0.22361
         -0.5      0.42361     0.023607     -0.67082
         -0.5     0.023607      0.42361      0.67082
         -0.5      0.37639     -0.82361      0.22361
D =
           34            0            0            0
            0       8.9443            0            0
            0            0      -8.9443            0
            0            0            0   2.0576e-15

As it turns out, you started with an initial vector that happened to be a true eigenvector of B. That blew this code out of the water. :) The best code will try to be robust to problems.

If, instead, we start with a random vector for x0, see what happens:

[sigma,x,iter] = rayqot(B,rand(4,1),ep,nmax)
sigma =
          33.9999999999865
x =
     1
     1
     1
     1
iter =
     3

So, your problem was you were unlucky in choosing x0. Or, we can say that the author was sloppy in writing the code, so that it was not robust to this very problem.

Must I go into more depth on this? Perhaps it might be illuminating if I did. I'll use the debugger.

At the very first iteration of this code when you pass in the actual eigenvector of B, we see

sigma = (x'*A*x)/(x'*x)
sigma =
    34

Remember that x0 was a vector of all ones, so simply a scaled version of the eigenvector we are looking for. sigma is returned on the very first iteration as 34, which is EXACTLY the desired eigenvalue. What else would you expect there?

So now, the very next line has you forming this matrix:

A-sigma * eye(n,n)
ans =
   -18     2     3    13
     5   -23    10     8
     9     7   -28    12
     4    14    15   -33
K>> rank(ans)
ans =
     3

Again, it will be singular. So backslash operation in that line will see MATLAB get upset, as it should!

xhat = (A-sigma * eye(n,n))\x;
> In rayqot (line 24) 
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.691450e-17. 

xhat
xhat =
   5.6295e+14
   5.6295e+14
   5.6295e+14
   5.6295e+14

Near infinite garbage, as you should expect here.

How should the code have been written to be more robust? A very simple solution would have been to test if the initial vector was indeed already an eigenvector, just stopping immediately.

So it is good that you were playing around with this code. Much can be learned by playing around in MATLAB.

[JR]

Thank you kindly for your help John! So for clarification,

1. the backlash operator in MATLAB is not ideal for matrices that are not full-rank? In the general case, how would one go about it?

[John D'Errico]

Ugh. The answer to your question could require a class in numerical linear algebra, coupled with a numerical analysis class. We would need to talk about vector spaces, about subspaces, about the column space of A, of shoes and ships and sealing wax, of cabbages and kings... Oops. Sorry. I started to channel Lewis Carroll. :)

It is true that if you try to solve the problem A*x=b using the operation A\b (when A is not full rank) then backslash is often not the best method. But "best" is a very subjective word.

The problem is the solution in that case is not well posed. Assume that A is square. Then you could write this using Gaussian elimination. At the end, when A is numerically singular, you would see a divide by zero. Depending on the problem, that could end up as a case of 0/0, or it could be a case of something/0. 0/0 is nan. k/0, for non-zero k will yield +/- inf.

The problem is, you generally don't see an exact zero there. In fact, that will generally not happen.

So you often see very large numbers as a result, on the order of 1/eps.

In general, backslash IS the proper tool to use for this, with the caveat that if it tells you the problem is nearly singular, then take a step back. That warning says you are trying to solve a problem that was probably not well posed. The thing is, unless the divisor is exactly zero, MATLAB does not know for sure that you have a truly singular problem, or something close to that, but one that you really want/need to solve. So MATLAB throws a warning.

So, when backslash burps at you, look at the problem you are trying to solve. First of all, ask yourself why is the problem numerically singular. Most of the time, that reflects something wrong in your problem formulation. It may be something simple or not. In some cases, you may have a mathematically non-singular problem, yet numerically, the problem is insolvable in floating point arithmetic. The thing is, when that happens, you almost never have data that is sufficiently accurate to produce a solution even if infinite precision were employed.

When you have concluded that your problem really merits a solution as is, and there is a good reason to need to solve the problem DESPITE the numerical singularity, then you can use pinv. Thus, instead of A\b, use pinv(A)*B.

Pinv has some different properties than backslash. One reason why pinv is not used as a default always is that it will be slower, sometimes significantly slower. Nobody wants slow code.

So would pinv be the "best" solution, better than backslash? Or perhaps lsqr, an alternative we have not discussed? Again, best is only definable in the eyes of the user.

Again, if backslash burps, only then should you worry. Or you could use use rank or cond to test the problem before throwing it at backslash.

In context of the code that you were trying to use, backslash was entirely reasonable to use, as long as x0 was not an eigenvector. Truly high quality numerical code is difficult to write. That requires the author thinking about every possible mode of error. What might the user have done to mess it up? As I suggested, better code might have done some input checks. Is x0 an eigenvector already? The point is, if it is, then the solution is to avoid the iterative computation.

Answer 2

[Jan]

ones(n,1) is a "bad" initial approximation to the eigen vector, because it is the result also. In the first iteration

sigma = (x'*A*x)/(x'*x);
xhat = (A-sigma * eye(n,n))\x;

solves a system, which is close to singular. When the initial guess is too good, the algorithm shows the warning, but still calculated the correct result.

Either choose a different initial guess or take into account that the first step has an inaccurate result.