sparsemat / sprs

sparse linear algebra library for rust
Apache License 2.0
381 stars 45 forks source link

added MulAcc<A,B> support for matrix csc & csr X dense vec formats #288

Closed experiment9123 closed 3 years ago

experiment9123 commented 3 years ago

2 functions that were easy to do. 1 minor change in each - flip the mul_acc operands to be matrix elem, vec elem (the order you see it in the expression Matrix X Vector) - and the tests ended up needing a couple of type annotations. besides that the callgraph is unaffected, everywhere else those fn's are used.. nothing changed.

mulimoen commented 3 years ago

I would like to see more of the code converted, just to see if type inference still works. Could you implement equally for all the smmp code too? This would make porting Mul<CsMat> trivial, so you could add this one as well.

experiment9123 commented 3 years ago

now I tried converting some more, and there is one thing where the inference starts asking for ascriptions here and there. in some of the tests, the result type is not inferable from the inside; hence I had to modify them eg "assert_eq!({let x:f32=blah;x},expected). I think this happened because these tests dont actually use or return the values.. they just calculate and assert.

the 'examples' code is unchanged.

Nonetheless if these needed ascriptions seem like a PITA, one idea would be to rename these fully generic "dot" "dot_infered" or something, and have a second "dot" (wrapper) with the original type assumptions( which then just calls dot_inferred)

Take a look at the measures I had to take in the latest version 7d7b923, and see what you think.

for the record, I think csr_mul_csvec which I just did is the inner function (called by Mul<...>??) that is needed for my use case, so we might be pretty close here

experiment9123 commented 3 years ago

I did another little experiment (not pushed) where I made a helper trait A:MulAddWith<B>{ type Accumulator } .. with impl such that A:MulAddWith<B, Accumulator=N> - allows binary operators to express equiv of c++ typedef decltype(A()B()+A()B()) N , "N=type_of(AB+AB)"

one fly in the ointment with that is the default impls might have conflicted with autoimpling MulAdd ->MulAcc1 (i think) i.e. once you have bothA:MulAddWith, N : MulAcc<A,B> yields the conflict with MulAdd

if this is the case, one workaround might be to manually impl the MulAccs using MulAdds for specific types in a macro (the primitive ints and floats)

i'm trying to think if I could build the MulAdd requirement into the MulAddWith trait .. but in the meantime i've converted the other easy functions

experiment9123 commented 3 years ago

right in the latest commit, i've extracted the Operator guts into a fully generic named function pub fn csmat_mul_csmat,(and genericised the smmp.rs fns). The Mul operator is left single-type. This free function means the functionality is available*, but less sadly not so discoverable (yet).

I'm uneasy about the idea of using Mul, Add<>::Output's to compute the types "forward" for the operator because these aren't related to a generic impl for MulAcc<A,B>. it's going to cause confusion between MulAcc, MulAdd, and simple A*B+A*B. I like the current design where so long as you impl MulAcc, everything works and the 'temporary products' needn't be exposed.

perhaps we can get the functionality and internals all working fine before we attack ideas like Lazy Expression templates ( I would also suggest that methods like mat.mul_mat(&rhs), mat.mul_vec(&rhs) could be almost as discoverable as the operator, when using IDEs) or whatever else it might take to get the operators working.

mulimoen commented 3 years ago

The extra type annotations are a bit worrying, and would make the code harder to use for first time users. Adding some free-standing functions which uses MulAcc internally would be a good first step, which is what we have in the smmp module. Maybe adding these accumulating methods on CsMat/CsVec is good enough for this specialized use, and would leave Mul for the scalar case?

mulimoen commented 3 years ago

For this case we must ensure functions requireing MulAcc takes care of the non-commutate operation. Are there no way to fix this as was done in #288?

experiment9123 commented 3 years ago

The extra type annotations are a bit worrying so far i've only seen it affect (i) the tests and (ii) benchmarks - which dont store or return anything .. i'm hoping it'll "just work" and these tests use the "internal" named functions. Currently all the operators still do N:MulAcc<N,N>, so their interface is unchanged.

Maybe adding these accumulating methods on CsMat/CsVec is good enough for this specialized use

