<p dir="auto">My issue is about high memory usage when using orthogonal indexing on sparse arrays, dependent on the shape of the indexers.</p>
<h4 dir="auto">Reproducing code example:</h4>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="%load_ext memory_profiler
from scipy import sparse
import numpy as np

s = sparse.random(10000, 10000, format=&quot;csr&quot;)
idx = np.arange(10000)

# This uses a reasonable amount of memory
%memit s[idx.reshape(-1, 1), idx]
# peak memory: 135.93 MiB, increment: 46.52 MiB

# This does not
%memit s[np.ix_(idx, idx)]
# peak memory: 2066.09 MiB, increment: 1930.15 MiB"><pre class="notranslate"><span class="pl-c1">%</span><span class="pl-s1">load_ext</span> <span class="pl-s1">memory_profiler</span>
<span class="pl-k">from</span> <span class="pl-s1">scipy</span> <span class="pl-k">import</span> <span class="pl-s1">sparse</span>
<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-s1">s</span> <span class="pl-c1">=</span> <span class="pl-s1">sparse</span>.<span class="pl-en">random</span>(<span class="pl-c1">10000</span>, <span class="pl-c1">10000</span>, <span class="pl-s1">format</span><span class="pl-c1">=</span><span class="pl-s">"csr"</span>)
<span class="pl-s1">idx</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">10000</span>)

<span class="pl-c"># This uses a reasonable amount of memory</span>
<span class="pl-c1">%</span><span class="pl-s1">memit</span> <span class="pl-s1">s</span>[<span class="pl-s1">idx</span>.<span class="pl-en">reshape</span>(<span class="pl-c1">-</span><span class="pl-c1">1</span>, <span class="pl-c1">1</span>), <span class="pl-s1">idx</span>]
<span class="pl-c"># peak memory: 135.93 MiB, increment: 46.52 MiB</span>

