# Performance of lazy wrappers applied to sparse matrices

**URL:** <https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431>\
**Category:** Internals & Design\
**Tags:** linearalgebra\
**Created:** [September 24, 2018, 5:38pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431 "2018-09-24T17:38:51Z")\
**Posts on this page:** 16\
**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:** [September 24, 2018, 5:38pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/1 "2018-09-24T17:38:51Z")

</div>

For motivation please look at the following suite of benchmarks:

```julia
julia> m, n , nz = 10000, 10000, 1000000
(10000, 10000, 1000000)

julia> A = spzeros(m, n)
10000×10000 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> @btime A' == A'
  2.415 s (4 allocations: 64 bytes)
true

julia> A = sprand(ComplexF64, m, n, nz / m / n);
julia> @btime A' == A'
  7.932 s (4 allocations: 64 bytes)
true

julia> x1 = A';
julia> x2 = (copy(A)');
julia> @btime x1 == x2
  12.005 s (2 allocations: 32 bytes)
true

julia> @btime x1.parent == x2.parent
  3.141 ms (0 allocations: 0 bytes)
true

```

Maybe you can share my feeling, that only the last figure is acceptable (0.003 seconds vs. 12 seconds).  
As `A' === Adjoint(A)`, there seems to be an issue with `Adjoint` and `==`. Nevertheless I think, the issue is with all of those wrappers, which were introduced to improve performance, namely:  
`Adjoint, Transpose, UpperTriangular, LowerTriangular, Symmetric, Hermitian`.  
And not only the `==` is problematic.

```julia
julia> @btime x = abs.(UpperTriangular(spzeros(m, n))');
  776.993 ms (8 allocations: 763.02 MiB)

```

is an example of combinations of those wrappers, which are not taking advantage of the sparsity of the underlying matrix.

My question is: are there efforts being undertaken to improve the performance of lazily wrapped sparse Matrices in general and specially of combined wrappers? I would like to start a project on this, but want first to see, if somebody else is working in this area.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [September 24, 2018, 6:40pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/2 "2018-09-24T18:40:56Z")

</div>

This seems like it’s just a missing method. Could you open an issue about this?

---

<div class="post-metadata">

