We experienced an issue with Matrix-Vector (or row-vector times column-vector) product. In the example the result should return zero: a = [1.3190624946305824, 1.3190624946305824] b = [-3.74 3.74] # Matrix-Vector-Product b * a # not always zero sum(b .* a) # zero On Machine 1, only the se…

[image] PSA: floating-point arithmetic Offtopic 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.2204460492503…

[image] judober: -3.583133801120832e-16 I find the deviation surprisingly large. Could you check with @code_native b * a if the same code is generated on both machines? If yes then maybe it is time to study CPU errata…

[image] judober: Given that matlab uses MKL we expected an issue there but using MKL.jl does not change the results. I think MKL.jl should give the same result as MKL in MATLAB (assuming it does use MKL, and the exact same version on the same processor). I think you’re not actually getting 0.…

[image] Palli: julia> isapprox(-3.583133801120832e-16, 0.0; atol=0.05) true But julia> 1 + -3.583133801120832e-16 0.9999999999999997 julia> nextfloat(nextfloat(nextfloat(1 + -3.583133801120832e-16))) 1.0 i.e. 3 ulp?

[image] Palli: while not: julia> -3.583133801120832e-16 ≈ 0.0 false https://twitter.com/MoseGiordano/status/1495208655653937152

[image] Palli: julia> -3.583133801120832e-16 ≈ 0.0 false From the isapprox documentation : Note that x ≈ 0 (i.e., comparing to zero with the default tolerances) is equivalent to x == 0 since the default atol is 0 . In such cases, you should either supply an appropriate atol (or use norm(x) ≤…

[image] goerch: deviation surprisingly large. Could you check with @code_native b * a if the same code is generated on both machines? If yes then maybe it is time to study CPU errata… I’m not sure that’s a reliable test. I.e. would you get exact same code (or very similar), on the Julia side,…

[image] Palli: If the deviation is actually considered too large (3 ulps) then could it rather be a software bug? AFAIU the cancellation error could be larger than 3 ulps. But that would mean the products should differ somewhere, which I don’t see: julia> reinterpret(UInt64, 1.31906249463058…

[image] Palli: You can use @Elrod ’s pure Julia code, but I’m not sure it’s NOT special-cased too (or if not yet, then might be later). Architectures with FMA will use it, those without will not. Block sizes differ, which will also result in different orderings for larger matrices.

And that’s what happening here it seems, I at least get the exact value with: julia> fma(1.3190624946305824, -3.74, 1.3190624946305824 * 3.74) 3.583133801120832e-16 So to be “expected”?

Inaccurate Matrix Multiplication

General Usage

tim.holy February 22, 2022, 11:31am 2

Topic		Replies	Views
Mul! generates machine noise General Usage matrices	21	718	February 21, 2024
X * y + z does not automatically use FMA instruction Performance fast-math	31	2915	July 27, 2023
Benchmark MATLAB & Julia for Matrix Operations Performance	148	20789	October 15, 2019
Why mul! is so fast? General Usage question , linearalgebra	7	7052	November 26, 2019
Vector - Matrix - Vector multiplication Performance	19	4648	March 14, 2021

Inaccurate Matrix Multiplication

Related topics