Minimize memory overhead in sparse matrix inverse
<p>As pretense, I am continuing development in Python 2.7 from a prior question: <a href="http://stackoverflow.com/questions/40050947/determining-a-sparse-matrix-quotient">Determining a sparse matrix quotient</a> </p>

<h2>My existing code:</h2>

<pre><code>import scipy.sparse as sp

k = sp.csr_matrix(([], ([],[])),shape=[R,R])

denom = sp.csc_matrix(denominator)
halfeq = sp.linalg.inv(denom)
k = numerator.dot(halfeq)
</code></pre>

<p>I was successful in calculating for the base <code>k</code> and <code>denom</code>. Python continued attempting calculation on <code>halfeq</code>. The process sat in limbo for aproximately 2 hours before returning an error </p>

<h2>Error Statement:</h2>

<pre><code>Not enough memory to perform factorization.
Traceback (most recent call last):
  File "&lt;myfilename.py&gt;", line 111, in &lt;module&gt;
    halfeq = sp.linalg.inv(denom)
  File "/opt/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.py", line 61, in inv Ainv = spsolve(A, I)
  File "/opt/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 151, in spsolve Afactsolve = factorized(A)
  File "/opt/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 366, in factorized return splu(A).solve
  File "/opt/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 242, in splu ilu=False, options=_options)
MemoryError
</code></pre>

<p>From the <a href="https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/dsolve/SuperLU/SRC/smemory.c" rel="nofollow">scipy/smemory.c sourcecode</a>, the initial statement from the error is found on line 256. I am unable to further analyze the memory defs to determine how to best reallocate memory usage sufficient for execution. </p>

<p>For reference, </p>

<p><code>numerator</code> has <code>shape: (552297, 552297)</code> with <code>stored elements: 301067607</code> calculated as <code>sp.csr_matrix(A.T.dot(Ap))</code></p>

<p><code>denominator</code> has <code>shape: (552297, 552297)</code> with <code>stored elements: 170837213</code> calculated as <code>sp.csr_matrix(A.T.dot(A))</code></p>

<p><strong>EDIT</strong>: I've found <a href="https://www.reddit.com/r/Python/comments/3c0m7b/what_is_the_most_precise_way_to_invert_large/" rel="nofollow">a related question on Reddit</a>, but cannot determine how I would change my equation from <code>numerator * inv(denominator) = k</code></p>

<p>No need to 'preallocate' <code>k</code>; this isn't a compiled language.  Not that this is costing anything.</p>

<pre><code>k = sp.csr_matrix(([], ([],[])),shape=[R,R])
</code></pre>

<p>I need to double check this, but I think the <code>dot/inv</code> can be replaced by one call to <code>spsolve</code>.  Remember in the other question I noted that <code>inv</code> is <code>spsolve(A, I)</code>; </p>

<pre><code>denom = sp.csc_matrix(denominator)
#halfeq = sp.linalg.inv(denom)
#k = numerator.dot(halfeq)
k = sp.linalg.spsolve(denom, numerator)
</code></pre>

<p>That said, it looks like the problem is in the <code>inv</code> part, the <code>factorized(denom)</code>.  While your arrays are sparse, (denom density is 0.00056), they still have a large number of values.</p>

<p>Maybe it would help to step back and look at:</p>

<pre><code>num = A.T.dot(Ap)
den = A.T.dot(A)
k = solve(den, num)
</code></pre>

<p>In other words, review the matrix algebra.</p>

<pre><code>(A'*Ap)/(A'*A)
</code></pre>

<p>I'm little rusty on this.  Can we reduce this?  Can we partition?  </p>

<p>Just throwing great big arrays together, even if they are sparse, isn't working.</p>

<p>How about providing small <code>A</code> and <code>Ap</code> arrays that we can use for testing?  I'm not interested in testing memory limits, but I'd like to experiment with different calculation methods.</p>

<p>The sparse linalg module has a number of iterative solvers.  I have no idea whether their memory use is greater or less.</p>

python scipy sparse-matrix