`solve_left` for sparse matrices extremely slow

[Peter Mueller]

Set:

n = 100000

A = matrix(QQ, n, 1, range(n))
B = A.sparse_matrix()
v = vector((0,))

The  trivial calculation

sol = A.solve_left(v)

instantly yields a solution (the all 0 vector), while

sol = B.solve_left(v)

takes forever ... !!! I haven't looked at the implementation, but it seems to me that there is a highly inefficient handling of sparse matrices.

-- Peter Mueller

[Thierry Dumont]

Sparse matrices are used by people who --like me-- solve PDEs
discretized by, say finite elements or finite volumes on irregular
meshes. In these case the matrix is sparse, that is to say it has O(n)
non zero terms, for a matrix of size n. In real computations, n can be
much more than 10^6 (10^8 or even more). So, you need something like a map:
(i,j)-> a_{ij}
to represent the matrix and store only the non zero terms. Have a look
at "csr" matrices for example.

These data structure cannot be avoided, but they result in slow
computations (because of indirection and low arithmetic intensity,
whatever you do). As a consequence, when solving PDEs, the less you
solve linear systems, the best it is !


If your matrix is banded you can compute faster (even if you need to
store some 0) using the banded matrices in lapack.

Consider also that:

- for full or banded matrix lapack and blas are extremely well optimized,

- but it seems that, for sparse matrices, Sage uses the scipy
implementation which is written in python, and... this is really not
good for speed ! So what you say is true : a slow implementation and
slow data structures...


Yours,
Thierry

> --
> You received this message because you are subscribed to the Google
> Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sage-support...@googlegroups.com
> <mailto:sage-support...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sage-support/a051a430-276b-445b-ae36-6afba9e6aecfn%40googlegroups.com <https://groups.google.com/d/msgid/sage-support/a051a430-276b-445b-ae36-6afba9e6aecfn%40googlegroups.com?utm_medium=email&utm_source=footer>.

[julian...@gmail.com]

Hi Thierry,

The sparse matrix actually uses a multimodular algorithm in this case. It's written in Cython but essentially runs at Python speed. However, that's not where the slowness comes from. Even when running the dense matrix with that same algorithm, it's still very fast.

I agree that fundamentally sparse matrices cannot be as fast. But in my experience they are often also not very optimized in SageMath. But here it seems that something is actually broken because what takes so long is computing the transpose of the input matrix. (To reduce a solve_left() to a solve_right().)

SageMath seems to have trouble with sparse matrices with many columns (actually not that many.) The underlying issue seems to be that this already takes 6 seconds on my machine (and this is one order of magnitude smaller than the example from Peter.):

A = matrix(QQ, 1, 10000, range(10000), sparse=True)

If somebody wants to dig deeper, I am attaching some profiler output that can be opened with https://www.speedscope.app/.

julian