Minimize memory overhead in sparse matrix inverse
Question: As pretense, I am continuing development in Python 2.7 from a prior question:
[Determining a sparse matrix
quotient](http://stackoverflow.com/questions/40050947/determining-a-sparse-
matrix-quotient)

## My existing 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)
    

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

## Error Statement:

    
    
    Not enough memory to perform factorization.
    Traceback (most recent call last):
      File "<myfilename.py>", line 111, in <module>
        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
    

From the [scipy/smemory.c
sourcecode](https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/dsolve/SuperLU/SRC/smemory.c),
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.

For reference,

`numerator` has `shape: (552297, 552297)` with `stored elements: 301067607`
calculated as `sp.csr_matrix(A.T.dot(Ap))`

`denominator` has `shape: (552297, 552297)` with `stored elements: 170837213`
calculated as `sp.csr_matrix(A.T.dot(A))`

**EDIT** : I've found [a related question on
Reddit](https://www.reddit.com/r/Python/comments/3c0m7b/what_is_the_most_precise_way_to_invert_large/),
but cannot determine how I would change my equation from `numerator *
inv(denominator) = k`


Answer: No need to 'preallocate' `k`; this isn't a compiled language. Not that this is
costing anything.

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

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

    
    
    denom = sp.csc_matrix(denominator)
    #halfeq = sp.linalg.inv(denom)
    #k = numerator.dot(halfeq)
    k = sp.linalg.spsolve(denom, numerator)
    

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

Maybe it would help to step back and look at:

    
    
    num = A.T.dot(Ap)
    den = A.T.dot(A)
    k = solve(den, num)
    

In other words, review the matrix algebra.

    
    
    (A'*Ap)/(A'*A)
    

I'm little rusty on this. Can we reduce this? Can we partition?

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

How about providing small `A` and `Ap` 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.

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

