sprs
sprs copied to clipboard
sparsemat
Implement alga traits
As motivated in this thread on Krylov subspace solvers, I think it may be a good idea to try and see, whether this library's types could implement the traits provided by the alga library.
This would, of course, require some work. I would be happy to do that experiment myself, if you think it makes sense – at least as an optional feature.
TODO
- [ ] By-value operators (Add, Sub)
- [ ] Assignment operators (AddAssign, SubAssign)
- [x] Universal zero issue
- [x] Additive group for vectors
- [ ] Scalar multiplication for vectors
- [ ] Module & vector space traits
- [ ] Everything for matrices…
Hello,
As I mentionned in the users.rust-lang.org thread, I'm not sure the traits provided by alga exactly match the requirements for Krylov subspace methods, but so far I've only looked at the matrix trait, so I could be very wrong, and you're very welcome to show me where I'm incorrect ;)
Still, implementing some interoperability trait seems like a good idea on its own, especially if it's done under an optional feature.
You might encounter some difficulties in the implementation though, from a quick glance scalars in alga require only Clone, but there are some spots in sprs where I relied on Copy for convenience. But changing that seems like a nice move too.
Hey, there! Thank you for the quick feedback :)
First of all, I completely agree that alga is not an exact fit of what we need for that specific problem. However, I do think that a proliferation of interoperability traits for each specific algorithm (or class of algorithms) may be a larger problem maintenance-wise than this slight mismatch. Thus, I would rather have a set of general algebraic abstractions. Alga does provide those to a certain extent. And I suppose it can be refined further, if that turns out to be necessary. After all, this use case perfectly fits its design goal, as far as I can tell.
For what I am doing practically, the type signature for my BiCGSTAB algorithm looks like this:
pub fn bicgstab<V, F>(guess: V, vector: V, matrix: F, tolerance: V::Field) -> V
where V: FiniteDimVectorSpace + NormedSpace + Clone,
V::Field: Identity<Multiplicative> + PartialOrd + Clone,
F: Fn(&V) -> V,
I actually ended up not using the Matrix trait at all, as it seemed more complicated than what I needed. All the algorithm needs is the transformation associated with the matrix. Thus, I decided to pass that in as a Fn(&V) -> V instead of a matrix.
Nonetheless, the vector space and field abstractions are extremely useful. They reduce the amount of generics boilerplate, as they contain all the vector space operations. So from my perspective, implementing those vector abstractions takes priority.
Afterwards it might make sense to find out, in how far the matrix abstractions and sprs could be made to work together. For the Krylov subspace method use case that issue is secondary though, if it is at all necessary.
It looks like the vector space abstraction is quite useful to abstract over dense vectors, that's something useful to most linalg solvers, not just Krylov subspace methods, so that clearly looks interesting.
Going with matrix as a simple function looks really nice, though I'm a bit surprised that you don't need a bit more information, eg the number of rows and colums.
I'd see some trait like:
pub trait LinalgOperator
{
type V; // should be constrained to be in a vector space
fn rows() -> usize;
fn cols() -> usize;
/// Matrix vector multiplication
fn mat_vec(&self, rhs: &V, out: &mut V);
// add more methods for transposed mat vec product, mat mat product, with
// a default implementation that does nothing, and add some methods to check
// the availability of such operations
}
I see you're concerned that this could lead to a proliferation of specific traits, but experience with scipy's LinearOperator shows that it's a very useful abstratction. Maybe that could make a good addition to alga if that proves useful?
You'll notice that the mat_vec method takes an output parameter, I think it's better in an iterative solver to avoid constantly re-allocating estimates.
I am sure the simple function approach does not work in general. The transpose will be required in many algorithms. But knowledge of the row and column numbers is unnecessary, as the algebraic operations encapsulate that quite well.
Comparing Matrix and LinalgOperator there are row, get, column and get_unchecked methods required by alga, that do not have a corresponding concept in your suggested trait. That should be easy enough to separate via a subtrait relation. Nonetheless, would it be a problem to implement these methods in sprs?
As a conclusion, I am going to start with the vector space integration behind a feature flag, if you have no further objection. It will take a few steps of PRs to get there.
Indeed, a subtrait relation looks like a good idea. Most of alga's matrix trait should be straightforward to implement, mostly methods to get a row or column may be tricky with compressed matrices.
I have no objection, please try, and do not hesitate to ask any question or raise any issue with sprs' internals along the way ;)
So, while looking into this, I found a bit of show-stopper for alga integration: the sprs vectors have too many zeros. A vector space has exactly one zero, but sprs vectors are actually infinitely many vector spaces combined, as differently sized vectors are still the same type. Without adding statically typed sizes, we can't implement a proper zero trait. And that comes with a lot of downsides as well.
Perhaps, alga's abstractions are too mathematically strict to be applied in this context. Hm… maybe we should have laxer abstractions for linear algebra that can cope with this after all.
Indeed, a statically typed size would be far too inconvenient. Looks like the same problem should arise for dense vectors, how does nalgebra handle this?
Hm… what would you think about allowing operations on vectors of different sizes?
Such that at least this would work:
CsVec::new(0, vec![], vec![]) + CsVec::new(10, vec![…], vec![…])
Maybe even this:
CsVec::new(5, vec![…], vec![…]) + CsVec::new(10, vec![…], vec![…])
Semantics would be to add as many zeros as necessary. The first special case would be fairly easy to implement as a special case, just to have some vector that serves the role of a universal zero. It would definitely solve this problem without creating a new set of abstractions.
Whether the general change should be made, is a more complex question. It would probably also allow certain kinds of errors to go unnoticed. But I do not think it would be required to solve this problem.
To add to the general idea of allowing operations on differently-sized vectors, here is a relevant exchange on the mathematical background. Basically, dynamically-sized vectors are isomorphic to polynomials as a vector space. So it would even be correct to turn them into a real vector space that way. Again, whether that is user-friendly is a different question.
This isomorphism to polynomials does make it sound, but I really fear that having no dimension checks would catch less errors. I'd really like to have an example where this feature would be really useful to get a better sense of its utility.
I may be over-worried, it's quite rare to directly manipulate and sum sparse vectors in my experience. But that would be a surprising behaviour.
Okay, that makes sense. I will look into how feasible the more conservative solution is then. That is, to use a zero-sized vector as a universal zero.