Sparse matrix implementation

[c410-f3r]

Hello,

Does anyone mind if I send a sparse matrix pull request? I am implement it and should be done in ~2 weeks.

Thanks

[AtheMathmo]

Hey!

Thanks for writing this issue. I would be really interested to see this but I'm not sure exactly what the requirements would be for it to be merged.

Here is some discussion on this issue on reddit.

I am not all that knowledgeable on sparse linear algebra and so will find it hard to recognize the quality of the implementation. With that said I think the following criteria are musts:

• Should be kept separate from existing Matrix/Slice/Mut implementation. I think at first these implementations shouldn't be affected at all.
• Should provide enough basic functionality to be usable in some applications (basic arithmetic and manipulation perhaps?)

If you are still up for implementing this I'll be very happy to read through the PR when it does land. I'll also be very happy to discuss your plan now if you would like some feedback in advance.

As one final note - even if it is decided that this shouldn't land in rulinalg I would be very excited to see it as a separate crate which rulinalg/rusty-machine can rely on.

[c410-f3r]

Thanks for the feedback. Currently I am not modifying any code but I think that we could have two implementations (dense and sparse) of a more abstract Matrix trait. There are many common corners that need to be discussed. Would be great to see runlinalg as a more complete linear algebra library but anyways, I learned a lot and I am doing my best starting with CSR, CSR slice and basic arithmetic. More questions are going to arrive soon :)

[Andlon]

It's great to hear that you are interested in implementing sparse linear algebra for rulinalg! I was part of the Reddit discussion that James linked to, and there I advocated for sparse linear algebra as a separate library. At the time I made that comment, rulinalg looked to be mostly a one-man library, and as such I figured it made more sense to focus primarily on dense linear algebra first. However, it seems it has picked up some traction and gained more contributors, so I wouldn't consider it a bad idea to investigate the implementation of sparse data structures and algorithms too.

As James pointed out though, I think it should be kept separate from the dense linear algebra. In fact, I can't recall ever having needed to mix sparse and dense matrices for any application, and the only link between them is the fact that they both operate on dense vectors. Hence, as far as I can see, it makes a lot of sense to simply delegate all sparse linear algebra code into its own module, without affecting the existing dense linear algebra code.

Figuring out the right style of API for sparse linear algebra can also grow to become very difficult, once you take into account various iterative algorithms and how to provide a somewhat uniform interface with regards to preconditioning, and on top of that how to make it all performant. However, that's a bit further into the future. I think we should rather make some baby steps, and instead be open to the idea that we'll likely have to iterate on different APIs and ideas quite a lot, until we get it right.

So for a first pass, I suggest avoiding the subject of iterative solvers for now, and make something usable with basic arithmetic. I suggest the following features for the start of an experimental sparse module:

• CSR matrix storage format
• CSRMatrix basic arithmetic (+, -)
• CSRMatrix matrix-vector multiplication
• CSRMatrix matrix-matrix multiplication
• In typical applications, you cannot efficiently construct a CSR matrix directly. One approach that is both flexible and reasonably performant is the so-called Coordinate List format. SciPy (scientific library for Python) has its own coo_matrix type which can be converted to CSR, while Eigen (C++ library) can construct a sparse matrix through a .setFromTriplets( ... ) function, taking in a vector of triplets. I can imagine that perhaps something like CSRMatrix::from_triplets(vector_of_triplets) would be idiomatic for Rust/rulinalg? Presumably, a triplet should be a trait that provides row, col indices as well as the value at the given position.
• I'm not sure exactly how the details of slicing should work. In this case we might need to involve CSC matrices for efficient slicing along the other axis. Would be interested to look at goals and potential problems here.

As for the solution of linear systems... While iterative solvers bring a lot of complication, direct solvers are also very challenging, in fact usually far more so than for dense matrices. As usual, Netlib is a decent starting point. I must admit that I don't know much about direct algorithms myself; I usually work with iterative ones.

Also, let's perhaps not worry if the arithmetic algorithms are well optimized in this first pass. Like I said, I think it's likely that a lot of changes to the API and internals must be made as more features are added. Let us perhaps instead focus on working implementations that are well backed up by tests which are easy to refactor?

[AtheMathmo]

Thank you @Andlon for such a detailed comment! I totally agree with everything you've said - and like last time I feel a lot better informed now too!

I agree that the goal should be working, well tested code that is easy to refactor and adjust.

[c410-f3r]

I just send a pull request so you guys can see (more or less) how the implementation looks like. It is not intended to be merged and there are much things left to do. If you guys are willing to include a sparse matrix implementation, I am committed to do my best and work hard for rulinalg.

[AtheMathmo]

@c410-f3r awesome! Thank you so much for getting that in. I'm really excited to check it out - though as I said before a lot of the details will be beyond me. Hopefully @Andlon will be around to provide some more technical feedback (if he doesn't mind).