# Getindex A\[i,j\] wrong complexity on SparseMatrixCSC?

**URL:** <https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277>\
**Category:** Numerics\
**Tags:** question, performance, linearalgebra, sparse\
**Created:** [April 14, 2021, 1:52pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277 "2021-04-14T13:52:14Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 1:52pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/1 "2021-04-14T13:52:14Z")

</div>

The complexity of extracting a sparse submatrix `A[i,j]` of size k x l from a large sparse matrix of size n with nnz entries should be O(k log k + l nz log k).

Accessing my sparse matrices seemed very slow and I checked the algorithms… They seem fairly optimized but the complexity does not match the one I stated above. I ran some benchmarks, keeping the number nonzero entries per column constant, as well as the k x l submatrix, while scaling the overall problem size. It seems that the problem scales linearly with the problem size, which should be avoided at all cost!

 ![Screenshot 2021-04-14 at 16.29.01](https://global.discourse-cdn.com/julialang/original/3X/e/6/e644f6f8846cdfa4b6725f0096e99c1a32e4891a.png)  
 ![Screenshot 2021-04-14 at 16.33.15](https://global.discourse-cdn.com/julialang/original/3X/1/6/169fad3a915678c1ac61df664d6b6d813fb43e96.png)

Is this a misunderstanding on my side/did I miss something? Could someone explain to me the dependency on the total size of the matrix? In my opinion this shouldn’t be happening with CSC matrices.

I have uploaded my tests online, you can run them as Pluto notebook or just as a script:

> <https://github.com/bonevbs/nla-notebooks/blob/main/sparse_getindex.jl>

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 14, 2021, 2:16pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/2 "2021-04-14T14:16:30Z")

</div>

Isn’t this the expected complexity? You said you expected `O(K*log(k) + L*nz*log(K))` which is more than linear in `K` and `L`.

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 2:19pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/3 "2021-04-14T14:19:57Z")

</div>

k and l are the size of I and J and nz is the number of nonzero entries per row. The plot shows the timing if the overall size of A is increased while keeping k, l, nz constant. So no, this is not what I would expect

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 14, 2021, 2:23pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/4 "2021-04-14T14:23:36Z")

</div>

Can you upload a version of the plot with labeled axes?

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 2:24pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/5 "2021-04-14T14:24:37Z")

</div>

good point, let me correct that

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 14, 2021, 2:29pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/6 "2021-04-14T14:29:13Z")

</div>

Here’s the relevant section of the notebook, since it wasn’t clear to me what you were measuring by your description alone:

> <https://github.com/bonevbs/nla-notebooks/blob/9aceb046281e978f5f6fb8df230a04d22e78d3be/sparse_getindex.jl#L36-L51>

In short, you’re looking specifically at `getindex(A::SparseMatrixCSC, I::Vector{Int}, J::Vector{Int})`, with `I` sometimes sorted (and contiguous).

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 2:31pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/7 "2021-04-14T14:31:50Z")

</div>

yes, I am trying to verify the complexity of getindex being quasilinear w.r.t. k, l, nz. In my understanding there should be no scaling with the problem size.

I can pre-sort `I` if necessary but that doesn’t seem to be the issue here. (I have also updated the plot)

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 14, 2021, 2:41pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/8 "2021-04-14T14:41:05Z")

</div>

Some of these indexing implementations use a cache the length of the entire column (and not just the subset), which does indeed seem suboptimal. That’s likely where this is coming from. Finding a way to limit that to the size of the index would be great.

> <https://github.com/JuliaLang/julia/blob/2fed7c34739cfe97610fa648b8ddc29b3c5141d6/stdlib/SparseArrays/src/sparsematrix.jl#L2318>

> <https://github.com/JuliaLang/julia/blob/2fed7c34739cfe97610fa648b8ddc29b3c5141d6/stdlib/SparseArrays/src/sparsematrix.jl#L2379-L2383>

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 2:51pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/9 "2021-04-14T14:51:39Z")

</div>

Oh wow - thank you for finding this. I was questioning my sanity.

