Julia: Reshape fails to exactly commute with matrix-vector product

[thisrod]

              _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.0-rc2.0 (2020-02-24)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> x = randn(8);

julia> M = randn(7,8);

julia> M*reshape(x,8,1) == reshape(M*x,7,1)
false

julia> M*reshape(x,8,1) - reshape(M*x,7,1)
7×1 Array{Float64,2}:
 0.0
 0.0
 2.220446049250313e-16
 0.0
 1.1102230246251565e-16
 0.0
 1.1102230246251565e-16

julia> versioninfo()
Julia Version 1.4.0-rc2.0
Commit b99ed72c95 (2020-02-24 16:51 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.0.0)
  CPU: Intel(R) Core(TM)2 Duo CPU     E7600  @ 3.06GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, penryn)

[dkarrasch]

This is because in one case (M*reshape(x,8,1)) you're performing a matrix-matrix product which is handled by BLAS-level-3 methods, and in the other (M*x) a matrix-vector product which is handled by BLAS-level-2 methods. The fact that the results differ by something on the order of machine precision is no surprise.

This has nothing to do with the reshape, except for that the reshape reinterprets the vector as a matrix.

[thisrod]

It might not surprise people who work on Julia's internal linear algebra implementation, but it surprised me. Would it be worth having a shortcut so that Julia uses BLAS-level-1 and BLAS-level-2 methods for row and column matrices?

[oscardssmith]

In general, you should never expect floating point arithmetic to be more accurate than 10e-14. The bigger issue isn't which blas method Julia uses, but that you expect floating point anything to be numerically perfect. This will (almost) never happen.

[dkarrasch]

It might not surprise people who work on Julia's internal linear algebra implementation

This is not about Julia's internals. I just wanted to point out that really two different algorithms are at work, for which, as @oscardssmith pointed out, you can't expect the exact same result in floating point numbers. I'm sure that any programming language that uses BLAS in the background for the heavy linalg lifting has the same "issue". In particular, the reshape function does not change values, which was the suspicion of the OP.

Would it be worth having a shortcut so that Julia uses BLAS-level-1 and BLAS-level-2 methods for row and column matrices?

I don't think we want to interfere the multiple dispatch process to check whether a "matrix" has only one column and then redirect to matrix-vector multiplication. Importantly, it's not true that one of the computations yields "the correct" the result and the other one is "wrong". Both are strictly speaking wrong, but both are reasonably good.

[timholy]

To demystify it a little bit, this comes down to the non-associativity of floating-point arithmetic:

julia> (0.1+0.2)+0.3
0.6000000000000001

julia> 0.1+(0.2+0.3)
0.6

You get the same result in Python & Matlab, although Matlab hides that fact from you by rounding when printing, making you vulnerable to hard-to-discover bugs if you happen to test result1 == result2 and then can't figure out why they aren't equal despite the fact that they print identically.

When you're multiplying matrices, you have lots of additions (corresponding to the different columns of the leading matrix), and the order in which you compute them matters for ulp-level accuracy. Different level BLAS algorithms make different choices. That's why most people ignore errors of this scale.

[thisrod]

Thanks for spending your time on this