Explosive GEMM!
A demonstration that arbitrarily large FP error can be incurred in the GEMM operation.
GEMM (GEneral Matrix Multiplication) computes:
For the purposes of this example, we only need to consider expressions
Suppose f and g are two GEMM computations equivalent up to associativity of arithmetic.
We will show that the difference f(α, A, B) - g(α, A, B) can be made arbitrarily large.
Run
; python explode.py
Matrix A:
[[1.0, 1.0, 1.0]]
Matrix B:
[[1e+16], [-1e+16], [1.0]]
α
1e+100
Left-associative result:
1e+100
Right-associative result:
0.0
How it works
Floating point is not associative. This means that in general,
(x + y) + z ≠ x + (y + z)
A basic example:
x = 1e16 y = -1e16 z = 1.0 print((x + y) + z) # prints 1.0 print(x + (y + z)) # prints 0.0
This gives an error of 1.0 between the two sums.
We can arbitrarily inflate this by multiplying by a large scalar,
leading to arbitrarily large error:
x = 1e16 y = -1e16 z = 1.0 s = 1e100 print(s * ((x + y) + z)) # 1e100 print(s * (x + (y + z))) # 0.0
We can then encode this into a single GEMM operation
using
α = 1e100
A = [1, 1, 1]
B = [x, y, z].T