Eigen: Demonstrate and fix issue with passing scipy.sparse matrices with unsorted indices

[ManifoldFR]

Hi,

I've just encountered an issue when the user passes a scipy sparse matrix with unsorted inner indices to Eigen.
This issue manifests in two ways:

• an eigen_assert() is triggered when compiling in Debug,
• or you can get incorrect values in your converted sparse matrix.

It seems until relatively recently it was not explicitly stated that Eigen::Map<SparseMatrix> objects expect the inner indices to be sorted. This issue really pops up when you try to evaluate the map into a plain sparse matrix object, like here:

nanobind/include/nanobind/eigen/sparse.h

Line 89 in 10eafbc

value = SparseMap(rows, cols, nnz, outer_indices.data(), inner_indices.data(), values.data());

Eigen's master branch has a nice a.sortInnerIndices() function and some helpers to check if the indices are sorted, for compressed storage types. For Eigen's released version (3.4.0 and below), for now I've added a call to SciPy's API to sort the indices of the passed sparse matrix.

This might be relevant for #782.
Also related: stack-of-tasks/eigenpy#538 (in Boost.Python world)

[wjakob]

Actually, I am wondering about something here: isn't it unsafe to use mapped_matrix.sortInnerIndices(); because it will modify the contents of the mapped index/value arrays? (which are owned by SciPy/NumPy). In that case, the safer alternative would seem to be to use the SciPy function.

The advantage of this is also that we can avoid sorting if the flag is True. It seems to me that the Eigen sortInnerIndices function in contrast doesn't know if the indices are sorted and will have to conservatively invoke a sort routine at some cost.

[ManifoldFR]

Actually, I am wondering about something here: isn't it unsafe to use mapped_matrix.sortInnerIndices(); because it will modify the contents of the mapped index/value arrays? (which are owned by SciPy/NumPy). In that case, the safer alternative would seem to be to use the SciPy function.

You might be right, we would be modifying SciPy's indices array from mapped data in C++, which SciPy itself might not expect, whereas calling the SciPy function would be safer (and would call whatever internal code the SciPy function has).

The advantage of this is also that we can avoid sorting if the flag is True. It seems to me that the Eigen sortInnerIndices function in contrast doesn't know if the indices are sorted and will have to conservatively invoke a sort routine at some cost.

I hadn't looked at the code that much in-depth. I think we can go with your suggestion.

[wjakob]

Thanks!