Precision of inverse matrix calculation

Hi,

I noticed that the two matrix operators 'divide by matrix' and 'multiply by inverse', which are theoretically identical, are NOT identical when calculated (for matrices with real numbers).

A / B  'might be almost but not equal to'  A * B^-1

By assuming the following two example matrices and evaluating both equations, results are equal for the first eleven digits following the decimal point and then deviate from each other. This may sound fussy and irrelevant. If some more calculations are however based on this one partial result, differences may increase exponentially.

A  = eye(2);
B  = [117, 822.2940998481383,
     822.2940998481383, 5783.818979511911];
C1 = A / B;
C2 = A * B^-1;

From my point of view, I would guess that two different algorithms are used. Anyone knows where this difference is?

Best, Chris

[Christoph]

Yeah, that's quite confusing. What I have in the back of my head is that the consensus was that A / B should be preferred.

[Roger Stafford]

The B matrix you have used is nearly singular. Take its determinant and see. That will make both kinds of computation prone to larger errors.

Answer 1

[John D'Errico]

Not confusing at all, or at least, it should not be so.

You don't want them to be the same. Computing the inverse and then multiplying by it is NOT the preferred choice. Instead, you should prefer to use slash (or backslash depending on the specific case). The slash code will be more efficient in general. It will be more numerically stable in general.

For a time comparison...

A = rand(10);
B = rand(10);
timeit(@() A*inv(B))
ans =
    3.242e-05

timeit(@() A/B)
ans =
    1.869e-05

Anyway, even if a subtly different order of additions and multiplies were used, you should expect a different result, due to floating point trash. For example, these two sets of operations yield different results in floating point arithmetic.

.3 - .2 - .1
ans =
  -2.7756e-17

-(.1 + .2 - .3)
ans =
  -5.5511e-17

So you cannot expect identical results from a pair of even slightly different algorithms in floating point arithmetic.

Floating point arithmetic is not truly mathematics, even though the two look a lot the same. With care, floating point arithmetic can approach the behavior you would expect from mathematics, but that does take care, and an understanding of what deviations one can expect.

[Chris]

Thanks for pointing that out!