Dense matrix algebra has the benefit that besides the matrix dimensions, the data itself has no impact on performance

I was very surprised when I first learned that's not strictly true, because of denormal numbers[1]. Here's a session in ipython --pylab (with MKL), demonstrating a 200x slow-down for matrix multiplication with tiny numbers. Crazy!

    In [1]: A = randn(1000, 1000)
    
    In [2]: %timeit A @ A
    25.9 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    
    In [3]: A *= 1e-160
    
    In [4]: %timeit A @ A
    5.52 s ± 21.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

You hit denormal numbers more quickly with single-precision floats. I have been bitten by this issue in the wild a couple of times now, and seen a couple of other people with it too. Sometimes denormals are created internally in algorithms, when you didn't think your input matrices were that small.

[1] https://en.wikipedia.org/wiki/Denormal_number#Performance_is...