# Non-sorted SparseMatrixCSC

**URL:** https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133
**Category:** Internals & Design
**Created:** [April 6, 2020, 11:21pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133 "2020-04-06T23:21:33Z")
**Posts on this page:** 4
**Page:** 2

<div class="post-metadata">

### Author: ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [August 28, 2020, 3:30pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/21 "2020-08-28T15:30:10Z")

</div>

Ok, this thread kind of went on in my head and I felt that one possibility was still missing…

I did some experiments with preventing setting zero entries in `jacobian.jl` . As a result, we get:  
 ![compare-prevent-zero-insertion](https://global.discourse-cdn.com/julialang/original/3X/0/f/0fd0ba3d240680bb18ac29e4764f2eefd6ac8cdf.png)

So this gives an improvement of assembly times by a constant factor, but we still see the O(N^2) behavior. This appears to be due to the fact that without further a priori information, ForwardDiff assumes that the jacobian is a full matrix and _calculates_ all these zeros, still leading to O(N^2) complexity.

For a known sparsity pattern (as we have for this case), the remedy [suggested](https://github.com/JuliaDiff/SparseDiffTools.jl) is to provide the sparsity pattern along with coloring information such that derivative calculations can be scheduled to handle only the nonzero entry cases.  
So instead of

```julia
ForwardDiff.jacobian!(jac, ffd, X, Y);

```

for `jac` with unknown sparsity pattern we do for `jac` with known sparsity pattern and nonzero entries:

```julia
using SparseDiffTools
...
colors = matrix_colors(jac)
forwarddiff_color_jacobian!(jac, ffd, X, colorvec = colors)

```

(for unpatched ForwardDiff) resulting in:  
 ![compare-colored-sparse](https://global.discourse-cdn.com/julialang/original/3X/b/2/b2b720a0c1e24db97882c02d0d69b4949cc6874e.png)

As we see, we now get O(N) behavior because now known zero jacobian entries are not calculated to begin with.

So we have two ways to handle the sparsity here: matrix coloring (once the sparsity pattern is known) or atomic assembly (which with ExtendableSparse performs reasonably well without a priori information about the sparsity pattern).

---

<div class="post-metadata">

### Author: ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [August 28, 2020, 9:24pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/22 "2020-08-28T21:24:38Z")

</div>

And then there is [SparsityDetection.jl](https://github.com/SciML/SparsityDetection.jl) which can detect the sparsity pattern automatically:

```julia
using SparsityDetection
input = rand(n)
output = similar(input)
sparsity_pattern = jacobian_sparsity(ffd,output,input)
jac = Float64.(sparse(sparsity_pattern))
colors = matrix_colors(jac)
forwarddiff_color_jacobian!(jac, ffd, X, colorvec = colors)

```

resulting in  
 ![compare-detected-colored-sparse](https://global.discourse-cdn.com/julialang/original/3X/c/7/c754af9feb72b4d7b10d46cbd6028d473f722835.png)

According to the test, the sparsity detection algorithm implemented exhibits O(N) behavior. In this example it is beaten by the atomic assembly. However, e.g. in an implementation of Newton’s algorithm and in many other cases one would perform the sparsity detection and matrix coloring only once and re-use the results during the computation of series of jacobians, so the relative overhead from the detection step would decrease.

I find the O(N) behavior of the sparsity detection step quite amazing. The algorithm behind this is described [here](https://openreview.net/pdf?id=rJlPdcY38B) and uses more deep Julia stuff…

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 29, 2020, 3:25am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/23 "2020-08-29T03:25:05Z")

</div>

Yeah our sparse automatic differentiation compiler tools are still pretty young but they should be pretty good for serious use now. That’s cool you did a verification of it on a big problem. It would be nice to have the code you’re testing on to turn this into a blog post.

---

<div class="post-metadata">

### Author: ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [September 1, 2020, 8:04pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/24 "2020-09-01T20:04:27Z")

</div>

Ok I will try to find some time for this in a week or two.

I also did some preliminary tests for 2D and saw the same O(N) behaviour. I assume there will be ways to beat the constant as well…

In fact I am thinking to integrate this into VoronoiFVM.jl to handle sparse couplings in systems of n PDEs - for flux and reaction terms I currently locally get full n\times 2n resp. n\times n matrices which often are sparse and test for zeros before putting values into the big matrix - this case could take advantage already of the current implementation.

[Previous page](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133.md?page=1)