I would certainly be happy with that. what I'm also seeing is that even "returning" named functions infer their types better. I think it's that the associated types in the operators is just slightly different to having a type-param list on the function, even though there's no difference in principle (eg the new func csmat_mul_csmat infers N,A,B, just fine). it should just be like "a function that happens to be called *".. but isn't exactly.

For this case we must ensure functions requireing MulAcc takes care of the non-commutate operation. Are there no way to fix this as was done in #288?

what I've seen of the code so far .. it should definitely be possible but it's intense to figure out the detail, i'll need to take a breather: the issue is the matrix ops flip the args and call multiples with transposed versions of themselves, and I think there's several places I'd need to add and propogate some kind of argflip lambda.

after a pause and learning abit more about this source base , new approaches might click .

If the layouts were compile-time known (like explicit struct SparseMatrixCsc, SparseMatrixCsr etc) the argflip problem goes away. but I'm guessing it does this because some of those permuations ("row vs col major" x "csc/csr") are just equivalent to eachother (what I have in my head is seperate single layout types, with simpler Mul overload implementations, then an uber-type which wraps them in an enum)

in the use case I personally have in mind, it's going to be a known format (csr or csc) X sparse vector - there's probably enough ported already to handle this now (I can call the format-specific functions).

A couple of motivations here.. the "experiment I had in mind" is pretty much doable now, but I also care about the "reputation of the language & ecosystem". compare the capabilities of Eigen/C++ with "the best sparse matrix lib I can find in rust" .. it would be nice to have something as versatile as possible.

experiment9123 commented 3 years ago

been paused on this for a couple of days.. in 2 minds about it if you merge this - you enable my use case (via fn csmat_mul_csmat) - but it does complicate the internals for the next person who tries to add a feature. The trait bounds get quite fiddly.

i'm increasingly tempted to just write the path+format i was after independently.. which is a shame as it doesn't feed back into the community as much as improving sprs. (it would take about 500-1000 lines max).

Haven't given up 100% though. This library could definitely do it, and it would be great to have something this versatile and compreheensive (to compete with C++ Eigen..)

mulimoen commented 3 years ago

If the layouts were compile-time known

I think this would be useful for reasoning about the cost of matrix operations, but it would possibly explode the number of Mul impls. How would this interact with CsMat, CsMatI, CsMatView, CsMatViewMut etc.?

experiment9123 commented 3 years ago

If the layouts were compile-time known

I think this would be useful for reasoning about the cost of matrix operations, but it would possibly explode the number of Mul impls. How would this interact with CsMat, CsMatI, CsMatView, CsMatViewMut etc.?

How it would interact with all those types - I dont even know my way around them fully enough to definitively answer that. You're right it's a concern.

I would not suggest supporting all of them - only a few permtations that have particular reason to exist (starting with my own particular use case) . I need to draw a roadmap ..

One possibility would be: [Step 1] - Operators only work with uniform types. Named functions for the multi-type generic matrix mul. Expose some internals for known type combinations,

[Step 2] expose a few known, easy combinations with operators for the multi-type cases. eg.. DenseVec X SparseMatrixCOO -> DenseVec has an easy impl. SparseVec X MatrixCSC -> DenseVec has an easy "Scatter" style multiply impl. DenseVec X MatixCSR has an easy "Gather" style multiply impl. ( for my use case, i''m guessing that SparseVec X MatrixCSC->Dense -> re-sparsify the result is worth trying)

[Step 3] if we fix "the param flip problem" the fully generic, all-types approach can be fully generalised.

experiment9123 commented 3 years ago

(just been reading through the requested changes, i'm ok with the suggestions, i'll have to figure out exactly what was happening with dot and the suggestion there .. are you suggesitng reverting some of the generics that break type inference for existing users, then adding new functions like "dst.dot_acc(&src1,&src2)" for full inference ? i'm fine with anything like that

mulimoen commented 3 years ago

Build failure is due to gitlab upgrade, will rerun and merge later on. Thanks @experiment9123 for this PR!

experiment9123 commented 3 years ago

Hi just letting you know I havent forgotten about this, i have paused whilst working on other things (still very much in rust). let me figure out how to grab your suggested changes. there is a risk that if i dont get this done this "round" of changes will diverge.

I agree it really wants the "lambda parameter flip" for the non-commutative behaviour for 100% versatility