LU errors out for sparse matrices over inexact fields

[maxale]

Steps To Reproduce

In the following code A.LU() works fine, but B.LU() errors out.

A = random_matrix(RDF,4)
print(A.LU())
B = A.sparse_matrix()
print(B.LU())

Expected Behavior

.LU() should work for both dense and sparse matrices.

Actual Behavior

The following error is produced:

TypeError: base ring of the matrix must be exact, not Real Double Field

which is quite weird given that this was not an issue for matrix A also defined over RDF.

Additional Information

No response

Environment

- **OS**: Ubuntu 22.04.3 LTS
- **Sage Version**: 10.2.beta5

Checklist

• I have searched the existing issues for a bug report that matches the one I want to file, without success.
• I have read the documentation and troubleshoot guide

[Shay2Shay]

I went through the code
inside : sage/src/sage/matrix

matrix2.pyx defines LU() which checks whether self.base_ring is exact or not : https://github.com/sagemath/sage/blob/439065ecf90026fb35d22edd44ef76d8d2ed9fd4/src/sage/matrix/matrix2.pyx#L13553C9-L13553C12

For dense matrix : inherits matrix2.pyx > matrix_double_dense.pyx
matrix_double_dense.pyx redefines (overloads) LU() function and uses scipy.linalg.lu :

sage/src/sage/matrix/matrix_double_dense.pyx

Line 1043 in 439065e

PM, LM, UM = scipy.linalg.lu(self._matrix_numpy)

But for sparse matrix : they inherit > matrix.pyx > matrix_sparse.pyx > matrix_generic_sparse.pyx > matrix_double_sparse.pyx
Neither of them redefine LU factorization for sparse matrix

matrix.pyx == matrix2.pyx

sage/src/sage/matrix/matrix.pxd

Line 1 in 439065e

from sage.matrix.matrix2 cimport Matrix

---

We need to implement LU() in matrix_sparse.pyx or matrix_generic_sparse.pyx or matrix_double_sparse.pyx
Scipy may help : https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html#scipy-sparse-linalg-splu

I would like to take it.

[Shay2Shay]

@tscrim @fchapoton @mkoeppe
I am sorry I did not know whom to tag.

[tscrim]

I am not sure this is a bug; although it certainly is surprising (at first) behavior. The issue is what does zero mean for an inexact ring. Is it something we know is zero or something that is not distinguishable from zero? This could cause subtle issues for sparse matrices that are not present for dense matrices.

It might be better to raise a NotImplementedError, but I am not convinced of that as it might not be something we can implement (with a reasonable specification).

For the RDF case, it might be good for the specific sparse implementation to just convert itself to a dense matrix and return the LU decomposition of that. Although one must be careful with this as making sparse matrices dense could use all the system's memory. However, this is done in some places, but we really should not do that IMO.

I am not sure what the best way forward is. It sometimes is important for users to know what type of object they are working with as those implementation details do matter at times. (In other words, we need to improve the documentation.)

[Shay2Shay]

Here is what I found on internet, please verify my understanding:

1. An exact ring is one whose elements are represented exactly in the system without approximation.
2. In an inexact ring, the value of zero is not always distinguishable from other values that are very close to zero.
   And an Example:
   Consider a sparse matrix with entries that are very small but not exactly zero. In a dense matrix, these entries would be represented as very small values, but in a sparse matrix, they would be represented as non-zero entries. This can lead to unexpected results when performing operations on sparse matrices.

I had a doubt: what does non-zero entries mean (in context of sage)?

Edited:

sage: RDF(0) == RDF(1e-323)
False
sage: RDF(0) == RDF(1e-324)
True
sage:

I tried different rings:

1. Exact ring: ZZ, QQ
2. Inexact ring: RDF, CDF

I believe exactness arises due to implementation aspect of Field elements.

[Shay2Shay]

As you mentioned there are 2 issues with sparse matrix LU decomposition implementation:

1. Meaning of zero, distinguishability and exactness
2. Memory overflow for larger sparse matrices
   And we also want for user to be more involved in this operation (LU decomposition of sparse matrix) we can go with your idea of Not Implemented error with message "convert to dense 1st then perform LU". We can override the LU in matrix_sparse.pyx or matrix_generic_sparse.pyx with doc string giving example.

@tscrim

[maxale]

The issue is what does zero mean for an inexact ring.

Over inexact fields, this issue already exists for conversion from dense matrices to sparse ones. Why then B = A.sparse_matrix() is somehow implemented, but B.LU() is not?

[maxale]

Also, it's worth mentioning that MATLAB does support LU factorization for sparse matrices (see this page for technical details). In the freeware world, such functionality is supported by GNU Octave via UMFPACK.

[Shay2Shay]

Over inexact fields, this issue already exists for conversion from dense matrices to sparse ones. Why then B = A.sparse_matrix() is somehow implemented, but B.LU() is not?

What problem does sparse inexact ring matrices raise that dense ones don't?

[maxale]

What problem does sparse inexact ring matrices raise that dense ones don't?

Haven't you already described the problems just a few messages above? Or please elaborate your question.

[Shay2Shay]

Docs of scipy.sparse.linalg.splu says:

1. Data will be represented as CSC format, if not then the input will be converted into CSC format. https://github.com/scipy/scipy/blob/5e4a5e3785f79dd4e8930eed883da89958860db2/scipy/sparse/linalg/_dsolve/linsolve.py#L343
2. It can only factor square matrices. https://github.com/scipy/scipy/blob/5e4a5e3785f79dd4e8930eed883da89958860db2/scipy/sparse/linalg/_dsolve/linsolve.py#L420C10-L420C10
3. splu uses superLU and at one point it allocates memory proportional to num columns 2 times https://github.com/scipy/scipy/blob/8e9d43838b7db8d40fc4567f3aafeed68466b464/scipy/sparse/linalg/_dsolve/_superluobject.c#L723

This was basic memory requirement review for scipy sparse LU decomposition and also square matrix limitation.