scipy: BUG: Memory Error in scipy.sparse.splu

Describe your issue.

Hi. I have been using scipy.sparse.linalg.splu for sparse LU factorization. I kept running into memory issues with problems that were nowhere near the RAM capacity of the system. After more experimenting, I am quite sure that the error has nothing to do with memory at all. I was able to reproduce the error on 2 different systems, one 128GPs of RAM and the other one 32 GB of RAM. The reproducing code example reliable leads to the error on the 32GP system, going nowhere near the available RAM. Thanks for any advice in advance.

Reproducing Code Example

from numpy.random import default_rng
import numpy as np
from scipy.sparse.linalg import splu
import scipy.sparse as sparse
import time
from scipy.sparse.linalg import use_solver

def add_to_diag(res, vector):
    diag = sparse.eye(len(vector), format="coo") ##make variance
    diag.setdiag(vector) ##make variance
    sparse_covariance = res + diag  ##add variance
    return sparse_covariance

def maked(n,m):
    rng1 = default_rng()
    rows = rng1.choice(n, size=m)
    rng2 = default_rng()
    cols = rng2.choice(n, size=m)
    data = np.random.rand(m)
    return rows,cols,data



use_solver(useUmfpack=False)

N = 10000
n = 2000000

matrix = sparse.coo_matrix((N,N))

r,c,d = maked(N,n)

res = sparse.coo_matrix((d,(r,c)), dtype = np.float64)
res = res @ res.T

vector = np.random.rand(N) * 1.0
res = add_to_diag(res,vector)

print("Number of Non-Zero Entries: ",res.count_nonzero())

st = time.time()
#########LU###########
A = res.tocsc()
print("Transferred after: ", time.time() - st)
LU = splu(A)

print(LU)
print("concluded after: ", time.time() - st)

Error message

---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Input In [2], in <module>
     43 A = res.tocsc()
     44 print("Transferred after: ", time.time() - st)
---> 45 LU = splu(A)
     47 print(LU)
     48 print("concluded after: ", time.time() - st)

File ~/VirtualEnvironments/fvgp_env/lib/python3.8/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:337, in splu(A, permc_spec, diag_pivot_thresh, relax, panel_size, options)
    334 if (_options["ColPerm"] == "NATURAL"):
    335     _options["SymmetricMode"] = True
--> 337 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
    338                       csc_construct_func=csc_construct_func,
    339                       ilu=False, options=_options)

MemoryError:

SciPy/NumPy/Python version information

1.7.3 1.22.1 sys.version_info(major=3, minor=8, micro=10, releaselevel=‘final’, serial=0)

About this issue

Original URL
State: open
Created 2 years ago
Comments: 34 (9 by maintainers)

Most upvoted comments

Hi @GregStroot. This is a SuperLU “issue”, meaning that it is partly expected behavior because of the structure of the matrices. I am working with the developer of SuperLU on trying different sorting mechanisms. I’ll post any results here.

MarcusMNoack on May 5, 2022

Maybe even some of the scipy patch can be upstreamed to SuperLU?

ev-br on Apr 14, 2023

So will that happen automatically?

No. Someone will have to do the work to replace the version that is vendored in SciPy (the README file says 5.2.x) with the latest release. We currently have a patch file (scipychanges.patch) containing changes that SciPy made to the vendored copy, so the update is not a trivial task.

I think we want to update, but someone will have to do the work to make it happen, and I don’t know if any of the current devs has this on their “to do” list at the moment. Anyone who takes on this task should have some familiarity with C, and be willing to dig into the history to understand the changes in the patch file.

WarrenWeckesser on Apr 13, 2023

The SuperLU developer fixed the issue. This can now be included in scipy. Please let me know if this will be done and, if so, when it’s done. Thank you all for figuring this out. The 6.0.0 release is here: https://github.com/xiaoyeli/superlu/releases

MarcusMNoack on Apr 7, 2023

@nicolasdonati The solution I found was to switch solvers (e.g. not use scipy). I am now using petsc4py, which allows you to access the MUMPS LU solver. I have successfully ran LU decompositions on matricies (N x N) of size N = 9e6 with complex double data-types, which took ~600GB of memory

Here is some code that will allow you to use it

        import scipy.sparse as spp
        import petsc4py 
        from petsc4py import PETSC

        def import_petsc(mat):
            mat = spp.csr_matrix(mat) #IMPORTANT. Must be in csr format
                                      #csc results in being imported transposed
            mat_petsc = PETSc.Mat().createAIJ(size=mat.shape,
                                              csr = (mat.indptr, mat.indices,
                                                     mat.data))
            mat_petsc.setUp()
            mat_petsc.assemble()
            return mat_petsc

        LHS_petsc = import_petsc(LHS)
        ksp = PETSc.KSP()
        ksp.create(PETSc.COMM_WORLD)
        ksp.setOperators(LHS_petsc,None)
        pc_ksp = ksp.getPC()
        pc_ksp.setType('lu')
        pc_ksp.setFactorSolverType('mumps') #Mumps runs in parallel and complex
        ksp.setFromOptions()

        v_in_rhs = np.random.rand(N) 
            #Create array to store answer and Convert to petsc form
            v_out = PETSc.Vec().createSeq(len(v_in_rhs))
            v_in_rhs = PETSc.Vec().createWithArray(v_in_rhs)
            ksp.solve(v_in_rhs, v_out)
            #Back to numpy
            v_out = v_out.getArray()

Feel free to send a message if you get stuck. The documentation is a bit unwieldy

GregStroot on Mar 7, 2023