ENH: sparse.linalg: lobpcg and iterative linear solvers: add array API support

[lobpcg]

Is your feature request related to a problem? Please describe.

Sparse lobpcg and linear iterative solvers are pure Python and ideal candidates for array API support.

Describe the solution you'd like.

No response

Describe alternatives you've considered.

No response

Additional context (e.g. screenshots, GIFs)

Is array API support in scipy already mature enough to be added to these solvers?

Maybe the effort for lobpcg can be coordinated with https://github.com/cupy/cupy/issues/7925 ?

[lucascolley]

See gh-18915 for discussion and gh-20190 for a first attempt towards an array API compatible namespace for scipy.sparse. I imagine that cupyx.scipy.sparse will follow in the footsteps after this work is done.

Adding array API standard support to these functions would make them interoperable with dense array libraries at the moment, but I don't know if that is all that valuable right now.

[lucascolley]

That said, I would be happy to review a PR for this if you would like to submit one!

[lobpcg]

Thanks for your comments! gh-20190 concerns specifically and only sparse matrix formats. LOBPCG and iterative linear solvers in sparse are all matrix-free instead require only a matvev function.

Linear solvers rely on LinearOperator, while LOBPCG goes further and allows calling matmat function directly without the LinearOperator structure so may be the easiest.

Not all dense np arrays in LOBPCG require changes at the first stage, only those of large sizes. Possibly helping could be the fact that there are already separate implementations available of LOBPCG, e.g., in cupy. The effort could be to come up with a single flexible implementation based on array API.

[lucascolley]

The effort could be to come up with a single flexible implementation based on array API.

Right, and what I'm saying is that the implementation based on the array API will only work for data structures with a compatible interface. So we could get cupy.ndarray support, but not cupy.scipy.sparse.csr_matrix. Of course, that could be added in addition to array API support by calling out to cupy, but that doesn't in itself help reduce the duplication of the implementations.

In other words, we could have the solvers in SciPy which support the array API standard, but we would need to keep implementations for sparse types separate until they can be captured via the same API, hence the link to gh-20190. So I wonder whether it is worth doing this now or waiting until the sparse libraries gel better with the array API.

[lobpcg]

LOBPCG and iterative linear solvers are in sparse only for historical reasons. In reality, they are all matrix-free require only a user-provided matvec function. So their array API implementation is independent of that for sparse arrays and actually rely only on array API for dense 1D and, for LOBPCG, 2D arrays.

[lucascolley]

scipy.sparse.linalg.lobpcg is documented to accept a sparse matrix for A. Unless you would propose to deprecate that, we would have to maintain support. In fact, the docs say "usually given by a sparse matrix".

[rkern]

No one is suggesting removing support for sparse matrix objects.

[lucascolley]

Okay, if support can be added whilst maintaining support for sparse data structures, without duplicating implementations, then feel free to submit a PR!

[lobpcg]

scipy.sparse.linalg.lobpcg is documented to accept a sparse matrix for A. Unless you would propose to deprecate that, we would have to maintain support. In fact, the docs say "usually given by a sparse matrix".

Yes, but the code only needs it to perform A @ X where X is a dense array, and only other X-like dense arrays are being modified in these matrix-free solvers and require array API rather than A that can remain in the traditional format.