Does scipy support multithreading for sparse matrix multiplication when using MKL BLAS?

Question

According to MKL BLAS documentation "All matrix-matrix operations (level 3) are threaded for both dense and sparse BLAS." http://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library

I have built Scipy with MKL BLAS. Using the test code below, I see the expected multithreaded speedup for dense, but not sparse, matrix multiplication. Are there any changes to Scipy to enable multithreaded sparse operations?

# test dense matrix multiplication
from numpy import *
import time    
x = random.random((10000,10000))
t1 = time.time()
foo = dot(x.T, x)
print time.time() - t1

# test sparse matrix multiplication
from scipy import sparse
x = sparse.rand(10000,10000)
t1 = time.time()
foo = dot(x.T, x)
print time.time() - t1

Countable

Answer 1

As far as I know, the answer is no. But, you can build your own wrapper around the MKL sparse multiply routines. You asked about the multiplying two sparse matrices. Below is some a wrapper code I used for multiplying one sparse matrix times a dense vector, so it shouldn't be hard to adapt (look at the Intel MKL reference for mkl_cspblas_dcsrgemm). Also, be aware of how your scipy arrays are stored: default is coo, but csr (or csc) may be better choices. I chose csr, but MKL supports most types (just call the appropriate routine).

From what I could tell, both scipy's default and MKL are multithreaded. By changing OMP_NUM_THREADS I could see a difference in performance.

To use the function below, if you havea a recent version of MKL, just make sure you have LD_LIBRARY_PATHS set to include the relevant MKL directories. For older versions, you need to build some specific libraries. I got my information from IntelMKL in python

def SpMV_viaMKL( A, x ):
 """
 Wrapper to Intel's SpMV
 (Sparse Matrix-Vector multiply)
 For medium-sized matrices, this is 4x faster
 than scipy's default implementation
 Stephen Becker, April 24 2014
 [email protected]
 """

 import numpy as np
 import scipy.sparse as sparse
 from ctypes import POINTER,c_void_p,c_int,c_char,c_double,byref,cdll
 mkl = cdll.LoadLibrary("libmkl_rt.so")

 SpMV = mkl.mkl_cspblas_dcsrgemv
 # Dissecting the "cspblas_dcsrgemv" name:
 # "c" - for "c-blas" like interface (as opposed to fortran)
 #    Also means expects sparse arrays to use 0-based indexing, which python does
 # "sp"  for sparse
 # "d"   for double-precision
 # "csr" for compressed row format
 # "ge"  for "general", e.g., the matrix has no special structure such as symmetry
 # "mv"  for "matrix-vector" multiply

 if not sparse.isspmatrix_csr(A):
     raise Exception("Matrix must be in csr format")
 (m,n) = A.shape

 # The data of the matrix
 data    = A.data.ctypes.data_as(POINTER(c_double))
 indptr  = A.indptr.ctypes.data_as(POINTER(c_int))
 indices = A.indices.ctypes.data_as(POINTER(c_int))

 # Allocate output, using same conventions as input
 nVectors = 1
 if x.ndim is 1:
    y = np.empty(m,dtype=np.double,order='F')
    if x.size != n:
        raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n))
 elif x.shape[1] is 1:
    y = np.empty((m,1),dtype=np.double,order='F')
    if x.shape[0] != n:
        raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n))
 else:
    nVectors = x.shape[1]
    y = np.empty((m,nVectors),dtype=np.double,order='F')
    if x.shape[0] != n:
        raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n))

 # Check input
 if x.dtype.type is not np.double:
    x = x.astype(np.double,copy=True)
 # Put it in column-major order, otherwise for nVectors > 1 this FAILS completely
 if x.flags['F_CONTIGUOUS'] is not True:
    x = x.copy(order='F')

 if nVectors == 1:
    np_x = x.ctypes.data_as(POINTER(c_double))
    np_y = y.ctypes.data_as(POINTER(c_double))
    # now call MKL. This returns the answer in np_y, which links to y
    SpMV(byref(c_char("N")), byref(c_int(m)),data ,indptr, indices, np_x, np_y ) 
 else:
    for columns in xrange(nVectors):
        xx = x[:,columns]
        yy = y[:,columns]
        np_x = xx.ctypes.data_as(POINTER(c_double))
        np_y = yy.ctypes.data_as(POINTER(c_double))
        SpMV(byref(c_char("N")), byref(c_int(m)),data,indptr, indices, np_x, np_y ) 

 return y

Engrossment

Thanks for this beautiful answer! I had to change c_char("N") to c_char(b'N')) (a byte string), but otherwise this worked very well out of the box. – Flora

Answer 2

(Idea for answer copied from https://mcmap.net/q/1773152/-is-it-possible-to-use-blas-to-speed-up-sparse-matrix-multiplication)

There is a Python wrapper that allows for multithreaded sparse matrix operations with the intel MKL library (project website: https://pypi.org/project/sparse-dot-mkl/). It can be installed using conda as well (follow instructions on https://anaconda.org/conda-forge/sparse_dot_mkl)!

Using this package, your code would change to

# test dense matrix multiplication
from numpy import *
import time    
x = random.random((10000, 10000))
t1 = time.time()
foo = dot(x.T, x)
print(time.time() - t1)

# test sparse matrix multiplication
from scipy import sparse
from sparse_dot_mkl import dot_product_mkl
x = sparse.rand(10000, 10000, format="csr")
t1 = time.time()
foo = dot_product_mkl(x.T, x)
print(time.time() - t1)

Note that you need to specify the format of the sparse matrix to either CSR or CSC to use dot_product_mkl. I changed the print statements to Python 3 syntax as well. The timings on my machine were 3.39s for the numpy code and 0.27s for the sparse multiplication.

Loisloise

While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review – Candi

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Krilov