<p dir="auto"><code class="notranslate">matmul</code> (aka <code class="notranslate">@</code>) is extremely slow when one of the two arrays is the real or imaginary part of a complex array extracted with <code class="notranslate">.real</code> or <code class="notranslate">.imag</code>.</p>
<p dir="auto">This was already observed in <a href="https://github.com/numpy/numpy/issues/13229" data-hovercard-type="issue" data-hovercard-url="/numpy/numpy/issues/13229/hovercard">this issue</a>, but was not marked as a bug, while I think this should be.</p>
<h3 dir="auto">Reproducing code example:</h3>
<p dir="auto">Consider the following benchmark.</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np
from timeit import repeat

A  = np.random.rand(200, 1000) + 0j
B = np.random.rand(1000, 50000)

def time(stmt):
    print(f'testing {stmt}')
    res = repeat(stmt, globals={'A':A,'B':B, 'np':np}, number=1, repeat=10)
    res = np.array(res)*1000
    print(f'{np.mean(res):.0f} ms ± {np.std(res):.0f} ms per loop (mean ± std. dev. of 10 runs, 1 loop each)')

time('A @ B')
time('A.real @ B')
time('A.imag @ B')
time('A.astype(np.float64()) @ B')
time('A.real.astype(np.float64()) @ B')
time('A.real.copy() @ B')"><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-k">from</span> <span class="pl-s1">timeit</span> <span class="pl-k">import</span> <span class="pl-s1">repeat</span>

<span class="pl-v">A</span>  <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">random</span>.<span class="pl-en">rand</span>(<span class="pl-c1">200</span>, <span class="pl-c1">1000</span>) <span class="pl-c1">+</span> <span class="pl-c1">0j</span>
<span class="pl-v">B</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">random</span>.<span class="pl-en">rand</span>(<span class="pl-c1">1000</span>, <span class="pl-c1">50000</span>)

<span class="pl-k">def</span> <span class="pl-en">time</span>(<span class="pl-s1">stmt</span>):
    <span class="pl-en">print</span>(<span class="pl-s">f'testing <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">stmt</span><span class="pl-kos">}</span></span>'</span>)
    <span class="pl-s1">res</span> <span class="pl-c1">=</span> <span class="pl-en">repeat</span>(<span class="pl-s1">stmt</span>, <span class="pl-s1">globals</span><span class="pl-c1">=</span>{<span class="pl-s">'A'</span>:<span class="pl-v">A</span>,<span class="pl-s">'B'</span>:<span class="pl-v">B</span>, <span class="pl-s">'np'</span>:<span class="pl-s1">np</span>}, <span class="pl-s1">number</span><span class="pl-c1">=</span><span class="pl-c1">1</span>, <span class="pl-s1">repeat</span><span class="pl-c1">=</span><span class="pl-c1">10</span>)
    <span class="pl-s1">res</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">array</span>(<span class="pl-s1">res</span>)<span class="pl-c1">*</span><span class="pl-c1">1000</span>
    <span class="pl-en">print</span>(<span class="pl-s">f'<span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">np</span>.<span class="pl-en">mean</span>(<span class="pl-s1">res</span>):.0f<span class="pl-kos">}</span></span> ms ± <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">np</span>.<span class="pl-en">std</span>(<span class="pl-s1">res</span>):.0f<span class="pl-kos">}</span></span> ms per loop (mean ± std. dev. of 10 runs, 1 loop each)'</span>)

<span class="pl-en">time</span>(<span class="pl-s">'A @ B'</span>)
<span class="pl-en">time</span>(<span class="pl-s">'A.real @ B'</span>)
<span class="pl-en">time</span>(<span class="pl-s">'A.imag @ B'</span>)
<span class="pl-en">time</span>(<span class="pl-s">'A.astype(np.float64()) @ B'</span>)
<span class="pl-en">time</span>(<span class="pl-s">'A.real.astype(np.float64()) @ B'</span>)
<span class="pl-en">time</span>(<span class="pl-s">'A.real.copy() @ B'</span>)</pre></div>
<p dir="auto">On my PC the results are:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="testing A @ B
1046 ms ± 73 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real @ B
21980 ms ± 4978 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.imag @ B
18140 ms ± 970 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.astype(np.float64()) @ B
&lt;string&gt;:6: ComplexWarning: Casting complex values to real discards the imaginary part
161 ms ± 16 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real.astype(np.float64()) @ B
166 ms ± 12 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real.copy() @ B
150 ms ± 11 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)"><pre class="notranslate"><code class="notranslate">testing A @ B
1046 ms ± 73 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real @ B
21980 ms ± 4978 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.imag @ B
18140 ms ± 970 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.astype(np.float64()) @ B
&lt;string&gt;:6: ComplexWarning: Casting complex values to real discards the imaginary part
161 ms ± 16 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real.astype(np.float64()) @ B
166 ms ± 12 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

