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 `AbstractArray`

or `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
`SparseMatrixCSC`

- 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 `Matrix`

or `SparseMatrixCSC`

depending on the underlying storage of `A`

. If `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 `SubArray`

s, and methods to transform transposed and adjoint wrappers of `Matrix`

and `SparseMatrixCSC`

.

Missing are methods for symmetric and hermitian matrices, for `Unit[Upper|Lower]Triangular`

, for `Diagonal`

, `Bidiagonal`

, etc. and for nested wrappers, e.g. `Adjoint(view(sparse(...)))`

.

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 (`[]`

, `getindex`

, `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 `transpose`

and `adjoint`

for this purpose. Also most constructors avoid to deliver duplicate things like `Symmetric(Symmetric(A))`

. Also `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 `SubArray`

and `ReshapedArray`

, we can reduce the potentially infinite number of combinations of wrappers to finite. So I propose new quasi-constructors like `symmetric`

, `upper_triangular`

, etc. besides `transpose`

and `adjoint`

.

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 (`Symmetric.uplo`

…).

In order to make these improvements available to the public, I would like to move them piece by piece to `stdlib/SparseArrays`

and `stdlib/LinearAlgebra`

. There seems to be a review bottleneck though, e.g. #30018.