Dense matrix algebra has the benefit that besides the matrix dimensions, the data itself has no impact on performance: the access pattern is entirely predictable and known as soon as you know the matrix sizes (side note: which is why just-in-time libraries such as libxsmm exist). For most cases, you can rely on MKL or Eigen to produce near-optimal performance for your dense linear algebra operations.

Sparse matrices do not have that luxury. Sparse matrix ops have access patterns that entirely depend on your data. You simply cannot make a library that covers every possible distribution of non-zeroes in an optimal way. At best you have a toolbox of a large number of representations and algorithms that you choose for each different use-case. Engineering simulations and solvers can have block or cluster patterns of non-zeroes, which can be exploited by doing dense operations for the dense blocks. Large real-world graphs often have an exponential non-zero distribution and are often extremely sparse, which might need an hash-based join algorithm for an axpy-kernel. Are you running a tight CG-solver loop? Then you might want to bundle multiple spmv-kernels to operate multiple vectors on the same matrix (for similar reasons as states in the original article). You are distributing with Spark? Unless your non-zero distribution is uniform, your matrix is really hard to partition evenly. If your matrix never changes, you may consider preprocessing using a partitioner (METIS, Mondriaan).

In short, it is highly unlikely that the library you are using is anywhere near optimal for your use case. Do explore alternatives based on your domain-insights, it is not uncommon to net double digit performance improvements.

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...

Yes! For some applications including machine learning you can set the option to flush denormals to zero so you don't hit this case.

Infinity and NaN can also be a performance problem but it depends on your processor. Though of course if you get those during matrix multiplication something has probably gone wrong.

Depends on what you mean by wrong.

The point of NaNs and infinities is that sometimes it is easier/faster/cleaner to just let those things propagate through the calculation than to do conditional logic that in the end just tends to just simulate the same effect.

True, that consideration is more valid for element-by-element operations rather than things like matrix multiplications. But depending on the operation it might make sense to rely on NaN propagation rather than trying to prescan you matrix for NaNs.

It's just incredible that there's still no way to force denormals to be rounded silently to zero on x86. You have to go out of your way to check for them, or face horrific performance penalties for a feature you didn't even want or need.

Huh? Set bits 6 (DAZ) and 15 (FZ) in MXCSR. Done.

The only instructions you can’t set to flush are the legacy x87 opcodes, which you shouldn’t be using in a performance-sensitive context anyway.

Sadly, some of us are stuck maintaining code that still needs to run on pre-SSE2 hardware.

This code:

    #define CSR_FLUSH_TO_ZERO         (1 << 15)
    unsigned csr = __builtin_ia32_stmxcsr();
    csr |= CSR_FLUSH_TO_ZERO;
    __builtin_ia32_ldmxcsr(csr);

from https://stackoverflow.com/a/8217313/ seems to work for me with gcc. I don't know how widely supported it would be.

The quote was referring to the distribution of the data in the matrix, not the specific details of data in the matrix (such as denormals).

Sparse matrix multiplication is generally considered to be a different problem. Sparse matrices need to be represented differently (See https://en.wikipedia.org/wiki/Sparse_matrix for an introduction). There has been recently a lot of research on how to do sparse matrix multiplication efficiently on manycore architectures. See this 2018 paper for example: https://arxiv.org/pdf/1804.01698.pdf. Note that sparse matrix-vector multiplication is also a different problem.

Thats fascinating, i wonder how much coordinate matrix multiplication in spark takes into account factors like choelesky decomposition. Noticed most of those libraries are written in C++ or fortran, is JVM based approaches generally discouraged in this realm? It fits well with our current stack based mostly around spark

I don't much about Java. But, for example, if you want to implement the technique from the link above in Java, you'd need LoadFloat8, BroadcastFloat8, and AdduFloat8. I'm not sure how they are implemented in Glow, but they can be implemented using inline assembly. You'd need something like SIMD intrinsics Java library to access these crucial x86 SIMD instructions.

There's a proposal to add SIMD intrinsics to Java (sort of), but i have no idea when it might happen:

https://bugs.openjdk.java.net/browse/JDK-8201271

Not quite the same, but similar, see the flowcharts at the end of this [1].

[1]: https://www.mathworks.com/help/matlab/ref/mldivide.html

Doesn't the multiplication of two sparse matrices increase fill-in? That cannot be a good idea in general.

Absolutely, but if that’s what you need to compute, that’s what you need to compute. Sometimes it is.