**Author:** ![Stephen\_Vavasis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephen_vavasis/32/3389_2.png) [@Stephen\_Vavasis](https://discourse.julialang.org/u/Stephen_Vavasis)\
**Post date:** [September 25, 2018, 12:44am UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/3 "2018-09-25T00:44:40Z")

</div>

I already opened an issue on this. I would say that the sparse matrix library needs a lot of work.

[https://github.com/JuliaLang/julia/issues/28432](https://github.com/JuliaLang/julia/issues/28432)

There are some other related issues:

[https://github.com/JuliaLang/julia/issues/29353](https://github.com/JuliaLang/julia/issues/29353)

[https://github.com/JuliaLang/julia/issues/28948](https://github.com/JuliaLang/julia/issues/28948)

---

<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:** [September 25, 2018, 3:58am UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/4 "2018-09-25T03:58:23Z")

</div>

I fear there is not only one method missing.

I agree. The `==` example was just a motivation. Actually it is a pain to work with sparse arithmetic; my project stalled, because after a Q-R decomposition, multiplying `R * x` with a dense vector `x` did not finish in a reasonable time. I made a PR about that issue some time ago [https://github.com/JuliaLang/julia/issues/28451](https://github.com/JuliaLang/julia/issues/28451), which will be reviewed soon, I hope.

I want to try to consolidate those related issues like yours. A similar effort is this one: [https://github.com/JuliaLang/julia/pull/28883](https://github.com/JuliaLang/julia/pull/28883)

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [September 25, 2018, 8:11am UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/5 "2018-09-25T08:11:09Z")

</div>

An unfinished effort to provide a relatively general fix can be found in [https://github.com/timholy/ArrayIteration.jl](https://github.com/timholy/ArrayIteration.jl). The strategy there is to create iterator types that allow you to synchronize your access of relevant entries in a matrix pair, and thereby write generic algorithms that only visit entries that need it. It’s unknown whether algorithms based on such iterators can be written without _any_ generalization penalty (indeed, that seems a bit unlikely to me), but I expect that it should be possible to provide something that’s at least a pretty good fallback definition for almost any combination of types.

The reason this makes sense is that if there are `N` specialized matrix types, even binary methods (e.g., `==`, `*`, `+`) result in `N choose 2` method specializations you need to write. Given that `N` is not small (and you’ve left out `Diagonal`, `Bidiagonal`, `Tridiagonal`, `SymTridiagonal`, `UnitLowerTriangular`, `UnitUpperTriangular`, and packages like `BandedMatrices`), this means a lot of specializations. In contrast, writing an iterator implementation for each matrix is an `O(N)` problem rather than an `O(N^2)` problem, and thus makes the problem far more manageable.

That package is not currently at a high point in my own priority list (which is quite long), so if folks are interested I’d be happy to see someone take ownership of it. Two notes:

- if memory serves you want to start from the `teh/stored2` branch rather than `master` (but I could be wrong)
- it’s possible this will only make sense once we can stack-allocate objects that contain a reference to the heap (see [performance - Unexpected memory allocation when using array views (julia) - Stack Overflow](https://stackoverflow.com/questions/47590839/unexpected-memory-allocation-when-using-array-views-julia/47607539#47607539)), because this package creates lots of wrappers and you don’t want that to allocate. However, I am hopeful that the improved optimizer we now have would elide the wrapper creations in most cases, so perhaps with careful work it won’t be necessary to wait.

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [September 25, 2018, 9:09am UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/6 "2018-09-25T09:09:38Z")

</div>

copy seems to help a lot for transpose

```julia
@time x' == x'
4.18 sec
@time copy(x') == copy(x')
0.00026 sec
```

---

<div class="post-metadata">

**Author:** ![Stephen\_Vavasis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephen_vavasis/32/3389_2.png) [@Stephen\_Vavasis](https://discourse.julialang.org/u/Stephen_Vavasis)\
**Post date:** [September 25, 2018, 1:02pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/7 "2018-09-25T13:02:41Z")

</div>

Is there anybody “in charge” of the sparse matrix library? Fixing all these missing methods with a new iteration protocol is not a small patch-- this requires a coordinated effort of several people.

---

<div class="post-metadata">

**Author:** ![Stephen\_Vavasis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephen_vavasis/32/3389_2.png) [@Stephen\_Vavasis](https://discourse.julialang.org/u/Stephen_Vavasis)\
**Post date:** [September 25, 2018, 1:14pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/8 "2018-09-25T13:14:41Z")

</div>

This pull request represents a very impressive effort! But it’s not clear to me-- are you attempting to fill in all triples (left operand, right operand, operator) of missing methods, or are you adopting the approach suggested by Tim Holy and others to institute a new iteration protocol in order to reduce the combinatorial explosion at the cost of some performance?

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [September 25, 2018, 1:18pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/9 "2018-09-25T13:18:28Z")

</div>

> Is there anybody “in charge” of the sparse matrix library?

Anyone who does the work. `git` shows that 29 people contributed to the SparseArrays stdlib in 2018 so far.

> Fixing all these missing methods with a new iteration protocol is not a small patch-- this requires a coordinated effort of several people.

Agreed, it’s not a small job, which is why it’s not done already. But it does easily break into discrete tasks. The biggest task is getting the infrastructure in ArrayIteration working properly, probably using a couple of array types as test beds. The second biggest task is to write iterators for all the (remaining) array types. The final task, writing a half-dozen generic methods, will be pretty simple by comparison: `==` might be less than a dozen lines, for example.

---

<div class="post-metadata">

**Author:** ![Stephen\_Vavasis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephen_vavasis/32/3389_2.png) [@Stephen\_Vavasis](https://discourse.julialang.org/u/Stephen_Vavasis)\
**Post date:** [September 25, 2018, 1:23pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/10 "2018-09-25T13:23:30Z")

</div>

Could you explain in more detail why it’s necessary to solve the “small-item-on-the-stack” issue before you can implement array iteration? In particular, the new iteration protocol in 0.7 seems to work fine even though the small-item-on-the-stack issue is still open.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [September 25, 2018, 1:31pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/11 "2018-09-25T13:31:15Z")

</div>

Algorithms like matrix multiplication can often be compactly represented in terms of column & row iterators, e.g., [https://github.com/timholy/ArrayIteration.jl/blob/6736ec22025e0b0c09d948a6ec1cecec7691e7bc/src/linalg\_couple.jl#L4-L12](https://github.com/timholy/ArrayIteration.jl/blob/6736ec22025e0b0c09d948a6ec1cecec7691e7bc/src/linalg_couple.jl#L4-L12) (wow, I got farther implementing this than I had remembered…). Here’s an example of a column iterator: [https://github.com/timholy/ArrayIteration.jl/blob/6736ec22025e0b0c09d948a6ec1cecec7691e7bc/src/sparse.jl#L50-L53](https://github.com/timholy/ArrayIteration.jl/blob/6736ec22025e0b0c09d948a6ec1cecec7691e7bc/src/sparse.jl#L50-L53). It holds a reference to the original array, and hence if the compiler doesn’t inline & elide everything then it will be allocated on the heap.

It’s not at all crazy to hope that the compiler will inline and elide everything, so this might be doable today. But _if_ you hit a roadblock, my prediction is that’s where it will happen.

---

<div class="post-metadata">

**Author:** ![Stephen\_Vavasis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephen_vavasis/32/3389_2.png) [@Stephen\_Vavasis](https://discourse.julialang.org/u/Stephen_Vavasis)\
**Post date:** [September 25, 2018, 2:27pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/12 "2018-09-25T14:27:50Z")

</div>

Allow me to summarize this situation:

- You (Tim) are proposing a system of coupled iterators (simultaneous iteration over both operands). The system is elegant and powerful, a bit tricky to understand, and may require some additional core-compiler optimizations. However, it has the potential to solve almost every current performance issue in the sparse matrix library. It is not at the top of your priority list.

- Others (including me) have proposed straightforward one-operand iterators that could solve a subset of the issues, mainly the low-hanging fruit, but would fail to address more complex matrix-matrix operations.

- The PR 28883 mentioned earlier in this thread proposes and partly implements fully optimized individually written codes for binary operators (+, \*) for many specific pairs of matrix operands.

Each of these approaches has its pros and cons, but to some extent they are competing, which makes it difficult for someone who wants to jump in.

This situation cries out for a leadership decision! Someone with more clout than me in the Julia community needs to post an overarching framework for fixing the sparse matrix library so that the rest of us can make PRs to implement the framework.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [September 25, 2018, 2:47pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/13 "2018-09-25T14:47:05Z")

</div>

That’s a great summary. From a logical standpoint these might appear to be competing strategies, but honestly I think in real-world development that’s unlikely. Suppose I or someone implemented my ambitious strategy: wouldn’t you want to benchmark it against something? Wouldn’t that something be a more straightforward implementation? If that more straightforward implementation already exists, doesn’t that save the “ambitious developer” a ton of time, not to mention the fact that meanwhile other users/developers have been enjoying the benefits of something that already works well in a subset of cases?

Consequently as usual I think leadership on this issue should come “from the bottom,” meaning the people who do the work. It seems that a number of good solutions are just around the corner, and I’d urge folks to push ahead with those and get them merged. Over the longer timescale one or more people might tackle the more ambitious approach, and then benchmark against what we have now. In a few places it’s possible the general strategy will be just as good as the handwritten one and can replace it, but when that’s not true then one can keep the hyper-optimized version and use the more generic fallback to handle any missing cases.

---

<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:** [September 25, 2018, 3:30pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/14 "2018-09-25T15:30:06Z")

</div>

If you are speaking about [8451](https://github.com/JuliaLang/julia/issues/28451), that is exclusively the multiplication and division of dense RHS by triangular sparse matrices. “triangular” means one of `UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular` literally and the `Transpose` and `Adjoint` of those. As for each matrix type multiplication and division is implemented, there are 24 cases to cover, which was implemented in 8 functions. In the pre-PR master branch only 6 of 24 are working at speed (only divisions, but not for `Unit...Triangular`).

---

<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:** [September 25, 2018, 4:14pm UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/15 "2018-09-25T16:14:41Z")

</div>

Inspired by your issue @Stephen_Vavasis  
[https://github.com/JuliaLang/julia/issues/28432#issuecomment-416672359](https://github.com/JuliaLang/julia/issues/28432#issuecomment-416672359)  
I just started the `nziterator` approach.

---

<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:** [October 1, 2018, 9:53am UTC](https://discourse.julialang.org/t/performance-of-lazy-wrappers-applied-to-sparse-matrices/15431/16 "2018-10-01T09:53:20Z")

</div>

[https://github.com/KlausC/SparseWrappers.jl](https://github.com/KlausC/SparseWrappers.jl)