When handling wrapped sparse matrices, for example
Symmetric(sparse(...)), it easily happens, the performance is not as expected, because you get caught in the “abstract array fallback” trap: while you are expecting to use an algorithm, which takes advantage of the sparse structure of your problem, a generic algorithm for argument type
AbstractMatrix takes control, while the same algorithm is available for
SparseMatrixCSC. This is not a new issue - see for example 1 or 2.
In my opinion it is also a severe competition drawback, compared to MATLAB, if we cannot handle really big sparse symmetric matrices in a better way.
I want to discuss a way out of the trouble.
- provide a unified way to transform any kind of wrapped sparse matrix to a
- provide specialized algorithms for selected wrapped types
- establish a “sparse matrix fallback”, which replaces the “abstract arrax fallback” for wrapped sparse types
- avoid deeply nested wrapped types by quasi constructors, which simplify combinations of wrappers
Ad 1: I would like to have something like
unwrap(A::AbstractArray), which returns an
SparseMatrixCSC depending on the underlying storage of
A has already the target type, return
A itself. Also nested wrappers must be supported.
Currently we have
tri[ul] to convert triangular matrices,
getindex to convert
SubArrays, and methods to transform transposed and adjoint wrappers of
Missing are methods for symmetric and hermitian matrices, for
Bidiagonal, etc. and for nested wrappers, e.g.
And there is currently no unique function to perform the conversions.
Ad 2: There has been already some effort in this area for a few common cases. Nevertheless there are some algorithms and there as quite a lot of useful (nested) wrappers. Potentially the number of nested wrappers is unlimited. So this approach can improve the situation at most important edges, but never be complete.
Ad 3: Some or all methods with argument type
AbstractMatrix, which implement an algorithm using the standard interface (
setindex) are affected. These methods obtain a new signature (rename or extra trait argument). At the same time a new method with the original signature is defined, which dispatches to the appropriate implementation method. The special cases of item 2. and the special implementations for
SparseMatrixCSC have to be shared by the dispatching process.
Ad 4: We have already
adjoint for this purpose. Also most constructors avoid to deliver duplicate things like
Hermitian(Symmetric(A)) = Hermitian(A), but these constructors do a half-heated job and
Symmetric(Hermitian(A)) = Symmetric(A) is not correct.
Also combinations with triangular matrices, and the structured types (
Diagonal, …) allow simplifications, which can reduce nesting depth. Actually, if we exclude
ReshapedArray, we can reduce the potentially infinite number of combinations of wrappers to finite. So I propose new quasi-constructors like
upper_triangular, etc. besides
I implemented most of those ideas in SparseWrappers.
It was not possible to implement the “sparse array fallback” pattern for real implementation, because that requires changes of the existing “abstract array fallback” implementations in
stdlib. The pattern is demonstrated in an example.
As a principle, all branches are evaluated at compile time by evaluating only type information. That was not possible at some areas where structural information is only stored in a data field (
In order to make these improvements available to the public, I would like to move them piece by piece to
stdlib/LinearAlgebra. There seems to be a review bottleneck though, e.g. #30018.