<span class="pl-c"># This does not</span>
<span class="pl-c1">%</span><span class="pl-s1">memit</span> <span class="pl-s1">s</span>[<span class="pl-s1">np</span>.<span class="pl-en">ix_</span>(<span class="pl-s1">idx</span>, <span class="pl-s1">idx</span>)]
<span class="pl-c"># peak memory: 2066.09 MiB, increment: 1930.15 MiB</span></pre></div>
<p dir="auto">As far as I can tell, the results are exactly the same. Only performance is different.</p>
<p dir="auto">It looks like this is (at least partially) a known issue (<a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="611576708" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/12015" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/12015/hovercard?comment_id=625854928&amp;comment_type=issue_comment" href="https://github.com/scipy/scipy/issues/12015#issuecomment-625854928">#12015 (comment)</a>).</p>
<p dir="auto">Initially discovered in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="621176807" data-permission-text="Title is private" data-url="https://github.com/scverse/anndata/issues/380" data-hovercard-type="issue" data-hovercard-url="/scverse/anndata/issues/380/hovercard" href="https://github.com/scverse/anndata/issues/380">scverse/anndata#380</a></p>
<h4 dir="auto">Scipy/Numpy/Python version information:</h4>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="# 1.4.1 1.18.4 sys.version_info(major=3, minor=7, micro=7, releaselevel='final', serial=0)"><pre class="notranslate"><code class="notranslate"># 1.4.1 1.18.4 sys.version_info(major=3, minor=7, micro=7, releaselevel='final', serial=0)
</code></pre></div>
<p dir="auto">A typical way to perform outer indexing can be unnecessarily and surprisingly super slow. It is illustrated in the in/out below.</p>
<p dir="auto">I think the issue could be fixed by adding another fast path for input as in my <code class="notranslate">In [8]</code> to <code class="notranslate">scipy._index.IndexMixin</code>, in a similar manner to how it's done when using the method of <code class="notranslate">In [10]</code>.</p>
<h4 dir="auto">Reproducing code example:</h4>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="In [1]: import scipy.sparse as ss, numpy as np                                                                                       

In [2]: matrix = ss.random(2000, 50000, format=&quot;csr&quot;)                                                                                

In [3]: idx0 = np.arange(1000)[:, None]                                                                                              

In [4]: idx1 = np.arange(1000, 2000)                                                                                                 

In [5]: imod = np.arange(25)[None, :]                                                                                                

In [6]: idx1 = (idx1[:, None] + 1000 * imod).T.reshape(-1, 1)                                                                        

In [8]: # The 'easy' way to index: also really slow 
   ...: %timeit out0 = matrix[idx0, idx1.T]                                                                                          
743 ms � 48 ms per loop (mean � std. dev. of 7 runs, 1 loop each)

In [9]: # Equivalent submatrix, but not so easy to write: Much faster 
   ...: %timeit out1 = matrix[idx0.ravel(), :][:, idx1.ravel()]                                                                      
3.57 ms � 89.6 �s per loop (mean � std. dev. of 7 runs, 100 loops each)

In [10]: # There is an existing fast path when the second index is 1D 
    ...: %timeit out2 = matrix[idx0, idx1.ravel()]                                                                                   
3.47 ms � 113 �s per loop (mean � std. dev. of 7 runs, 100 loops each)"><pre class="notranslate"><span class="pl-v">In</span> [<span class="pl-c1">1</span>]: <span class="pl-s1">import</span> <span class="pl-s1">scipy</span>.<span class="pl-s1">sparse</span> <span class="pl-k">as</span> <span class="pl-s1">ss</span>, <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>                                                                                       

<span class="pl-v">In</span> [<span class="pl-c1">2</span>]: <span class="pl-s1">matrix</span> <span class="pl-c1">=</span> <span class="pl-s1">ss</span>.<span class="pl-en">random</span>(<span class="pl-c1">2000</span>, <span class="pl-c1">50000</span>, <span class="pl-s1">format</span><span class="pl-c1">=</span><span class="pl-s">"csr"</span>)                                                                                

<span class="pl-v">In</span> [<span class="pl-c1">3</span>]: <span class="pl-s1">idx0</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">1000</span>)[:, <span class="pl-c1">None</span>]                                                                                              

<span class="pl-v">In</span> [<span class="pl-c1">4</span>]: <span class="pl-s1">idx1</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">1000</span>, <span class="pl-c1">2000</span>)                                                                                                 

<span class="pl-v">In</span> [<span class="pl-c1">5</span>]: <span class="pl-s1">imod</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">25</span>)[<span class="pl-c1">None</span>, :]                                                                                                

<span class="pl-v">In</span> [<span class="pl-c1">6</span>]: <span class="pl-s1">idx1</span> <span class="pl-c1">=</span> (<span class="pl-s1">idx1</span>[:, <span class="pl-c1">None</span>] <span class="pl-c1">+</span> <span class="pl-c1">1000</span> <span class="pl-c1">*</span> <span class="pl-s1">imod</span>).<span class="pl-v">T</span>.<span class="pl-en">reshape</span>(<span class="pl-c1">-</span><span class="pl-c1">1</span>, <span class="pl-c1">1</span>)                                                                        

<span class="pl-v">In</span> [<span class="pl-c1">8</span>]: <span class="pl-c"># The 'easy' way to index: also really slow </span>
   ...: <span class="pl-c1">%</span><span class="pl-s1">timeit</span> <span class="pl-s1">out0</span> <span class="pl-c1">=</span> <span class="pl-s1">matrix</span>[<span class="pl-s1">idx0</span>, <span class="pl-s1">idx1</span>.<span class="pl-v">T</span>]                                                                                          
<span class="pl-c1">743</span> <span class="pl-s1">ms</span> � <span class="pl-c1">48</span> <span class="pl-s1">ms</span> <span class="pl-s1">per</span> <span class="pl-en">loop</span> (<span class="pl-s1">mean</span> � <span class="pl-s1">std</span>. <span class="pl-s1">dev</span>. <span class="pl-s1">of</span> <span class="pl-c1">7</span> <span class="pl-s1">runs</span>, <span class="pl-c1">1</span> <span class="pl-s1">loop</span> <span class="pl-s1">each</span>)

<span class="pl-v">In</span> [<span class="pl-c1">9</span>]: <span class="pl-c"># Equivalent submatrix, but not so easy to write: Much faster </span>
   ...: <span class="pl-c1">%</span><span class="pl-s1">timeit</span> <span class="pl-s1">out1</span> <span class="pl-c1">=</span> <span class="pl-s1">matrix</span>[<span class="pl-s1">idx0</span>.<span class="pl-en">ravel</span>(), :][:, <span class="pl-s1">idx1</span>.<span class="pl-en">ravel</span>()]                                                                      
<span class="pl-c1">3.57</span> <span class="pl-s1">ms</span> � <span class="pl-c1">89.6</span> �<span class="pl-s1">s</span> <span class="pl-s1">per</span> <span class="pl-en">loop</span> (<span class="pl-s1">mean</span> � <span class="pl-s1">std</span>. <span class="pl-s1">dev</span>. <span class="pl-s1">of</span> <span class="pl-c1">7</span> <span class="pl-s1">runs</span>, <span class="pl-c1">100</span> <span class="pl-s1">loops</span> <span class="pl-s1">each</span>)

<span class="pl-v">In</span> [<span class="pl-c1">10</span>]: <span class="pl-c"># There is an existing fast path when the second index is 1D </span>
    ...: <span class="pl-c1">%</span><span class="pl-s1">timeit</span> <span class="pl-s1">out2</span> <span class="pl-c1">=</span> <span class="pl-s1">matrix</span>[<span class="pl-s1">idx0</span>, <span class="pl-s1">idx1</span>.<span class="pl-en">ravel</span>()]                                                                                   
<span class="pl-c1">3.47</span> <span class="pl-s1">ms</span> � <span class="pl-c1">113</span> �<span class="pl-s1">s</span> <span class="pl-s1">per</span> <span class="pl-en">loop</span> (<span class="pl-s1">mean</span> � <span class="pl-s1">std</span>. <span class="pl-s1">dev</span>. <span class="pl-s1">of</span> <span class="pl-c1">7</span> <span class="pl-s1">runs</span>, <span class="pl-c1">100</span> <span class="pl-s1">loops</span> <span class="pl-s1">each</span>)</pre></div>
<h4 dir="auto">Scipy/Numpy/Python version information:</h4>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="In [40]: import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)                                     
1.3.0 1.16.4 sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)"><pre class="notranslate"><span class="pl-v">In</span> [<span class="pl-c1">40</span>]: <span class="pl-s1">import</span> <span class="pl-s1">sys</span>, <span class="pl-s1">scipy</span>, <span class="pl-s1">numpy</span>; <span class="pl-en">print</span>(<span class="pl-s1">scipy</span>.<span class="pl-s1">__version__</span>, <span class="pl-s1">numpy</span>.<span class="pl-s1">__version__</span>, <span class="pl-s1">sys</span>.<span class="pl-s1">version_info</span>)                                     
<span class="pl-c1">1.3</span>.<span class="pl-c1">0</span> <span class="pl-c1">1.16</span>.<span class="pl-c1">4</span> <span class="pl-s1">sys</span>.<span class="pl-en">version_info</span>(<span class="pl-s1">major</span><span class="pl-c1">=</span><span class="pl-c1">3</span>, <span class="pl-s1">minor</span><span class="pl-c1">=</span><span class="pl-c1">7</span>, <span class="pl-s1">micro</span><span class="pl-c1">=</span><span class="pl-c1">3</span>, <span class="pl-s1">releaselevel</span><span class="pl-c1">=</span><span class="pl-s">'final'</span>, <span class="pl-s1">serial</span><span class="pl-c1">=</span><span class="pl-c1">0</span>)</pre></div>