I had completely overlooked this. This seems quite problematic to me. The whole point of sparse matrices is to avoid any scaling related to the overall problem-size. The stated complexity is for `getindex_I_sorted_bsearch_I`… Using a cache like this seems to undermine the entire implementation of the algorithm…

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 14, 2021, 2:55pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/10 "2021-04-14T14:55:00Z")

</div>

Yup, definitely should be improved. Want to take a crack at a pull request? I’m happy to help along the way if it’s new to you.

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 3:00pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/11 "2021-04-14T15:00:20Z")

</div>

Definitely, glad to accept your help as well as I have never done this before. I will have a crack at it!

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 14, 2021, 3:45pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/12 "2021-04-14T15:45:16Z")

</div>

Just to get started, a really helpful interactive development/debugging technique is to copy the existing implementation into your favorite IDE/REPL, edit it, and then update the definition with an `@eval SparseArrays` to re-evaluate it in the context of the SparseArrays module. You may also want to remove `@inbounds` while developing to make errors more obvious.

Once you have an implementation that makes you happy, you can just edit it in the browser on GitHub directly on the page I linked above.

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 14, 2021, 3:46pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/13 "2021-04-14T15:46:57Z")

</div>

Thank you so much - this is very helpful! I had copied the function into the above notebook but I think I will follow your advice 🙂

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 17, 2021, 9:25pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/14 "2021-04-17T21:25:26Z")

</div>

Hi, I had some time to have a go at it and I wrote new routines for `getindex_I_sorted_linear` and `getindex_I_sorted_bsearch_I`. As pointed out above, they had suboptimal complexities with dependency on the overall problem size, which can be quite catastrophic, when one does scaling studies for instance. As a consequence, it seems to me that my new routines are not quite as optimised (I have little experience in optimising Julia code using `@simd` and `@inbounds`), so I would be grateful if someone could have look at it.

That being said, the complexities are now correct, as can be verified in the repository that I have linked above. I have also attached some plots that prove this.

 ![Screenshot 2021-04-17 at 23.22.49](https://global.discourse-cdn.com/julialang/original/3X/4/a/4a722cfbea9c18c08a91b30560e8e12fc6ea19cb.png)

The code can be found here: [https://github.com/bonevbs/nla-notebooks/blob/main/mygetindex.jl](https://github.com/bonevbs/nla-notebooks/blob/main/mygetindex.jl)

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 17, 2021, 9:34pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/15 "2021-04-17T21:34:41Z")

</div>

I should note that `getindex_I_sorted_bsearch_I` is the function which I would like tips on how to optimise it as I feel the algorithm is the correct one and satisfactory to be put into a pull request.

For `getindex_I_sorted_linear`, I am currently using a `Dict{Int,Int}`, which effectively implements a hashmap. This seems to be quite slow though, even though at some point it would win. For index sets that are relatively compressed one could stick with the current implementation and just reduce the cache to the envelope of the index set. I am not entirely sure whether I should write a custom hashmap algorithm for this one…

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 18, 2021, 5:22am UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/16 "2021-04-18T05:22:27Z")

</div>

I think you should probably make the pull request now. There is probably a bunch of improvements, but making a PR is a pretty decent way to get eyes on it.

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 18, 2021, 7:51am UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/17 "2021-04-18T07:51:24Z")

</div>

ok, thanks for the input 🙂 I made the pull-request here [Change getindex algorithms on sparse matrices by bonevbs · Pull Request #40519 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/40519)

---

<div class="post-metadata">

**Author:** ![viralbshah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/viralbshah/32/54_2.png) [@viralbshah](https://discourse.julialang.org/u/viralbshah)\
**Post date:** [April 19, 2021, 8:48pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/18 "2021-04-19T20:48:23Z")

</div>

Can you also try very large n (\> 10^6) and \< 1 nonzero/row?

-viral

---

<div class="post-metadata">

**Author:** ![bonevbs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bonevbs/32/24029_2.png) [@bonevbs](https://discourse.julialang.org/u/bonevbs)\
**Post date:** [April 19, 2021, 11:30pm UTC](https://discourse.julialang.org/t/getindex-a-i-j-wrong-complexity-on-sparsematrixcsc/59277/19 "2021-04-19T23:30:12Z")

</div>

this should be covered by `getindex\_I\_sorted\_bsearch\_A, right? I will have some time on the weekend to look at it