6 ms·

How would this change if you were dealing with more sparse matrices? Im used to using spark mllib distributed to do matrix multiplication on large sparse coord

How would this change if you were dealing with more sparse matrices? Im used to using spark mllib distributed to do matrix multiplication on large sparse coordinate matrices and it works fine for me but i dont think i need the level of optimality this approach would demand.

knlji 8y ago
Doesn't the multiplication of two sparse matrices increase fill-in? That cannot be a good idea in general.
- stephencanon 8y ago
  Absolutely, but if that’s what you need to compute, that’s what you need to compute. Sometimes it is.
hadibrais 8y ago
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 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 https://arxiv.org/pdf/1804.01698.pdf. Note that sparse matrix-vector multiplication is also a different problem.
- nautilus12 8y ago
  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
  - joshuamorton 8y ago
    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 https://www.mathworks.com/help/matlab/ref/mldivide.html
  - hadibrais 8y ago
    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.
    - twic 8y ago
      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 https://bugs.openjdk.java.net/browse/JDK-8201271
yvdriess 8y ago
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.
- imurray 8y ago
  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_issues https://en.wikipedia.org/wiki/Denormal_number#Performance_is...

knlji 8y ago

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

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

hadibrais 8y ago

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 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 https://arxiv.org/pdf/1804.01698.pdf. Note that sparse matrix-vector multiplication is also a different problem.

nautilus12 8y ago
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
- joshuamorton 8y ago
  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 https://www.mathworks.com/help/matlab/ref/mldivide.html
- hadibrais 8y ago
  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.
  - twic 8y ago
    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 https://bugs.openjdk.java.net/browse/JDK-8201271

yvdriess 8y ago

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.

imurray 8y ago
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_issues https://en.wikipedia.org/wiki/Denormal_number#Performance_is...