-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[WIP] Optimizing {+,-,*} for structured matrices #28883
[WIP] Optimizing {+,-,*} for structured matrices #28883
Conversation
Thank you!!! How can I/others best help with this? |
@chriscoey There are a few places other than what I have already started (which I encourage people to check out and make more elegant/efficient!). The two that come to mind the most is just checking out other combinations of structured matrix types that aren't covered here and seeing if there is a way to speed them up and working with Adjoint/Transpose matrices. For the former, it seems like just a manual process of trying all the combinations and making sure they output the best type possible and do it in the best way possible (for multiply, the specialized methods are pretty efficient but there is some work that can be done on +/-). As for Adjoint/Transpose, I have precisely 0 experience working with those in Julia so any help at all on that front is greatly appreciated. The other thing that I briefly mentioned above is removing the matrix multiplication methods from SparseArrays. There are a few combinations that output sparse but unstructured matrices (at least not structured like anything supported in the stdlib, the ones that come to mind are things like upper bidiag * tridiag, which would be a great candidate for BandedMatrices but is currently returning a sparse matrix). Currently, the multiplication methods for these are in SparseArrays, not LinearAlgebra. Someone suggested having a fallback method inside of LinearAlgebra that outputs a suboptimal type and then having it overwritten in SparseArrays which seems reasonable but I don't know how the core contributors feel about that. One upside of this approach comes from the discussion on Discourse about removing SparseArrays from the stdlib: https://discourse.julialang.org/t/future-of-sparse-matrices-in-base-stdlib/11914 Looking forward to working with anyone who is interested! |
…nto optimize_binop_structured_matrices
Thanks! And what about ldiv/rdiv methods for structured matrices? I opened up an issue on that yesterday #28864 |
Ah right, in my head that was a totally different beast. If you have some progress there and want to combine it into this one that would be fine. Or they can proceed independently. Don't the solvers/div methods usually end up calling C or Fortran code? I thought most aren't implemented in pure Julia unlike the + or * implementations. |
I have been looking at div with a triangular matrix, and there are already pure Julia implementations in use for that (they have been there for years), eg julia/stdlib/LinearAlgebra/src/triangular.jl Line 1348 in d55b044
julia/stdlib/LinearAlgebra/src/triangular.jl Line 1200 in d55b044
so I think there is as much to gain from structured div methods as for the +,-,* methods you listed. and it is probably worth considering them at the same time. can anyone else chime in to agree/disagree? |
Thanks for the info. In that case, I think they should be done at the same time. If you have some progress done on these we can continue from there. |
… from binops of special.jl so that broadcasting takes over. Cleaned up some of the tests for special.jl
…tridiag, symtridiag, diag)
…These should go in another PR so this one can be merged more quickly. Revert "added sparse multiplication and division for triangular matrices. Fix JuliaLang#28451" This reverts commit 11c1d1d.
So as to keep this PR simple, I tried to remove the commits related to sparse matrix multiplication and division by @KlausC but I messed it up. I could use a bit of help with it. I think division and sparse matrix operations should be in a separate (maybe more than one) PR. The benchmarks have been updated and all of the regressions for structured +/- have been fixed. There are a few more small optimizations that can be made (making special +/- for some of the triangular matrices, similar to how it is done for multiplication), but they aren't major performance increases. |
Hi Marco,
ther is already PR #28507 under review, which includes only the changes for sparse strangular matrix multiplications. As it changes only one file and adds few new methods, it should be easy to back out the changes made there.
Klaus
|
@andreasnoack @Sacha0 Sorry to ping you but I believe this is ready for review. There is one part in particular that I think may have a better solution which I discuss in #28883 (comment) Thanks! |
I have a deadline tomorrow but will try to find time to review later this week. |
This was discussed in Slack but there is a class of errors in this PR that are related to JuliaLang/LinearAlgebra.jl#562.
so this is not ready for review just yet. Note that this also happens in v1. |
…ing diagonals are of different types. This still fails when the representation is a range and we get a step size of 0
…n eltype <: AbstractArray See PR 27289
Sorry for the long delay here. Great work. |
I forgot to add that I don't think that the error you mention should hold this back further as they are already errors. I'd prefer to get this PR in and then proceed from there. |
This is a work in progress to get binary operations on structured matrix types to work as intended. Currently there are several inefficiencies (usually due to fall backs or, the algorithm chosen to perform the binary operation, or some undesirable conversion mid operation) or missing/broken methods.
Current Status:
In the tables below, ❌ is for errors (where the operation is not even evaluated) and ✔️ are for problems that have a PR open to fix them.
First, we look at the current state of matrix addition and multiplication:
Below is a reduced table, where only things to be fixed are listed:
These errors can divided into roughly 3 categories (I am working on a list of examples for each):
Type 1 often happens when a catch all fallback is hit. Many of the multiplication cases have very efficient implementations (thanks to #15505) but convert to Array or SparseMatrixCSC somewhere which eliminates all performance. See bidiag.jl for most of those.
Type 2 usually is with very highly structured upper/lower matrices. I am compiling a list of these, but I noticed some in addition for diagonal-ish matrices.
Type 3 seems to be the most pressing, since the others still produce correct results but not as quickly. These are often found when we have symmetric types. For example:
Progress:
Currently, I have three related open PRs: #28343, #28345, #28534. My hope is that this PR subsumes them. I also have a patch that fixes most of the multiplication issues, but it fails for a symmetric matrices when the eltype is BigFloat/BigInt. This seems to be a problem with the matrix constructors and calling undefined values.
Related Issues and PRs:
#28507, #28534, #28343, #28345, JuliaLang/LinearAlgebra.jl#548, #27405, JuliaLang/LinearAlgebra.jl#533, JuliaLang/LinearAlgebra.jl#515, JuliaLang/LinearAlgebra.jl#525, JuliaLang/LinearAlgebra.jl#136, #28451, #29045
Benchmarks:
Matrices are
1000x1000
matrices ofFloat64
.To do: