Estimating precision of matlabs matrix solver

Hi friends of the MATLAB community,

I'm interested in knowing the precision of the matrix solution using matlabs

"x = A\b" solver.

The "A" I am working with is a non-singular, symmetric n x n matrix (that is typically 100x100, but I am interested in the general case). Elements in this matrix are either 0 or some number between 10^-8 and 10^-4. There are typically around 7n non-zero elements (3n in upper diagonal, 3n in lower diagonal, all diagonal elements are non-zero).

The "b" I am working with is a vector of length n with all elements zero except one or two elements which are numbers between -5 and 5 (of some arbitrary order of magnitude).

Does anyone have any ideas of where I can look or a method that I could use to calculate this?

[Bruno Luong]

Read literature about "condition number"

Answer 1

[John D'Errico]

The precision? Double precision. It works in double precision. But that is not your real question. The issue is how good are the estimates of your solution? And that is not a precision issue, but of understanding first linear algebra, as well as understanding floating point arithmetic, and perhaps most important, is to understand the condition number of your matrix and what that tells you.

What am I talking about? When you solve a linear system of equations, there are numerical errors that arise. First, your right hand side vector (b) is not known exactly. It is stored in floating point arithmetic, so there are errors in those numbers that are on the order of eps(b). That is, the least significant bit in those numbers may be an issue. Remember that a number is stored in double precision using BINARY storage. For example, what is 0.05?

>> sprintf('%0.55f',0.05)
ans =
    '0.0500000000000000027755575615628913510590791702270507812'
>> eps(0.05)
ans =
   6.9389e-18

It turns out that most decimal fractions are not exactly representable in binary, just as 1/3 is not exactly representable in a finite number of decimal digits, the same is true for many of the numbers you will encounter in MATLAB.

So when you create a vector b in MATLAB, in fact, you store only an approximation to b. And that means, when you try to solve a problem in the form A*x=b, you are actually trying to solve the related problem

(A + A_delta)*x = (b+b_delta)

where A_delta and b_delta are a matrix and vector composed of tiny numbers, but are in general both non-zero. (And this is part of the reason why you can benefit from classes in numerical methods, as well as classes in linear algebra.)

So now what happens? It turns out that the solution x will not be exact. In fact, the solution can amplifiy those tiny errors. But worse, often you do not know the true value of the right hand side vector. It may be composed of data that has true noise in it. As I said though, the solution is a noise amplification process, and the amplification factor can be as large as cond(A). That is, the condition number of your matrix will tell you much.

For example, suppose I just create a perfectly random matrix A?

>> A = rand(100,100);
>> cond(A)
ans =
       1860.5

In this case, the condition number was 1860.5. Not huge. But if the matrix was singular, then the condition number would be theoretically infinite. In floating point arithmetic, that means it will be on the order of 1e16.

>> A(2,:) = A(1,:);
>> cond(A)
ans =
   1.6594e+17

So when I replicated the first row of A into the second row, the matrix is now singular. Trying to solve a singular system of equations will essentially amplify the error so much that the results will always be numerical garbage.

That means, what matters is not how precisie is the matrix solver, but it tells you that what really matters is the condition number of your matrix A.

>> A = rand(100,100);
>> cond(A)
ans =
       3161.6
>> xtrue = rand(100,1);
>> b = A*xtrue;

Now I've recreated a new, non-singular random matrix A, and use it to solve the system of equations. We know that b should be exactly correct, right? (NOT so, because b is composed of numbers in floating point, created from floating point computations. b is not exactly correct, even if you think it is.)

>> xsol = A\b;
>> norm(xtrue - xsol)
ans =
   2.7452e-13
>> eps
ans =
   2.2204e-16

So, the solution amplified even the tiny (virtually random) errors in b by a factor that was almost as large as cond(A). This is NOT the precision of the solver, but of basic mathematics and linear algebra.

As I said above, you can and should take entire courses in numerical methods, numerical linear algebra. And this only scratches the surface. But the number you will often care to know here is what is cond(A), for YOUR problem. We don't know that number, since you never gave us A.

Note that when A is perfectly well-behaved, perfectly well conditioned, then the condition number of A will be 1.

>> cond(eye(100))
ans =
     1

It cannot get any better than that. You will rarely, if ever see a condition number that small however.

[Joel Hochstetter]

Hi John,

Thanks for your detailed answer. I understand the points you are trying to make.

For reference the condition number of my matrix is typically ~10^9.

"So when you create a vector b in MATLAB, in fact, you store only an approximation to b. And that means, when you try to solve a problem in the form A*x=b, you are actually trying to solve the related problem

(A + A_delta)*x = (b+b_delta)"

I have an estimate on the size of the uncertainty of each element in A and b, that is larger than the limitations of double precision. I have an estimate on the uncertainty of each element in b.

I would like a naive estimate on the uncertainty of x, given A_delta and b_delta (there are absolute errors).

Is it fair to say that

i.e. Condition number linearly scales the percentage uncertainty in A.

Then I can take maximal percentage uncertainties in A and B to get a naive estimate on the percentage error of x