Matrix multiplication inverse rounding error in Julia

Question

julia> A = [1 2; 3 4]
2x2 Array{Int64,2}:
 1  2
 3  4

julia> A * inv(A)
2x2 Array{Float64,2}:
 1.0          0.0
 8.88178e-16  1.0

Why is the result in the lower left corner 8.88178e-16 rather than 0.0 like in the upper right? How can I perform the operation so that I get [1.0 0.0; 0.0 1.0] as the result.

Is there a way to display the result as desired, but actually perform the matrix operation symbolically as one would be able to in muPad or some other symbolic evaluator?

That's how floating point numbers work. 64 bits of IEEE doubles gives 16 digits of precision. https://en.wikipedia.org/wiki/IEEE_floating_point — duffymo, Oct 01 '15 at 00:24
Welcome to the world of floating-point computations. This behaviour is normal. You may implement your own matrix-operations with a decimal-like number object, but there is a reason, why floating-point numbers are dominating the scientific world — sascha, Oct 01 '15 at 00:24
You can use `SymPy` if you want symbolic math in `Julia`: `using SymPy; A = [Sym(1) 2; 3 4]; A * inv(A)`. (This just passes the task to Python's SymPy module.) — jverzani, Oct 01 '15 at 00:54
Possible duplicate of [What mistake did I make in this matrix multiplication in Julia?](http://stackoverflow.com/questions/30552853/what-mistake-did-i-make-in-this-matrix-multiplication-in-julia) — mbauman, Oct 01 '15 at 03:26
@jverzani: Alternatively, Arthur could use the native Rational type: ``A = Rational[1 2; 3 4]; A*inv(A)`` — Moritz Schauer, Oct 01 '15 at 09:22
@mschauer: I tried your suggestion but got this error (ERROR: no promotion exists for Rational{T<:Integer} and Float64) — Arthur Collé, Oct 06 '15 at 23:29
@ArthurCollé Does it work with ``A = Rational{Int}[1 2; 3 4]; Ainv=float(inv(A))``? — Moritz Schauer, Oct 07 '15 at 09:53
Keep in mind the option of using rounding after performing the operations in floating point, which may be satisfactory in many use cases and will be pretty efficient computationally. — Michael Ohlrogge, Aug 26 '16 at 16:20

score 4 · Answer 1 · edited Aug 26 '16 at 16:13

As mentioned in the comments, the best way is to use Rational type. E.x., A = Rational{Int}[1 2;3 4], to get the output in your desired format use, A*float(inv(A)). Credits to @mschauer.

To elaborate a little, the floating point numbers typically used in computation have limited precision, and cannot represent all real numbers. Hence, you get these kinds of errors when operating with floating point numbers. See links in the comments for details on this issue.

Julia however makes it very easy to define other numeric types, and operate on them if you want. One of them is a rational number type, which can represent fractions exactly. We can create the array as rationals:

julia> A = Rational{Int}[1 2;3 4] 2x2 Array{Rational{Int64},2}: 1//1 2//1 3//1 4//1

Then, operations on this array will return exact rational results

julia> A*(inv(A)) 2x2 Array{Rational{Int64},2}: 1//1 0//1 0//1 1//1

If we want the results as floating point numbers, they can be converted at the end of the computation.

julia> float(A*(inv(A))) 2x2 Array{Float64,2}: 1.0 0.0 0.0 1.0

It is important to note that the reason floating point numbers are used by default is performance. CPU's are optimised to operate on floating point numbers, and using rationals as above is going to be much slower. Julia however gives you the tools to allow you to make that choice yourself, based on your own needs.

Matrix multiplication inverse rounding error in Julia

1 Answers1