Why matlab returns inaccurate response for rank of my matrix?

This is my matrix:

 A=
6.07787530038946e+15  0  0  0  0  0  0  0
0  6.07787530038946e+15  0  0  0  0  0  0
0  0  6.07787530038946e+15  0  0  0  0  0
0  0  0  1215578336877.89  0  0  0  0
0  0  0  0  1215578336877.89  0  0  0
0  0  0  0  0  0  0  0
0  0  0  0  0  0  -1.47264000000000e-12  0
0  0  0  0  0  0  0  -1.96352000000000e-16

When ask matlab to say the rank of A it returns 5, while it is 7 as you see. why it returns 5, and how to fix it?

Answer 1

[Walter Roberson]

If you have the symbolic toolbox then

rank(sym(A))

[Mehdi]

this was an example, the main matrix have very big dimentions (2000*2000), so sym is not good way to now the rank of such matrices.

Answer 2

[Jan]

If you have a 2000x2000 matrix with elements on the diagonal, which vary by a factor of 10^30 on the diagioal, and do not want to process the data symbolically, you are simply lost. This is the wrong method to solve your problem. It can neither be done symbolically nor numerically. The only solution will be to do it mathematically, which means using your brain.

You have to determine the meaning of the values during the creation already. You need a method to decide, if 1e-15 is a rounding error or a meaningful value. Note that this is not trivial, because then you have to consider 0 also as potential meaningful value, when it was created e.g. by 1e17 + 1 - 1e17.

A standard approach in numerics is scaling the values: Together with e.g. a vector of positions, another vector is created carrying scaling factors such that the numbers get about 1. Then e.g. a variation of the values is much easier, or the calculation of the Hessian matrix. Scaling is used in optimization and simulation of dynamic systems, example: The integration of the equations of motion of a lunar orbiter. To control vibrations of the engines, millimeters matter. Right at the start the position is known in the magnitud of centimeters, but when the lunar orbit is reached, such a precision would be nonsense. A scaling of the different components of the trajectory allows the integrator to consider them with a physically meaningful precision.

Summary: Think twice. Your problem cannot be solved by standard methods.

[Mehdi]

my matrix is not diagonal in general. What do you mean by scaling a matrix that have very big and very small elements. Please give an example how you scale those matrices.

Here I attached my matix

[Jan]

@Mehdi: You ask for an example? Did you read my answer?

Posting the array does not help here. As I've written before: You need to know the meaning of the elements in your case. If you want to distinguish a 1e-15 from a 0.0, you have to know how it was created and what the physical meaning is. Example:

sin(pi)

This is not 0.0, but 1.22e-16. This is not a bug, but the expected behavior considering the limited precision. Now you ask for a method, which can overcome the effects of the limited precision magically, but this cannot work. You tell the result of the rank() command "inaccurate". But if you want a different result, you have to defined an algorithm to find it. Currently based on the provided information, the result of rank() is correct, but your expectations are inaccurate.

Answer 3

[Star Strider]

The result of the rank function is correct. With respect to the first 5 values of diag(A), the last three are vanishingly small and evaluate to zero. For all intents and purposes, rank(A) is 5.

[Mehdi]

I know they are very small but I want matlab account them in rank(A). Is it possible to say matlab do not ignore them??? How?

[Star Strider]

With respect to matrix inversion, the small numbers are essentially zero. There is no way for MATLAB to ‘ignore’ that they are essentially zero. They are overwhelmed by the magnitudes of the other numbers.

[Steven Lord]

If the number in the upper-left corner represented Bill Gates's net worth, the number in the lower-right corner would represent something like a (very) small fraction of a penny. Is that lower-right number exactly 0? No. Is it, for almost all intents and purposes, 0? Yes.

Is it possible to force rank to treat those extremely small numbers (relative to the other entries in the matrix) as nonzero? Yes. Is that a good idea? Probably not. If you tell us WHY those very small numbers are important, we may be able to offer guidance on how to properly handle them.

[Walter Roberson]

If the number in the upper-left corner represents the size of the observable universe, then the number in the lower-right corner would represent about 30 microns.

Answer 4

[John D'Errico]

As others have said, just wanting magic to happen will not suffice. In MATLAB, you need either to use symbolic tools for that computation, or numerical tools. The numerical tools in MATLAB CANNOT be used to give you higher than double precision accuracy, because they are written in double precision. Your problem is wildly beyond the domain of double precision arithmetic.

Of course, if your problem is too large, then symbolic tools will also not suffice, as they will take far too long for any possible solution.

There is no law that says that ALL problems are solvable. If you choose to create a problem that is too big and with too wide of a dynamic range, then expect your problem to fall into the set that do not have a solution, at least in MATLAB. And tomorrow, we would probably find out that you have now decided to make your problem even more difficult, with a wider dynamic range yet, because if it was easy then anybody can do it.

Nothing stops you from writing or finding linear algebra code in some other language that supports a highly extended precision arithmetic. Even an implementation of quad precision would probably be insufficient here. But but be careful, as infinite precision tools (like the JAVA tools) tend to do poorly on problems like this, as those computations generate numbers with more and more digits at each step of the process.

A solution could be to write a version of rank that used my HPF tool, it would probably be fine on your problem if you decided to work in 50 or 100 decimal digits of precision. You could write a simple version of rank then, based on a column pivoted Q-less QR algorithm. Or use a rref variation, that will work on HPF inputs. (SVD would be slower, and take a bit more effort. rref would be easy enough to hack to allow HPF. I might even write one to put in HPF if I have some free time.) But at least using a variable precision tool like HPF, you can choose the precision to work in.

(Actually, rref is REALLY easy to hack to work with HPF inputs).

[John D'Errico]

I did write an HPF version of rank, and am testing it on your array. Not terribly fast, but it will work.