# Tackling the AbstractArray fallback trap - wrapped sparse matrices

**URL:** <https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592>\
**Category:** Internals & Design\
**Tags:** proposal, linearalgebra\
**Created:** [December 12, 2018, 11:35am UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592 "2018-12-12T11:35:49Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [December 12, 2018, 11:35am UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/1 "2018-12-12T11:35:49Z")

</div>

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](https://discourse.julialang.org/t/slow-sparse-matrix-vector-product-with-symmetric-matrices/3968) or [2](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431).  
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.

1. provide a unified way to transform any kind of wrapped sparse matrix to a `SparseMatrixCSC`
2. provide specialized algorithms for selected wrapped types
3. establish a “sparse matrix fallback”, which replaces the “abstract arrax fallback” for wrapped sparse types
4. 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](https://github.com/KlausC/SparseWrappers.jl).  
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](https://github.com/JuliaLang/julia/pull/30018).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 12, 2018, 12:13pm UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/2 "2018-12-12T12:13:35Z")

</div>

There are two not easily separable concepts here: storage (dense, sparse, block), and matrix properties (diagonal, unit lower triangular, symmetric, hermitian). Currently, wrappers like `Symmetric` accept all abstract matrices (that satisfy some properties, eg are square), and just impose symmetry by just using the elements above/below the diagonal, but it is my impression that this works ideally only when the underlying storage is dense.

Moreover, a lot of LinearAlgebra routines use modifying versions (eg `mul!` as a fallback), which do not necessarily mesh well with anything but dense matrices (and some similar layouts).

I run into these things all the time, eg [#28869](https://github.com/JuliaLang/julia/issues/28869), but I don’t know what a good solution would be. Perhaps a disjoint set of types instead of wrappers, and _traits_ to query matrix properties (eg a matrix can be diagonal, lower triangular, upper triangular, and sparse at the same time), which would allow some (possibly not 100% optimal) fallbacks until all cases are implemented.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [December 12, 2018, 12:56pm UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/3 "2018-12-12T12:56:33Z")

</div>

> [@Tamas\_Papp](#):
>
> lot of LinearAlgebra routines use modifying versions (eg `mul!` as a fallback), which do not necessarily mesh well with anything but dense matrices

I do not agree. The `mul!(::StridedArray, X, ::StridedArray)` are a counter example, where the sparsity of `X` can be used efficiently. Implementations for `Transpose(sparse)`, `UpperTriangular(sparse)`, `Transpose(LowerTriangular(sparse))` etc. are already merged in master, while `Symmetric(sparse)` etc. exist as reviewable PR.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [December 12, 2018, 12:59pm UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/4 "2018-12-12T12:59:27Z")

</div>

> [@Tamas\_Papp](#):
>
> I run into these things all the time, eg [#28869](https://github.com/JuliaLang/julia/issues/28869)

You get

```julia
julia> using SparseWrappers

julia> A = lower_triangular(Diagonal([1, 2]))
2×2 Diagonal{Int64,Array{Int64,1}}:
 1 ⋅
 ⋅ 2

julia> A * A'
2×2 Diagonal{Int64,Array{Int64,1}}:
 1 ⋅
 ⋅ 4

```

That is a point for item 4. in the original post.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [December 12, 2018, 1:04pm UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/5 "2018-12-12T13:04:13Z")

</div>

> [@Tamas\_Papp](#):
>
> two not easily separable concepts here: storage (dense, sparse, block), and matrix properties

That is true. My statement is not, that it is easy, but that it is possible and has been done.

> [@Tamas\_Papp](#):
>
> … but it is my impression that this works ideally only when the underlying storage is dense.

When the underlying storage is sparse, there is a gap in the current release of `Julia`. This proposal is going to close this gap.

---

<div class="post-metadata">

**Author:** ![klacru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klacru/32/27890_2.png) [@klacru](https://discourse.julialang.org/u/klacru)\
**Post date:** [January 5, 2023, 6:01pm UTC](https://discourse.julialang.org/t/tackling-the-abstractarray-fallback-trap-wrapped-sparse-matrices/18592/6 "2023-01-05T18:01:43Z")

</div>

There is now a new traits based idea for this problem. See [README](https://github.com/KlausC/AbsArrFallbackTrap.jl).  
Comments are welcome.