Multi-threaded integer matrix multiplication in NumPy/SciPy
<p>Doing something like</p>

<pre><code>import numpy as np
a = np.random.rand(10**4, 10**4)
b = np.dot(a, a)
</code></pre>

<p>uses multiple cores, and it runs nicely.</p>

<p>The elements in <code>a</code>, though, are 64-bit floats (or 32-bit in 32-bit platforms?), and I'd like to multiply 8-bit integer arrays. Trying the following, though:</p>

<pre><code>a = np.random.randint(2, size=(n, n)).astype(np.int8)
</code></pre>

<p>results in the dot product not using multiple cores, and thus running ~1000x slower on my PC.</p>

<pre><code>array: np.random.randint(2, size=shape).astype(dtype)

dtype    shape          %time (average)

float32 (2000, 2000)    62.5 ms
float32 (3000, 3000)    219 ms
float32 (4000, 4000)    328 ms
float32 (10000, 10000)  4.09 s

int8    (2000, 2000)    13 seconds
int8    (3000, 3000)    3min 26s
int8    (4000, 4000)    12min 20s
int8    (10000, 10000)  It didn't finish in 6 hours

float16 (2000, 2000)    2min 25s
float16 (3000, 3000)    Not tested
float16 (4000, 4000)    Not tested
float16 (10000, 10000)  Not tested
</code></pre>

<p>I understand NumPy uses BLAS, which doesn't support integers, but if I use the SciPy BLAS wrappers, ie.</p>

<pre><code>import scipy.linalg.blas as blas
a = np.random.randint(2, size=(n, n)).astype(np.int8)
b = blas.sgemm(alpha=1.0, a=a, b=a)
</code></pre>

<p>the computation <em>is</em> multi-threaded. Now, <code>blas.sgemm</code> runs with exactly the same timing as <code>np.dot</code> for float32's, but for non-floats it converts everything to <code>float32</code> and outputs floats, which is something <code>np.dot</code> doesn't do. (In addition, <code>b</code> is now in <code>F_CONTIGUOUS</code> order, which is a lesser issue). </p>

<p>So, if I want to do integer matrix multiplication, I have to do one of the following:</p>

<ol>
<li>Use NumPy's painfully slow <code>np.dot</code> and be glad I get to keep the 8-bit integers.</li>
<li>Use SciPy's <code>sgemm</code> and use up 4x memory.</li>
<li>Use Numpy's <code>np.float16</code> and only use up 2x memory, with the caveat that <code>np.dot</code> is much slower on float16 arrays than on float32 arrays, more so than int8.</li>
<li>Find an optimized library for multi-threaded integer matrix multiplication (actually, <strong>Mathematica</strong> does this, but I'd prefer a Python solution), ideally supporting 1-bit arrays, although 8-bit arrays is also fine... (I'm actually aiming to do multiplication of matrices over the finite field Z/2Z, and I know I can do this with <strong>Sage</strong>, which is quite Pythonic, but, again, is there something strictly Python?)</li>
</ol>

<p>Can I follow option 4? Does such a library exist?</p>

<p>Disclaimer: I'm actually running NumPy + MKL, but I've tried a similar test on vanilly NumPy, with similar results.</p>