testing A.real.copy() @ B
150 ms ± 11 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)
</code></pre></div>
<p dir="auto">Note how the most trivial/elegant/obvious solution <code class="notranslate">A.real @ B</code> (and <code class="notranslate">A.imag @ B</code> depending on what you have to do) is 120 times slower than <code class="notranslate">A.real.copy() @ B</code></p>
<p dir="auto">Also <code class="notranslate">A@B</code> is much slower, but I suspect this has to be with the fact that in this case <code class="notranslate">A</code> is complex bringing <code class="notranslate">matmul</code> on a different path, even if <code class="notranslate">A.imag</code> is all 0.</p>
<h3 dir="auto">NumPy/Python version information:</h3>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import sys, numpy; print(numpy.__version__, sys.version)
1.18.1 3.7.4 (default, Aug  9 2019, 18:34:13) [MSC v.1915 64 bit (AMD64)]"><pre class="notranslate"><code class="notranslate">import sys, numpy; print(numpy.__version__, sys.version)
1.18.1 3.7.4 (default, Aug  9 2019, 18:34:13) [MSC v.1915 64 bit (AMD64)]
</code></pre></div>
<h3 dir="auto">Describe the issue:</h3>
<p dir="auto">Hi there,</p>
<p dir="auto">Let me go straight to the point. I have a list of arrays which I want to apply np.argsort() to. I was finding wrong answers using the <code class="notranslate">axis</code> option, so I tried to debug it with a simple loop, printing everything.</p>
<p dir="auto">The results clearly don't match expectations in some rows, and frankly it's very weird. I am pasting the code (and output in a comment) below.</p>
<p dir="auto">Thanks in advance.</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="mat = np.array([[0.00514059, 0.00506887, 0.00574913, 0.00531327, 0.00518162],
       [0.00494017, 0.00489627, 0.00566709, 0.00454714, 0.00499841],
       [0.00493342, 0.00488984, 0.00564926, 0.00505881, 0.00494269],
       [0.00491196, 0.00488021, 0.00528805, 0.00511658, 0.00481696],
       [0.00487965, 0.00483828, 0.00538671, 0.004878  , 0.00493741]])

for m in mat:
    print(m)
    print(np.argsort(m)

'''
[0.00514059 0.00506887 0.00574913 0.00531327 0.00518162]
[1 0 4 3 2]
[0.00494017 0.00489627 0.00566709 0.00454714 0.00499841]
[3 1 0 4 2]
[0.00493342 0.00488984 0.00564926 0.00505881 0.00494269]
[1 0 4 3 2]
[0.00491196 0.00488021 0.00528805 0.00511658 0.00481696]
[4 1 0 3 2]
[0.00487965 0.00483828 0.00538671 0.004878   0.00493741]
[1 3 0 4 2]
'''"><pre class="notranslate"><span class="pl-s1">mat</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">array</span>([[<span class="pl-c1">0.00514059</span>, <span class="pl-c1">0.00506887</span>, <span class="pl-c1">0.00574913</span>, <span class="pl-c1">0.00531327</span>, <span class="pl-c1">0.00518162</span>],
       [<span class="pl-c1">0.00494017</span>, <span class="pl-c1">0.00489627</span>, <span class="pl-c1">0.00566709</span>, <span class="pl-c1">0.00454714</span>, <span class="pl-c1">0.00499841</span>],
       [<span class="pl-c1">0.00493342</span>, <span class="pl-c1">0.00488984</span>, <span class="pl-c1">0.00564926</span>, <span class="pl-c1">0.00505881</span>, <span class="pl-c1">0.00494269</span>],
       [<span class="pl-c1">0.00491196</span>, <span class="pl-c1">0.00488021</span>, <span class="pl-c1">0.00528805</span>, <span class="pl-c1">0.00511658</span>, <span class="pl-c1">0.00481696</span>],
       [<span class="pl-c1">0.00487965</span>, <span class="pl-c1">0.00483828</span>, <span class="pl-c1">0.00538671</span>, <span class="pl-c1">0.004878</span>  , <span class="pl-c1">0.00493741</span>]])

<span class="pl-k">for</span> <span class="pl-s1">m</span> <span class="pl-c1">in</span> <span class="pl-s1">mat</span>:
    <span class="pl-en">print</span>(<span class="pl-s1">m</span>)
    <span class="pl-s1">print</span>(<span class="pl-s1">np</span>.<span class="pl-en">argsort</span>(<span class="pl-s1">m</span>)

<span class="pl-s">'''</span>
<span class="pl-s">[0.00514059 0.00506887 0.00574913 0.00531327 0.00518162]</span>
<span class="pl-s">[1 0 4 3 2]</span>
<span class="pl-s">[0.00494017 0.00489627 0.00566709 0.00454714 0.00499841]</span>
<span class="pl-s">[3 1 0 4 2]</span>
<span class="pl-s">[0.00493342 0.00488984 0.00564926 0.00505881 0.00494269]</span>
<span class="pl-s">[1 0 4 3 2]</span>
<span class="pl-s">[0.00491196 0.00488021 0.00528805 0.00511658 0.00481696]</span>
<span class="pl-s">[4 1 0 3 2]</span>
<span class="pl-s">[0.00487965 0.00483828 0.00538671 0.004878   0.00493741]</span>
<span class="pl-s">[1 3 0 4 2]</span>
<span class="pl-s">'''</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.20.3 3.9.7 (default, Aug 31 2021, 13:28:12)<br>
[GCC 11.1.0]</p>