Why MATLAB/Numpy/Scipy performance is slow and doesn't reach CPU capabilities (flops)?
<p>First of all, I know there are multiple threads that touch this issue, however I could not get a straight answer and encountered some flops miscalculations.</p>

<p>I have prepared a MATLAB and Python benchmark of element-wise multiplication. This is the simplest and most forward way one can easily calculate the flop count. </p>

<p>It uses NxN array (matrix) but does not do a matrix multiplication, rather an element wise multiplication. This is important because when using the matrix multiplication the number of operation is not N^3 !!!</p>

<p>The <a href="https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm" rel="nofollow">lower level algorithm that performs matrix multiplication</a> does it in less than N^3 operations.</p>

<p>The execution of element-wise multiplication of randomly generated numbers however, will have to be performed in N^2 operations</p>

<p>I have an intel i7-4770 (which I think have 4 physical cores and 8 virtual cores) @ 3.5GHz. So if assuming 4flops per cycle that should be 14 GFLOPS per core!</p>

<p>MATLAB/Numpy/Scipy don't get anywhere near.</p>

<p>Why?</p>

<p>MATLAB:</p>

<pre><code>%element wise multiplication benchmark
N = 10^4;
nOps = N^2;
m1 = randn(N);
m2 = randn(size(m1));
m = randn(size(m1));

m1 = single(m1);
m2 = single(m2);

% clear m
tic
m1 = m1 .* m2;
t = toc;
gflops = nOps/t*1e-9;
t_gflops = [t gflops]

% clear m
tic
m1 = m1.*m2;
t = toc;
gflops = nOps/t*1e-9;
t_gflops = [t gflops]


% clear m
tic
m1 = m1.*m2;
t = toc;
gflops = nOps/t*1e-9;
t_gflops = [t gflops]

version('-blas')
version('-lapack')
</code></pre>

<p>Which results in:</p>

<pre><code>t_gflops =

    0.0978    1.0226


t_gflops =

    0.0743    1.3458


t_gflops =

    0.0731    1.3682


ans =

Intel(R) Math Kernel Library Version 11.1.1 Product Build 20131010 for Intel(R) 64 architecture applications



ans =

Intel(R) Math Kernel Library Version 11.1.1 Product Build 20131010 for Intel(R) 64 architecture applications
Linear Algebra PACKage Version 3.4.1
</code></pre>

<p>And now Python:</p>

<pre><code>import numpy as np
# import gnumpy as gnp
import scipy as sp
import scipy.linalg as la
import time

if __name__ == '__main__':
    N = 10**4
    nOps = N**2
    a = np.random.randn(N,N).astype(np.float32)
    b = np.random.randn(N,N).astype(np.float32)

    t = time.time()
    c = a*b
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = np.multiply(a, b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = sp.multiply(a, b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = sp.multiply(a, b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    a = np.random.randn(N,1).astype(np.float32)
    b = np.random.randn(1,N).astype(np.float32)

    t = time.time()
    c1 = np.dot(a, b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = np.dot(a, b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = la.blas.dgemm(1.0,a,b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    t = time.time()
    c = la._fblas.dgemm(1.0,a,b)
    dt = time.time()-t
    gflops = nOps/dt*1e-9
    print("dt = ", dt, ", gflops = ", gflops)

    print("numpy config")
    np.show_config()
    print("scipy config")
    sp.show_config()
    # numpy
</code></pre>

<p>Which result in:</p>

<pre><code>dt =  0.16301608085632324 , gflops =  0.6134364136022663
dt =  0.16701674461364746 , gflops =  0.5987423610209003
dt =  0.1770176887512207 , gflops =  0.5649152957845881
dt =  0.188018798828125 , gflops =  0.5318617107612401
dt =  0.151015043258667 , gflops =  0.6621856858903415
dt =  0.17201733589172363 , gflops =  0.5813367558659613
dt =  0.3080308437347412 , gflops =  0.3246428142959423
dt =  0.39503931999206543 , gflops =  0.253139358385916
</code></pre>

<p>numpy config</p>

<p>mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library\\include']
</code></pre>

<p>lapack_mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library\\include']
</code></pre>

<p>lapack_opt_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library\\include']
</code></pre>

<p>blas_opt_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library\\include']
</code></pre>

<p>openblas_lapack_info:</p>

<p>NOT AVAILABLE
blas_mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library\\include']
</code></pre>

<p>scipy config</p>

<p>mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library', 'C:\\Minonda\\envs\\_build\\Library\\include', 'C:\\Minonda\\envs\\_build\\Library\\lib']
</code></pre>

<p>lapack_mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library', 'C:\\Minonda\\envs\\_build\\Library\\include', 'C:\\Minonda\\envs\\_build\\Library\\lib']
</code></pre>

<p>lapack_opt_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_lapack95_lp64', 'mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library', 'C:\\Minonda\\envs\\_build\\Library\\include', 'C:\\Minonda\\envs\\_build\\Library\\lib']
</code></pre>

<p>blas_opt_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library', 'C:\\Minonda\\envs\\_build\\Library\\include', 'C:\\Minonda\\envs\\_build\\Library\\lib']
</code></pre>

<p>openblas_lapack_info:</p>

<p>NOT AVAILABLE
blas_mkl_info:</p>

<pre><code>define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
libraries = ['mkl_core_dll', 'mkl_intel_lp64_dll', 'mkl_intel_thread_dll']
library_dirs = ['C:\\Minonda\\envs\\_build\\Library\\lib']
include_dirs = ['C:\\Minonda\\envs\\_build\\Library', 'C:\\Minonda\\envs\\_build\\Library\\include', 'C:\\Minonda\\envs\\_build\\Library\\lib']
</code></pre>

<p>Process finished with exit code 0</p>

<p>Well, in this case you're being limited by memory bandwidth, not CPU power. Assuming:</p>

<ul>
<li><a href="https://en.wikipedia.org/wiki/DDR3_SDRAM#JEDEC_standard_modules" rel="nofollow">PC3-12800
RAM</a>
in dual channel mode; </li>
<li>for every multiplication (single precision) 12 bytes need to be transferred between CPU and RAM;</li>
</ul>

<p>The theoretical maximum sustained performance would be about 2 GFLOPS. I calculated this number as <code>peak DDR3 transfer rate</code> * <code>number of RAM channels</code> / <code>bytes tranferred per FLOP</code>.</p>

<p>BTW, in numpy elementwise operations are not accelerated by BLAS. I'm not sure about MATAB.</p>

python performance matlab numpy scipy