What mistake did I make in this matrix multiplication in Julia?

Question

In Julia:

In  [1]: M1 = [1 3 4;
              45 64 33;
              456 3 454;]

Out [1]: 3x3 Array{Int64,2}:
   1   3    4
  45  64   33
 456   3  454

In  [2]: M1 * inv(M1)

Out [2]: 3x3 Array{Float64,2}:
  1.0           6.93889e-18  -8.67362e-19
  0.0           1.0          -2.08167e-17
 -1.42109e-14  -8.88178e-16   1.0

M1 * inv(M1) is supposed to get the Identity matrix by definition. What's wrong?

I tried the same thing in Matlab:

>> M1 = [1 3 4;
         45 64 33;
        456 3 454;]
M1 =
     1     3     4
    45    64    33
   456     3   454
>> inv(M1)
ans =
  -0.280088987764182   0.013057987135465   0.001518595540939
   0.052057842046719   0.013251438796731  -0.001421869710306
   0.280978865406007  -0.013203075881414   0.000686753397495
>> M1 * inv(M1)
ans =
   1.000000000000000   0.000000000000000  -0.000000000000000
                   0   1.000000000000000  -0.000000000000000
  -0.000000000000014  -0.000000000000001   1.000000000000000
>>

Matlab returns the right result here. I guess Julia will not make a mistake here. So what's wrong with my calculation / notation?

Edit

The problem is caused by number of digits in floating point result. I should have asked, how to set result digits precision in Julia?

Related question: http://stackoverflow.com/questions/27786864/why-am-i-getting-incorrect-results-for-these-simple-operations-in-julia — Vincent Zoonekynd, May 31 '15 at 01:46

score 0 · Accepted Answer · answered May 31 '15 at 01:42

0

Julia and Matlab actually give the same result (for instance, the bottom-left element is -1.4e-14 in both cases): it is not exactly the identity matrix because floating point arithmetic is not exact.

You can explicitly round the result before displaying it.

M1 = [
  1    3   4;
  45  64  33;
  456  3 454
]
round( M1 * inv(M1), 6 )
# 3x3 Array{Float64,2}:
#  1.0   0.0  -0.0
#  0.0   1.0  -0.0
#  0.0  -0.0   1.0

If you want an exact result, you can also use rationals.

M1 = [
  1//1  3   4;
  45   64  33;
  456   3 454
]
M1 * inv(M1)
# 3x3 Array{Rational{Int64},2}:
#  1//1  0//1  0//1
#  0//1  1//1  0//1
#  0//1  0//1  1//1

answered May 31 '15 at 01:42

Vincent Zoonekynd

31,893
5
69
78

The usual reference for floating point arithmetic is [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) (D. Goldberg, 1991), but I found [Practically Accurate Floating-Point Math](http://jeapostrophe.github.io/home/static/tm-cise2014.pdf) (N. Toronto and J. McCarthy, 2014) easier to read. – Vincent Zoonekynd May 31 '15 at 03:08
The `Practically Accurate Floating-Point Math` looks very neat, and the author John McCarthy is well-known in the computer science field. Definitely a good read! Thank you! – Nick May 31 '15 at 03:16
Oh, it's Jay McCarthy, not John McCarthy. The doc is still great. – Nick May 31 '15 at 03:19

What mistake did I make in this matrix multiplication in Julia?

1 Answers1

Linked