<p dir="auto">I have a freshly-built numpy on OS X 10.9, and am running into pretty regular crashes related to umath.so. <a href="https://gist.github.com/fonnesbeck/31d372805647da7af6e3">Here is a Gist with the crash report</a>. It seems to crash in a variety of ways, but I can reliably get it to crash by importing a module like Theano, which has numpy as a dependency. I have re-built all of my packages since upgrading to OSX 10.9.</p>
<h3 dir="auto">Describe the issue:</h3>
<p dir="auto">I have replaced the last 3 lines in this code</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="A = np.ones((150000,3,3))
j = np.ones((150000))
s = np.ones((150000))
jDuplicated = np.repeat(j, 3**2).reshape((len(elements), dim, dim))
sDuplicated = np.repeat(s, 3**2).reshape((len(elements), dim, dim))
res= sDuplicated * A @ np.swapaxes(A,1,2) * jDuplicated"><pre class="notranslate"><span class="pl-v">A</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>,<span class="pl-c1">3</span>,<span class="pl-c1">3</span>))
<span class="pl-s1">j</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>))
<span class="pl-s1">s</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>))
<span class="pl-s1">jDuplicated</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">repeat</span>(<span class="pl-s1">j</span>, <span class="pl-c1">3</span><span class="pl-c1">**</span><span class="pl-c1">2</span>).<span class="pl-en">reshape</span>((<span class="pl-en">len</span>(<span class="pl-s1">elements</span>), <span class="pl-s1">dim</span>, <span class="pl-s1">dim</span>))
<span class="pl-s1">sDuplicated</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">repeat</span>(<span class="pl-s1">s</span>, <span class="pl-c1">3</span><span class="pl-c1">**</span><span class="pl-c1">2</span>).<span class="pl-en">reshape</span>((<span class="pl-en">len</span>(<span class="pl-s1">elements</span>), <span class="pl-s1">dim</span>, <span class="pl-s1">dim</span>))
<span class="pl-s1">res</span><span class="pl-c1">=</span> <span class="pl-s1">sDuplicated</span> <span class="pl-c1">*</span> <span class="pl-v">A</span> @ <span class="pl-s1">np</span>.<span class="pl-en">swapaxes</span>(<span class="pl-v">A</span>,<span class="pl-c1">1</span>,<span class="pl-c1">2</span>) <span class="pl-c1">*</span> <span class="pl-s1">jDuplicated</span></pre></div>
<p dir="auto">with</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="res = np.einsum('i,i,ijk,ilk-&gt;ijl',s,j,A,A)"><pre class="notranslate"><span class="pl-s1">res</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">einsum</span>(<span class="pl-s">'i,i,ijk,ilk-&gt;ijl'</span>,<span class="pl-s1">s</span>,<span class="pl-s1">j</span>,<span class="pl-v">A</span>,<span class="pl-v">A</span>)</pre></div>
<p dir="auto">The single line looks much better than duplicating a lot of data so that I can perform element-wise multiplication. However the top version is around 10% faster on my pc.</p>
<p dir="auto">I would expect einsum to be faster, since it does not replicate a lot of data. Does anybody have an idea why this is not the case? Maybe the repeated data causes less cache misses?</p>
<h3 dir="auto">Reproduce the code example:</h3>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np

A = np.ones((150000,3,3))
j = np.ones((150000))
s = np.ones((150000))
# the numbers might need to be random to avoid simplification
res = np.einsum('i,i,ijk,ilk-&gt;ijl',s,j,A,A)"><pre class="notranslate"><span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>

<span class="pl-v">A</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>,<span class="pl-c1">3</span>,<span class="pl-c1">3</span>))
<span class="pl-s1">j</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>))
<span class="pl-s1">s</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">150000</span>))
<span class="pl-c"># the numbers might need to be random to avoid simplification</span>
<span class="pl-s1">res</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">einsum</span>(<span class="pl-s">'i,i,ijk,ilk-&gt;ijl'</span>,<span class="pl-s1">s</span>,<span class="pl-s1">j</span>,<span class="pl-v">A</span>,<span class="pl-v">A</span>)</pre></div>
<h3 dir="auto">Error message:</h3>
<p dir="auto"><em>No response</em></p>
<h3 dir="auto">NumPy/Python version information:</h3>
<p dir="auto">1.22.1 3.10.2 (main, Jan 15 2022, 18:03:19) [GCC 7.5.0]</p>