<p dir="auto">I noticed that in many of my codes, seemingly harmless lines like</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="numpy.add.at(target, idx, vals)"><pre class="notranslate"><span class="pl-s1">numpy</span>.<span class="pl-s1">add</span>.<span class="pl-en">at</span>(<span class="pl-s1">target</span>, <span class="pl-s1">idx</span>, <span class="pl-s1">vals</span>)</pre></div>
<p dir="auto">take a large share of the runtime. I investigated and found that one gets a speed-up of a factor of 40 (!) by simply moving the critical code to C++.</p>
<p dir="auto"><a target="_blank" rel="noopener noreferrer nofollow" href="https://user-images.githubusercontent.com/181628/40485655-6e5a98a0-5f5f-11e8-9e5a-13d9b02eae86.png"><img src="https://user-images.githubusercontent.com/181628/40485655-6e5a98a0-5f5f-11e8-9e5a-13d9b02eae86.png" alt="f" style="max-width: 100%;"></a></p>
<p dir="auto">(Not sure what's actually done in <code class="notranslate">numpy.ufunc</code>.) I put this into a very simple module, <a href="https://github.com/nschloe/fastfunc">https://github.com/nschloe/fastfunc</a>, so feel free to take some code out of there.</p>
<p dir="auto">numpy 1.18.4 and 1.18.5 test errors</p>
<p dir="auto">run numpy.test()</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="
_______________ TestRandomDist.test_multivariate_normal[eigh] _________________

self = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40&gt;
method = 'eigh'

    @pytest.mark.parametrize(&quot;method&quot;, [&quot;svd&quot;, &quot;eigh&quot;, &quot;cholesky&quot;])
    def test_multivariate_normal(self, method):
        random = Generator(MT19937(self.seed))
        mean = (.123456789, 10)
        cov = [[1, 0], [0, 1]]
        size = (3, 2)
        actual = random.multivariate_normal(mean, cov, size, method=method)
        desired = np.array([[[-1.747478062846581,  11.25613495182354  ],
                             [-0.9967333370066214, 10.342002097029821 ]],
                            [[ 0.7850019631242964, 11.181113712443013 ],
                             [ 0.8901349653255224,  8.873825399642492 ]],
                            [[ 0.7130260107430003,  9.551628690083056 ],
                             [ 0.7127098726541128, 11.991709234143173 ]]])
    
&gt;       assert_array_almost_equal(actual, desired, decimal=15)
E       AssertionError: 
E       Arrays are not almost equal to 15 decimals
E       
E       Mismatched elements: 6 / 12 (50%)
E       Max absolute difference: 3.7418697
E       Max relative difference: 2.2477228
E        x: array([[[ 1.994391640846581, 11.25613495182354 ],
E               [ 1.243646915006621, 10.342002097029821]],
E       ...
E        y: array([[[-1.747478062846581, 11.25613495182354 ],
E               [-0.996733337006621, 10.342002097029821]],
E       ...

actual     = array([[[ 1.99439164, 11.25613495],
        [ 1.24364692, 10.3420021 ]],

       [[-0.53808839, 11.18111371],
        [-0.64322139,  8.8738254 ]],

       [[-0.46611243,  9.55162869],
        [-0.46579629, 11.99170923]]])
cov        = [[1, 0], [0, 1]]
desired    = array([[[-1.74747806, 11.25613495],
        [-0.99673334, 10.3420021 ]],

       [[ 0.78500196, 11.18111371],
        [ 0.89013497,  8.8738254 ]],

       [[ 0.71302601,  9.55162869],
        [ 0.71270987, 11.99170923]]])
mean       = (0.123456789, 10)
method     = 'eigh'
random     = Generator(MT19937) at 0x7F77D1CD1D60
self       = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40&gt;
size       = (3, 2)

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1264: AssertionError
________ TestRandomDist.test_multivariate_normal_basic_stats[cholesky] _________

self = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0&gt;
method = 'cholesky'

    @pytest.mark.parametrize(&quot;method&quot;, [&quot;svd&quot;, &quot;eigh&quot;, &quot;cholesky&quot;])
    def test_multivariate_normal_basic_stats(self, method):
        random = Generator(MT19937(self.seed))
        n_s = 1000
        mean = np.array([1, 2])
        cov = np.array([[2, 1], [1, 2]])
        s = random.multivariate_normal(mean, cov, size=(n_s,), method=method)
        s_center = s - mean
        cov_emp = (s_center.T @ s_center) / (n_s - 1)
        # these are pretty loose and are only designed to detect major errors
        assert np.all(np.abs(s_center.mean(-2)) &lt; 0.1)
&gt;       assert np.all(np.abs(cov_emp - cov) &lt; 0.2)
E       AssertionError: assert False
E        +  where False = &lt;function all at 0x7f7861f6b3a0&gt;(array([[0.55018097, 0.10366964],\n       [0.10366964, 0.51195134]]) &lt; 0.2)
E        +    where &lt;function all at 0x7f7861f6b3a0&gt; = np.all
E        +    and   array([[0.55018097, 0.10366964],\n       [0.10366964, 0.51195134]]) = &lt;ufunc 'absolute'&gt;((array([[2.55018097, 0.89633036],\n       [0.89633036, 1.48804866]]) - array([[2, 1],\n       [1, 2]])))
E        +      where &lt;ufunc 'absolute'&gt; = np.abs

cov        = array([[2, 1],
       [1, 2]])
cov_emp    = array([[2.55018097, 0.89633036],
       [0.89633036, 1.48804866]])
mean       = array([1, 2])
method     = 'cholesky'
n_s        = 1000
random     = Generator(MT19937) at 0x7F77D1CD1A90
s          = array([[-0.7576799 ,  3.53844484],
       [-0.34235607,  2.41886531],
       [ 2.77073967,  3.44656296],
       ...,
       [-0.51773839,  0.49648485],
       [ 0.57610204,  3.15188553],
       [ 1.10139309,  1.47737657]])
s_center   = array([[-1.7576799 ,  1.53844484],
       [-1.34235607,  0.41886531],
       [ 1.77073967,  1.44656296],
       ...,
       [-1.51773839, -1.50351515],
       [-0.42389796,  1.15188553],
       [ 0.10139309, -0.52262343]])
self       = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0&gt;

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1335: AssertionError
=========================== short test summary info ============================
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal[eigh]
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal_basic_stats[cholesky]
2 failed, 10264 passed, 68 skipped, 180 deselected, 17 xfailed, 3 xpassed in 69.53s (0:01:09)
False"><pre class="notranslate"><code class="notranslate">
_______________ TestRandomDist.test_multivariate_normal[eigh] _________________

self = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40&gt;
method = 'eigh'

    @pytest.mark.parametrize("method", ["svd", "eigh", "cholesky"])
    def test_multivariate_normal(self, method):
        random = Generator(MT19937(self.seed))
        mean = (.123456789, 10)
        cov = [[1, 0], [0, 1]]
        size = (3, 2)
        actual = random.multivariate_normal(mean, cov, size, method=method)
        desired = np.array([[[-1.747478062846581,  11.25613495182354  ],
                             [-0.9967333370066214, 10.342002097029821 ]],
                            [[ 0.7850019631242964, 11.181113712443013 ],
                             [ 0.8901349653255224,  8.873825399642492 ]],
                            [[ 0.7130260107430003,  9.551628690083056 ],
                             [ 0.7127098726541128, 11.991709234143173 ]]])
    
&gt;       assert_array_almost_equal(actual, desired, decimal=15)
E       AssertionError: 
E       Arrays are not almost equal to 15 decimals
E       
E       Mismatched elements: 6 / 12 (50%)
E       Max absolute difference: 3.7418697
E       Max relative difference: 2.2477228
E        x: array([[[ 1.994391640846581, 11.25613495182354 ],
E               [ 1.243646915006621, 10.342002097029821]],
E       ...
E        y: array([[[-1.747478062846581, 11.25613495182354 ],
E               [-0.996733337006621, 10.342002097029821]],
E       ...

actual     = array([[[ 1.99439164, 11.25613495],
        [ 1.24364692, 10.3420021 ]],

       [[-0.53808839, 11.18111371],
        [-0.64322139,  8.8738254 ]],

       [[-0.46611243,  9.55162869],
        [-0.46579629, 11.99170923]]])
cov        = [[1, 0], [0, 1]]
desired    = array([[[-1.74747806, 11.25613495],
        [-0.99673334, 10.3420021 ]],

       [[ 0.78500196, 11.18111371],
        [ 0.89013497,  8.8738254 ]],

       [[ 0.71302601,  9.55162869],
        [ 0.71270987, 11.99170923]]])
mean       = (0.123456789, 10)
method     = 'eigh'
random     = Generator(MT19937) at 0x7F77D1CD1D60
self       = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40&gt;
size       = (3, 2)

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1264: AssertionError
________ TestRandomDist.test_multivariate_normal_basic_stats[cholesky] _________

self = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0&gt;
method = 'cholesky'

    @pytest.mark.parametrize("method", ["svd", "eigh", "cholesky"])
    def test_multivariate_normal_basic_stats(self, method):
        random = Generator(MT19937(self.seed))
        n_s = 1000
        mean = np.array([1, 2])
        cov = np.array([[2, 1], [1, 2]])
        s = random.multivariate_normal(mean, cov, size=(n_s,), method=method)
        s_center = s - mean
        cov_emp = (s_center.T @ s_center) / (n_s - 1)
        # these are pretty loose and are only designed to detect major errors
        assert np.all(np.abs(s_center.mean(-2)) &lt; 0.1)
&gt;       assert np.all(np.abs(cov_emp - cov) &lt; 0.2)
E       AssertionError: assert False
E        +  where False = &lt;function all at 0x7f7861f6b3a0&gt;(array([[0.55018097, 0.10366964],\n       [0.10366964, 0.51195134]]) &lt; 0.2)
E        +    where &lt;function all at 0x7f7861f6b3a0&gt; = np.all
E        +    and   array([[0.55018097, 0.10366964],\n       [0.10366964, 0.51195134]]) = &lt;ufunc 'absolute'&gt;((array([[2.55018097, 0.89633036],\n       [0.89633036, 1.48804866]]) - array([[2, 1],\n       [1, 2]])))
E        +      where &lt;ufunc 'absolute'&gt; = np.abs

cov        = array([[2, 1],
       [1, 2]])
cov_emp    = array([[2.55018097, 0.89633036],
       [0.89633036, 1.48804866]])
mean       = array([1, 2])
method     = 'cholesky'
n_s        = 1000
random     = Generator(MT19937) at 0x7F77D1CD1A90
s          = array([[-0.7576799 ,  3.53844484],
       [-0.34235607,  2.41886531],
       [ 2.77073967,  3.44656296],
       ...,
       [-0.51773839,  0.49648485],
       [ 0.57610204,  3.15188553],
       [ 1.10139309,  1.47737657]])
s_center   = array([[-1.7576799 ,  1.53844484],
       [-1.34235607,  0.41886531],
       [ 1.77073967,  1.44656296],
       ...,
       [-1.51773839, -1.50351515],
       [-0.42389796,  1.15188553],
       [ 0.10139309, -0.52262343]])
self       = &lt;numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0&gt;

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1335: AssertionError
=========================== short test summary info ============================
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal[eigh]
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal_basic_stats[cholesky]
2 failed, 10264 passed, 68 skipped, 180 deselected, 17 xfailed, 3 xpassed in 69.53s (0:01:09)
False
</code></pre></div>
<p dir="auto">OS is ubuntu, same failure on python-3.7.6 and 3.8.2.<br>
numpy installed via "pip3 install  --no-binary numy numpy" compiled from source with MKL.<br>
But same test failures persist if link numpy to OpenBLAS (version 0.3.9)<br>
numpy-1.17.5 (installed from source the same way) produces no test error in both.</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; np.show_config()
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; np.show_config()
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/opt/intel/mkl/lib/intel64/']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/intel/mkl/include']

</code></pre></div>
<p dir="auto">One of the failures is the same as <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="562885885" data-permission-text="Title is private" data-url="https://github.com/numpy/numpy/issues/15546" data-hovercard-type="issue" data-hovercard-url="/numpy/numpy/issues/15546/hovercard" href="https://github.com/numpy/numpy/issues/15546">#15546</a> for numpy-1.18.1<br>
but this time it fails also in python 3.7.6</p>