<p dir="auto">I am trying to solve a linear system as effeciently as possible, and one of the methods I have been trying is by reordering the matrix beforehand with reverse_cuthill_mckee.<br>
This works well for small matrices, but when I move to a large size this breaks. (60000x60000)</p>
<p dir="auto">The code to reproduce this is really simple:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="perm = scipy.sparse.csgraph.reverse_cuthill_mckee(H, symmetric_mode=True)
I,J = np.ix_(perm,perm)
H_reordered = H[I,J]"><pre class="notranslate"><code class="notranslate">perm = scipy.sparse.csgraph.reverse_cuthill_mckee(H, symmetric_mode=True)
I,J = np.ix_(perm,perm)
H_reordered = H[I,J]
</code></pre></div>
<p dir="auto">Where H is a symmetric 60000x60000 sparse csc matrix with some nonzero elements. In my case it has about 7 million nonzero elements</p>
<p dir="auto">When I run the above I get:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="Traceback (most recent call last):
  File &quot;C:\Program Files\JetBrains\PyCharm Community Edition 2020.1\plugins\python-ce\helpers\pydev\pydevd.py&quot;, line 1438, in _exec
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File &quot;C:\Program Files\JetBrains\PyCharm Community Edition 2020.1\plugins\python-ce\helpers\pydev\_pydev_imps\_pydev_execfile.py&quot;, line 18, in execfile
    exec(compile(contents+&quot;\n&quot;, file, 'exec'), glob, loc)
  File &quot;D:/Dropbox/ActiveLearning/Pytorch/Active-Learning/pcg.py&quot;, line 37, in &lt;module&gt;
    H_ref = H[I,J]
  File &quot;D:\Dropbox\ActiveLearning\Pytorch\Active-Learning\venv\lib\site-packages\scipy\sparse\_index.py&quot;, line 75, in __getitem__
    return self._get_arrayXarray(row, col)
  File &quot;D:\Dropbox\ActiveLearning\Pytorch\Active-Learning\venv\lib\site-packages\scipy\sparse\compressed.py&quot;, line 664, in _get_arrayXarray
    csr_sample_values(M, N, self.indptr, self.indices, self.data,
ValueError: could not convert integer scalar"><pre class="notranslate"><code class="notranslate">Traceback (most recent call last):
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2020.1\plugins\python-ce\helpers\pydev\pydevd.py", line 1438, in _exec
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2020.1\plugins\python-ce\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "D:/Dropbox/ActiveLearning/Pytorch/Active-Learning/pcg.py", line 37, in &lt;module&gt;
    H_ref = H[I,J]
  File "D:\Dropbox\ActiveLearning\Pytorch\Active-Learning\venv\lib\site-packages\scipy\sparse\_index.py", line 75, in __getitem__
    return self._get_arrayXarray(row, col)
  File "D:\Dropbox\ActiveLearning\Pytorch\Active-Learning\venv\lib\site-packages\scipy\sparse\compressed.py", line 664, in _get_arrayXarray
    csr_sample_values(M, N, self.indptr, self.indices, self.data,
ValueError: could not convert integer scalar
</code></pre></div>
<p dir="auto">And the reason as far as I can see is that the reordering involves making a 60000x60000 ndarray, which has about 3.6*10^9 elements and hence is too big int32.</p>
<p dir="auto">Does anyone know of a workaround for this problem? or a fix?</p>
<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>