# 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:** 20
**Page:** 1
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 6, 2020, 11:21pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/1 "2020-04-06T23:21:33Z")
Adding elements to a sparse one after another is slow due to sorting, allocating and copying. Is there another `AbstractSparseMatrix` that is more friendly to bit-by-bit construction. Or even one that can leverage the fact that elements are added in the correct (sorted) sequence.
My quick and dirty prototyping led me to a performance bottleneck that is appending to existing vectors. Can someone point me to a Julia vector data structure that grows in chunks, maybe even without the need to copy old elements?
---
### Author:  [@platawiec](https://discourse.julialang.org/u/platawiec)
#### Post date: [April 7, 2020, 1:05am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/2 "2020-04-07T01:05:55Z")
Can you show a little bit more about how you are appending these vectors?
`push!` and `append!` do grow the underlying vector storage in chunks. Perhaps `sizehint!` will fit your needs. It is difficult to say without more context.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 7, 2020, 1:16am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/3 "2020-04-07T01:16:44Z")
> [@MaximilianJHuber](#):
>
> Adding elements to a sparse one after another is slow due to sorting, allocating and copying. Is there another `AbstractSparseMatrix` that is more friendly to bit-by-bit construction.
Assuming you don’t need to _use_ the matrix until it is fully constructed, then you could add your data to `(I,J,V)` arrays element-by-element in an arbitrary order and then call the `sparse(I, J, V)` constructor.
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 7, 2020, 1:32am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/4 "2020-04-07T01:32:58Z")
> Assuming you don’t need to _use_ the matrix until it is fully constructed, then you could add your data to `(I,J,V)` arrays element-by-element in an arbitrary order and then call the `sparse(I, J, V)` constructor.
That was my plan.
@platawiec thanks for the advice:
```julia
using SparseArrays
mutable struct MySparseMatrix{Tv,Ti<:Integer} <: SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
I::Array{Ti,1} # row indices
J::Array{Ti,1} # column indices
V::Array{Tv,1} # value
Imax::Ti # number of rows
Jmax::Ti # number of columns
count::Ti # number of non-zero
maxcount::Ti # counter that triggers the allocation of more memory
end
import Base.setindex!
function setindex!(A::MySparseMatrix, x, i1::Int64, i2::Int64, I::Int64...)
if abs(x) > 0
if A.count == A.maxcount
A.maxcount += 1000000
sizehint!(A.I, A.maxcount)
sizehint!(A.J, A.maxcount)
sizehint!(A.V, A.maxcount)
end
append!(A.I, i1)
append!(A.J, i1)
append!(A.V, i1)
A.count += 1
end
end
```
I have a test that allocates 7 million entries. It spends most time in `sizehint!` and `resize!` and takes 1.6s. Do I use `sizehint!` correctly?
What I had was:
```julia
using SparseArrays
mutable struct MySparseMatrix{Tv,Ti<:Integer} <: SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
I::Array{Ti,1}
J::Array{Ti,1}
V::Array{Tv,1}
Imax::Ti
Jmax::Ti
count::Ti
maxcount::Ti
end
import Base.setindex!
function setindex!(A::MySparseMatrix, x, i1::Int64, i2::Int64, I::Int64...)
@inbounds if abs(x) > 1e-12
if A.count == A.maxcount
A.maxcount += 1000000
append!(A.I, Array{Int64}(undef, 1000000))
append!(A.J, Array{Int64}(undef, 1000000))
append!(A.V, Array{Float64}(undef, 1000000))
end
A.I[A.count] = i1
A.J[A.count] = i2
A.V[A.count] = x
A.count += 1
end
end
```
Which takes 1s and spends a lot of time in `copyto!`.
Hence, my question about more clever vector data structures.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 7, 2020, 1:36am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/5 "2020-04-07T01:36:03Z")
Note that you can already construct `SparseMatrix` instances with the `I,J,V` you dont’ need to reinvent the wheel here.
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 7, 2020, 1:44am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/6 "2020-04-07T01:44:45Z")
A litte context: I want to provide `ForwardDiff.jacobian!` with a sparse matrix instead of a dense matrix, because the Jacobian is very large but sparse.
In this [line](https://github.com/JuliaDiff/ForwardDiff.jl/blob/920717a5960433488043263b774c6028308fa8bd/src/jacobian.jl#L126) the elements of the resulting Jacobian matrix are set one-by-one.
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [April 7, 2020, 2:47am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/7 "2020-04-07T02:47:11Z")
Sounds like [GitHub - JuliaDiff/SparseDiffTools.jl: Fast jacobian computation through sparsity exploitation and matrix coloring](https://github.com/JuliaDiff/SparseDiffTools.jl) might be useful
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 7, 2020, 2:54am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/8 "2020-04-07T02:54:36Z")
True! But what I had in mind does not rely on any a priori sparsity, that is sparsity of the Jacobian that holds globally over the space of input variables.
Image a function which has a Jacobian that is a priori dense, because every element of the Jacobian has a spot in the space of input variables where it is non-zero. But ex post, when evaluating the function’s Jacobian at a certain point in the state space, it comes out extremely sparse.
Maybe the dense Jacobian does not fit into memory, but the chunks that `ForwardDiff` work on do. Or the extreme sparsity gives a performance advantage when using the Jacobian for some calculation.
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [April 7, 2020, 7:17am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/9 "2020-04-07T07:17:48Z")
> [@MaximilianJHuber](#):
>
> Do I use `sizehint!` correctly?
It seems you are trying to re-implement the already built-in array growing yourself. There is no need. Julia itself already does that, just use `append!` on its own. You can use a single `sizehint!` if you know approximately how many non-zero entries the final arrays will have but calling it like you are doing now is pointless (in fact it can be pessimising).
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 7, 2020, 2:01pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/10 "2020-04-07T14:01:01Z")
I see, but what happens when I append beyond the size hint given? Do I read you correctly that Julia grows (allocate, copy, maybe garbage collect) the vector by a large chunk when necessary so that the next thousands of `append!`s are cheap?
I do not know how many elements the vector may have.
With the vectors becoming very large I imaging that it is increasingly difficult to store them [continuously](https://stackoverflow.com/questions/34757689/is-julias-vectorvectort-stored-contiguously-in-memory) in memory. So I will experiment with `Vector{Vector{Int64}}`. Give a size hint whenever I add another inner vector and add another when it is full. Does that make sense?
---
### Author:  [@rdeits](https://discourse.julialang.org/u/rdeits)
#### Post date: [April 7, 2020, 2:04pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/11 "2020-04-07T14:04:58Z")
> [@MaximilianJHuber](#):
>
> I see, but what happens when I append beyond the size hint given? Do I read you correctly that Julia grows (allocate, copy, maybe garbage collect) the vector by a large chunk when necessary so that the next thousands of `append!` s are cheap?
Yes, except that the amount of over-allocation is proportional to the current size of the array, rather than being some fixed amount. This turns out to be important to achieve [amortized constant time](https://en.wikipedia.org/wiki/Amortized_analysis) for `push!`. This is the same behavior of pretty much any other dynamic array type (like `std::vector` in C++).
> [@MaximilianJHuber](#):
>
> Give a size hint whenever I add another inner vector and add another when it is full. Does that make sense?
This seems like unnecessary complexity. Try just using `push!` or `append!` first and see if that is sufficient.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 8, 2020, 8:42am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/12 "2020-04-08T08:42:51Z")
Hi you may have a look at Extendable Sparse.jl.
Disclaimer: I’m the author…
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 8, 2020, 3:20pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/13 "2020-04-08T15:20:07Z")
This is perfect! [ExtendableSparse.jl](https://github.com/j-fu/ExtendableSparse.jl) is about as fast as my prototype in a test where I pre-allocate vectors that fit all elements that are successively added!
When this [issue](https://github.com/j-fu/ExtendableSparse.jl/issues/4) is resolved, one should be able to provide an `ExtendableSparseMatrix` to `ForwardDiff.jacobian!`, right?
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 8, 2020, 8:28pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/14 "2020-04-08T20:28:41Z")
In the moment I doubt this. I’ll investigate this as it indeed would be intriguing if it would work.
The way I use it is in the assembly of coupled nonlinear PDEs based on edge callbacks (FEM with cell callbacks should work in a similar way). From each edge or cell I get a dense local matrix from ForwardDiff.jacobian! and sort these entries into the global matrix of type ExtendableSparse. So I don’t make any assumptions on sparsity handling in ForwardDiff . In a way this is not optimal if I have sparse couplings between the PDEs.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 10, 2020, 8:57pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/15 "2020-04-10T20:57:16Z")
Ok, with v0.3.0 ExtendableSparse cooperates well with ForwardDiff.jl . Thanks for asking - I never found the time to try this out before.
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 11, 2020, 3:04am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/16 "2020-04-11T03:04:08Z")
I can confirm that a `ExtendableSparseMatrix` works as Jacobian argument in a `ForwardDiff.jacobian!` call. However, the performance is much worse than when using a `SparseMatrix`:
```julia
function f!(dx, x)
dx .= x[1]
end
using ForwardDiff, SparseArrays, ExtendableSparse, BenchmarkTools
n = 100
m = 100000
jac_ext = ExtendableSparseMatrix(m, n);
@btime ForwardDiff.jacobian!(jac_ext, f!, zeros(m), zeros(n));
jac = spzeros(m, n);
@btime ForwardDiff.jacobian!(jac, f!, zeros(m), zeros(n));
all(jac .== jac_ext)
```
Using the code from `ExtendableSparse`s [README](https://github.com/j-fu/ExtendableSparse.jl) yields the same results.
This is despite a `ExtendableSparseMatrix` being much faster than a `SparseMatrix` in a test where I successively add elements, like ForwardDiff [does](https://github.com/JuliaDiff/ForwardDiff.jl/blob/52f4be7d5a8d858edb9f8f08f9574b97a08f46f1/src/jacobian.jl#L126).
@j-fu have you done performance tests?
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 11, 2020, 9:29am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/17 "2020-04-11T09:29:40Z")
Well, not for this kind of use with ForwardDiff. In fact, the performance very much depends on the way things are called internally in ForwardDiff which I did not investigate.
Your function is kind of the worst case for the linked list internal format used in ExtendableSparse, as it generates a full row in the matrix (will have to mention this caveat in the readme: I assume that “sparse” means that all rows/columns have \<\< n entries). But it seems that this is not the problem here.
However, this brings me back to my previous post: “atomic” assembly of the finite difference operator without too much assumptions on the inner workings of ForwardDiff. So let me extend your example a bit (VoronoiFVM is my package):
```julia
using ExtendableSparse
using ForwardDiff
using LinearAlgebra
using SparseArrays
using VoronoiFVM
#
# Finite difference operator for
# -\Delta u^2 + u^2 = 1
#
function ffd(y,x)
n=length(x)
h=1/(n-1)
y[1]=(x[1]^2-1)*0.5*h
for i=2:n-1
y[i]=(x[i]^2-1.0)*h
end
y[end]=(x[end]^2-1)*0.5*h
for i=1:n-1
dx=(x[i+1]^2-x[i]^2)/h
y[i+1]+=dx
y[i]-=dx
end
end
function flux!(f,u,edge)
f[1]=u[1,1]^2-u[1,2]^2
end
## Storage term
function storage!(f,u,node)
f[1]=u[1]
end
function reaction!(f,u,node)
f[1]=u[1]^2-1.0
end
function runtest(n)
jac_ext = ExtendableSparseMatrix(n, n);
@time begin
ForwardDiff.jacobian!(jac_ext, ffd, ones(n), ones(n));
flush!(jac_ext)
end
jac = spzeros(n, n);
@time ForwardDiff.jacobian!(jac, ffd, ones(n), ones(n));
all(jac .== jac_ext)
# Set up stuff for VoronoiFVM which uses atomic assembly
h=1.0/convert(Float64,n-1)
X=collect(0:h:1)
grid=VoronoiFVM.Grid(X)
physics=VoronoiFVM.Physics(flux=flux!, storage=storage!,reaction=reaction!)
sys=VoronoiFVM.System(grid,physics,unknown_storage=:dense)
enable_species!(sys,1,[1])
solution=unknowns(sys,inival=1.0)
oldsol=unknowns(sys,inival=1.0)
residual=unknowns(sys,inival=1.0)
tstep=1.0e100 # large timestep aproximates stationary problem
@time begin
VoronoiFVM.eval_and_assemble(sys,solution,oldsol,residual,tstep)
flush!(sys.matrix)
end
all(sys.matrix.≈jac_ext)
end
```
second run:
```julia
julia> runtest(20000)
5.806884 seconds (32 allocations: 7.174 MiB)
3.925591 seconds (41 allocations: 6.275 MiB)
0.025753 seconds (117.03 k allocations: 3.849 MiB)
true
```
I intentionally use @time here, as @btime probably won’t be dominated by
the initial structure buildup phase. This is the kind of things I also benchmarked
before (see the benchmarking stuff with fdrand!() in ExtendableSparse).
In VoronoiFVM, I call ForwardDiff with flux!, reaction!, storage! and get local matrices which I then assemble into the global matrix.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 11, 2020, 10:27am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/18 "2020-04-11T10:27:50Z")
Just an add-on here: here is a timing plot:
Seemingly, ForwardDiff seems to consider all possible matrix entries due to being O(n^2), and I guess this can be improved only with a priori information on sparsity.
OTOH, “atomic” assembly in VoronoiFVM automatically takes the sparsity into account in an implicit manner.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [April 11, 2020, 10:43am UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/19 "2020-04-11T10:43:18Z")
The explantion of the worse performance of ExtenableSparse with ForwardDiff ist then the following:
```julia
function Base.setindex!(ext::ExtendableSparseMatrix{Tv,Ti}, v, i,j) where{Tv,Ti<:Integer}
k=findindex(ext.cscmatrix,i,j)
if k>0
ext.cscmatrix.nzval[k]=v
else
if ext.lnkmatrix==nothing
ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
end
ext.lnkmatrix[i,j]=v
end
end
```
I assume the matrix is always a sum of a cscmatrix and a lnkmatrix. This allows to update entries and flush! multiple times. If we try to insert a zero into a non-existing slot, we have to look up both, so searching twice. This IMHO explains the performance difference between ExtendableSparse and SparseMatrixCSC if we try to insert mostly zero values.
Inserting a zero into SparseMatrixCSC is connected with a search as well, see
[sparsematrix.jl](https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/sparsematrix.jl)
I guess if this case would be catched in ForwardDiff, we would see a different picture. As far as I can see in [jacobian.jl](https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/src/jacobian.jl) this is not done.
---
### Author:  [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)
#### Post date: [April 11, 2020, 1:26pm UTC](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133/20 "2020-04-11T13:26:09Z")
Thank you, for that detailed explanation!
I think the only way forward is to make [adjustments](https://github.com/JuliaDiff/ForwardDiff.jl/issues/452) in `ForwardDiff`. And it seems an to be an easy fix!
[Next page](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133.md?page=2)