How to tell if a sparse matrix can't be solved

Question

I'm writing a program in Python using scipy's spsolve to solve a linear equation using a sparse matrix (csr_matrix). The matrices are fairly large (M=90826x90826, b=90826x1) and are hard to check by hand.

The problem I've encountered is that for a small fraction of the matrices I'm creating, scipy.sparse.linalg.spsolve(M,b) undergoes such a catastrophic failure that it crashes the entire program. Even wrapping the line in a try/except and the program itself in another try/except, doesn't help. I don't even get a response in the Exception field.

I verified manually that the data I'm sending it is incorrect for what I want, but I can't really check it anywhere else. My program involves face detection, and the problem may be caused by bad detection. In this case, the "face" was found to be in the cheek of the actual face. However, manually verifying that the detection is correct before proceeding is not really an option (the final product will be used by a non-technical client). And automatic verification is outside the scope of the project.

Anyway, if I can just detect that the matrices will cause spsolve to crash, it's fine to simply skip the image. But I can't seem to find literature about how to keep spsolve from crashing.

The values I'm sending of type float64 and can be positive or negative. From what I've seen of the matrices that can be solved, M is generally filled with 4's and -1's, while b can be almost any positive or negative number.

Answer 1

Without knowing the exact error, it's hard to say what's going wrong. I'm not overly familiar with scipy, but I suspect if there was no solution to these problems due to an inconsistent system, you would get a meaningful error.

My best guess would be a memory issue. During Gaussian elimination, sparse matrices can undergo a large amount of fill-in, where zeros become non-zeros. For example, one non-zero along the top row could potentially result in all the zeros below it up to the diagonal being filled in (set to non-zero). Have you experimented with different node reorderings (permc_spec parameter)? The default should be a pretty good job, but given there are a few in-built options I think it would be silly not to try them out.

There's a very good description of how this works (with pictures!) here (though this is mathworks site, so any implementation will be different to scipy).

Alternatively, if you can accept a close-but-not-exact answer, there are plenty of iterative methods that can get an approximate answer in a fraction of the time and memory requirements. Without knowing more about the nature of the matrix, it's hard to say which would be best - but you can experiment with any of the functions under 'Iterative methods for linear equation systems' listed here.