# Matrix multiplication precision issue
**URL:** https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172
**Category:** General Usage
**Created:** [January 28, 2019, 7:09am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172 "2019-01-28T07:09:08Z")
**Posts on this page:** 14
**Page:** 1
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 28, 2019, 7:09am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/1 "2019-01-28T07:09:08Z")
Why A1 and A2 have different precision and A1 != A2, though mathematically, they should have the same results. I’m wondering is there anything special about matrix slicing that would affect the precision? or the results are caused by something esle?
Example(Julia 0.6.4):
A = randn(Float64, 200, 100)
A = [ones(Float64, 200) A]
B = randn(Float64, 200, 1)
ApB = A’B
idx = [1, 2]
A1 = ApB[idx]
A2 = (A[:, idx]’ \* B)[:]
A1 == A2
A2[1] == sum(B)
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [January 28, 2019, 8:13am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/2 "2019-01-28T08:13:58Z")
> [@Caesar\_KG](#):
>
> mathematically, they should have the same results
> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678):
>
> Sometimes people are surprised by the results of floating-point calculations such as julia\> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia\> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…
Also, please quote code
> [@PSA: how to quote code with backticks](https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530):
>
> This is a short post on how to use backticks (` and ```) to quote code, so it is easy to read. This post can be linked to new users who are confused by this feature. Why By quoting your code, you get a monospaced font (which preserves indentation) and syntax highlighting. This makes it easier to read your code and help you. Displayed code looks like this: function displayed\_code(x::Int, y::Int) if x \< y println("x is smaller") else println("or not") end end Inline…
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [January 28, 2019, 8:18am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/3 "2019-01-28T08:18:42Z")
I think you’re just hitting rounding errors (as already said above). You may find it interesting to try the same thing using `DoubleFloats` (see below). Note that if you use `≈` instead of `==` in your example it also “works” but that doesn’t answer your question I guess.
```julia
julia> using DoubleFloats
julia> using Random
julia> Random.seed!(123);
julia> A = hcat(ones(200), randn(200, 100));
julia> B = randn(200);
julia> # let's declare objects in double double precision
julia> A2 = convert(Matrix{Double64}, A);
julia> B2 = convert(Vector{Double64}, B);
julia> # Now let's do the computations as you did
julia> ApB = A'B;
julia> idx = [1, 2];
julia> A_1 = ApB[idx];
julia> A_2 = (A[:, idx]' * B)[:];
julia> # the results are not identical (floating point error)
julia> A_1 - A_2
2-element Array{Float64,1}:
0.0
1.9539925233402755e-14
julia> # let's do the same thing in double double...
julia> ApB2 = A2'B2;
julia> A_12 = ApB2[idx];
julia> A_22 = (A2[:, idx]' * B2)[:];
julia> # the rounding error goes away
julia> A_12 - A_22
2-element Array{Double64,1}:
0.0
0.0
```
[link to the awesome DoubleFloats package repo](https://github.com/JuliaMath/DoubleFloats.jl)
---
### Author:  [@bennedich](https://discourse.julialang.org/u/bennedich)
#### Post date: [January 28, 2019, 8:58am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/4 "2019-01-28T08:58:36Z")
Rounding issues, as others have said. Use this to calculate the relative norm (“error”):
```julia
julia> norm(A1-A2) / max(norm(A1), norm(A2))
2.0277114478533147e-15
```
Seems small enough. And here’s what you can use instead of `==`:
```julia
julia> isapprox(A1, A2; rtol = 1e-10)
true
julia> isapprox(A1, A2; rtol = 1e-15)
false
```
Or, shorter but perhaps less readable:
```julia
julia> ≈(A1, A2; rtol = 1e-10)
true
julia> ≈(A1, A2; rtol = 1e-15)
false
```
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 2:28am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/5 "2019-01-29T02:28:19Z")
I’ve read the previous post already. But I still don’t understand. Aren’t the numbers in the first two column of matrix A used for calculation the same in both cases? If they are same, the results should be the same following the floating point arithmetic rules.
BTW, thanks for the tip regarding code quotation.
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 2:59am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/6 "2019-01-29T02:59:49Z")
Thanks for the tip on DoubleFloats. Really helps!
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [January 29, 2019, 3:39am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/7 "2019-01-29T03:39:49Z")
I think the reason is to be found behind ordering of operations. Julia may re-order operations to get optimal performance and if you sum floating point numbers in different orders you may get different results.
Small experiment below:
```julia
using Random
using LinearAlgebra
Random.seed!(125)
A = hcat(ones(200), randn(200, 100));
B = randn(200);
ApB = A'B;
idx = [1, 2];
# -- let's unroll the operations partially
ApB = Vector{Float64}(undef, 101)
@inbounds for i ∈ 1:101
ApB[i] = dot(view(A, :, i), B)
end
A2 = Vector{Float64}(undef, 2)
@inbounds for i ∈ idx
A2[i] = dot(view(A, :, i), B)
end
ApB[idx] == A2 # true as expected
# -- now let's do the same thing but using
# -- a different ordering for ApB and for A2
rp1 = randperm(200)
rp2 = randperm(200)
ApB = Vector{Float64}(undef, 101)
@inbounds for i ∈ 1:101
ApB[i] = dot(view(A, rp1, i), B[rp1])
end
A2 = Vector{Float64}(undef, 2)
@inbounds for i ∈ idx
A2[i] = dot(view(A, rp2, i), B[rp2])
end
ApB[idx] == A2 # false !
ApB[idx] ≈ A2 # still true
```
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 6:16am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/8 "2019-01-29T06:16:54Z")
Thanks for the help, man!
That may explain the issue. It’s related to the Julia implementation of floating point arithmetic.
I tried it on python and A1==A2 is true.
```python
import numpy as np
import numpy.random as nprd
In [27]: A = np.c_[np.ones((200, 1)), nprd.randn(200, 100)]
In [28]: B = nprd.randn(200, 1)
In [29]: ApB = np.dot(A.T, B)
In [30]: idx = [1, 2]
In [31]: A1 = ApB[idx]
In [32]: A2 = np.dot(A[:, idx].T, B)
In [33]: A1 == A2
Out[33]:
array([[True],
[True]])
In [34]: type(A1[0][0])
Out[34]: numpy.float64
```
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [January 29, 2019, 6:29am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/9 "2019-01-29T06:29:24Z")
> [@Caesar\_KG](#):
>
> That may explain the issue. It’s related to the Julia implementation of floating point arithmetic.
> I tried it on python and A1==A2 is true.
that’s misleading (or you got lucky). I just had a shot and without surprise
```julia
>>> A1 - A2
array([[-1.77635684e-15],
[-5.32907052e-15]])
```
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 6:44am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/10 "2019-01-29T06:44:49Z")
OK. That’s weird.
Mine is:
```python
In [27]: A = np.c_[np.ones((200, 1)), nprd.randn(200, 100)]
In [28]: B = nprd.randn(200, 1)
In [29]: ApB = np.dot(A.T, B)
In [30]: idx = [1, 2]
In [31]: A1 = ApB[idx]
In [32]: A2 = np.dot(A[:, idx].T, B)
In [33]: A1 == A2
Out[33]:
array([[True],
[True]])
In [34]: type(A1[0][0])
Out[34]: numpy.float64
In [36]: A1 - A2
Out[36]:
array([[0.],
[0.]])
```
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [January 29, 2019, 7:18am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/11 "2019-01-29T07:18:52Z")
Have you tried several times? I set the seed to `1515` (at least doing that here gives me a standard epsilon machine error) anyway.
Another thing to note is that your Python might use MKL as a BLAS backend (mine does) where your Julia would use OpenBLAS, slight differences would then be expected.
PS: as a side note, I’m definitely not an expert but I’m pretty sure the Julia devs didn’t implement their own sauce for floating point arithmetics. There is a standard ([ieee754](https://en.wikipedia.org/wiki/IEEE_754)) which is most likely scrupulously followed. However the ordering of how operations are effectively done in a matrix-vector multiplication can depend (I think) on what BLAS backend you use.
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 7:35am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/12 "2019-01-29T07:35:19Z")
Yes, tried several times, the results are same in Python.
BTW, if the different results are caused by the order of computation, does it mean the implementations of whole matrix multiplication and block matrix multiplication are different in Julia?
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [January 29, 2019, 7:39am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/13 "2019-01-29T07:39:08Z")
No, that’s what I was trying to write (probably clumsily) above. These core linear algebra things (typically matrix-matrix and matrix-vector as far as I know) will be done by a BLAS backend. See also [the LinearAlgebra doc](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/), the “BLAS Functions” part.
See for instance in [the LinearAlgebra code](https://github.com/JuliaLang/julia/blob/6175bd9e3996f55f3e2a7bf8bf350bb7c96980ad/stdlib/LinearAlgebra/src/matmul.jl) you’ll see lines like `GC.@preserve x y BLAS.dot` etc. These call your BLAS backend and results might differ on what implementation of BLAS you’re using.
---
### Author:  [@Caesar\_KG](https://discourse.julialang.org/u/Caesar_KG)
#### Post date: [January 29, 2019, 8:35am UTC](https://discourse.julialang.org/t/matrix-multiplication-precision-issue/20172/14 "2019-01-29T08:35:20Z")
OK. I got it. Thx man!