# Taking advantage of rowwise-constant sparse matrix

**URL:** <https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366>\
**Category:** Performance\
**Tags:** linearalgebra, sparse\
**Created:** [August 30, 2021, 3:53pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366 "2021-08-30T15:53:57Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [August 30, 2021, 3:53pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/1 "2021-08-30T15:53:57Z")

</div>

Hi,

This is a follow-up of [this previous thread](https://discourse.julialang.org/t/faster-matrix-product-with-a-specifically-structured-matrix/67355/19) where we discussed with @Sukera the possibility to improve matrix-vector product in the case where the matrix is sparse, but also has only one possible non-zero value per line.

Mathematically, this predicate on the matrix structure should allow a factorization in the multiplication algorithm and hence enhance performance. I managed to code a simple MWE :

```julia
using SparseArrays
n,m = 1000,1000
x = rand(m)
P = sparse((rand(n,m).>0.97) .* rand(m));

# check predicate and decompose: 
function decompose(P)
    @assert all([sum(unique(P[i,:]).!=0)<=1 for i in 1:size(P,1)])
    b = [sum(unique(P[i,:])) for i in 1:size(P,1)]
    M = sparse(BitMatrix(P .!= 0))
    return b,M
end

b,M = decompose(P)
@inline my_prod(b,M,x) = b.*(M*x)

function my_prod(r,b,M,x)
    @inbounds for i in size(M,1)
        r[i] = b[i]*sum(x[M[i,:]])
    end
    return r
end
r = P*x
# check : 
my_prod(r,b,M,x) ≈ P*x

# test : 
using BenchmarkTools
@btime P*$(x);
@btime my_prod(r,b,M,$x);

```

Which gives :

```julia
julia> @btime P*$(x);
  59.700 μs (1 allocation: 7.94 KiB)

julia> @btime my_prod(r,b,M,$x);
  21.500 μs (11 allocations: 1.16 KiB)

julia> 

```

So a 1/2 time. I was hoping for more, but I start wandering if my testing is fair…

Did I make something wrong ? I notice that my code is allocating a lot, could this be fixed ?

In my application, the matrix `P` is held fixed and the vector `x` changes from one iteration to the other, hence my concerns about whether this particular structure could be exploited a little more.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [August 30, 2021, 4:23pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/2 "2021-08-30T16:23:33Z")

</div>

Apologies, I was cooking my dinner. That’s why I didn’t answer right away in the other thread.

> [@lrnv](#):
>
> ```julia
> function my_prod(r,b,M,x)
> @inbounds for i in size(M,1)
> r[i] = b[i]*sum(x[M[i,:]])
> end
> return r
> end
> 
> ```

This function doesn’t do what you think it does.

- `size(M,1)` is a single integer - the loop is only running for a single iteration because numbers are iterable in julia and thus the result must be wrong. You probably meant to write `axes(M, 1)` or `1:size(M,1)`. The former gives an iterable over the first axis, the latter constructs the same manually. Your check still passes because _only the last entry is overwritten_ and `r` still holds the correct result from `P*x`.

```julia
julia> r_ = zeros(size(r)); 
                            
julia> my_prod(r_,b,M,x) ≈ r
false                       
                            
julia> findall(!=(0), r_)   
1-element Vector{Int64}:    
 1000                       

```

- Since `M[i,:]` is a vector, `x[M[i,:]]` will slice (and thus allocate). This doesn’t hurt too bad for one iteration, but it will _kill_ performance for all appropriate iterations:

```julia
julia> function my_prod(r,b,M,x)         
           for i in 1:size(M,1) # fixed iterations         
               r[i] = b[i]*sum(x[M[i,:]])
           end                           
           return r                      
       end                               

my_prod (generic function with 1 method) 
julia> @btime P*x;                        
  29.900 μs (1 allocation: 7.94 KiB)      
                                          
julia> @btime my_prod(r,b,M,$x);          
  39.710 ms (7038 allocations: 952.67 KiB)

julia> r_ = zeros(size(r)); 
                            
julia> my_prod(r_,b,M,x) ≈ r
true                        

```

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [August 30, 2021, 4:28pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/3 "2021-08-30T16:28:29Z")

</div>

Indeed, a silly mistake… And hence the result is 1000x slower… Thanks for catching it.

The inline version still has the same performance as the direct one, with 2 allocations instead of one.

let’s forget about the unrolled version, wich I cannot make a lot faster than 500x times slower than the direct approach. What about the inline version :

```julia
using SparseArrays
n,m = 1000,1000
x = rand(m)
P = sparse((rand(n,m).>0.97) .* rand(m));

# check predicate and decompose: 
function decompose(P)
    @assert all([sum(unique(P[i,:]).!=0)<=1 for i in 1:size(P,1)])
    b = [sum(unique(P[i,:])) for i in 1:size(P,1)]
    M = sparse(BitMatrix(P .!= 0))
    return b,M
end

b,M = decompose(P)
@inline my_prod(b,M,x) = @inbounds b.*(M*x)
my_prod(b,M,x) ≈ P*x

# test : 
using BenchmarkTools
@btime P*$(x);
@btime my_prod(b,M,$x);

```

Gives:

```julia
julia> @btime P*$(x);
  59.500 μs (1 allocation: 7.94 KiB)

julia> @btime my_prod(b,M,$x);
  60.600 μs (2 allocations: 15.88 KiB)

```

Maybe the `SparseMatrixCSC{Bool, Int64}` product already takes adventage of the `Bool` eltype ?

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [August 30, 2021, 4:37pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/4 "2021-08-30T16:37:30Z")

</div>

I think you really need to exploit the sparsity of this `M`, or of the original `P`. Here’s a fast version of the idea to slice `M` and multiply by the dense `x`, and it’s still many times slower than the original:

```julia
julia> function my_prod_dot2(r,b,MT,x)
           @inbounds for i in axes(MT,2)
               r[i] = b[i] * dot(x, @view MT[:,i])
           end
           return r
       end
my_prod_dot2 (generic function with 1 method)

julia> MT = Matrix{Float64}(M');

julia> @btime my_prod_dot2($r,$b,$MT,$x);
  406.125 μs (0 allocations: 0 bytes)

julia> @btime $P*$(x);
  18.792 μs (1 allocation: 7.94 KiB)

julia> function my_prod_dense(r,b,MT,x)
         mul!(r, MT', x)
         r .*= b
         r
       end;

julia> @btime my_prod_dense($r,$b,$MT,$x);
  49.333 μs (0 allocations: 0 bytes)

julia> r ≈ P * x
true

```

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [August 30, 2021, 4:38pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/5 "2021-08-30T16:38:34Z")

</div>

> [@lrnv](#):
>
> The inline version still has the same performance as the direct one, with 2 allocations instead of one.

> [@lrnv](#):
>
> ` M = sparse(BitMatrix(P .!= 0))`

That’s because you’re doing a sparse multiplication again:

> [@lrnv](#):
>
> `b.*(M*x)`

The sparse matmul `M*x` allocates its result `res`, and `b .* res` allocates its result. You’re not going to beat sparse matrix multiplication by using sparse matrix multiplication.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [August 30, 2021, 4:39pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/6 "2021-08-30T16:39:37Z")

</div>

> [@Sukera](#):
>
> You’re not going to beat sparse matrix multiplication by using sparse matrix multiplication.

Indeed 😉

I guess my next move is to go see how this sparse matrix multiplication is implemented.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [August 30, 2021, 4:44pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/7 "2021-08-30T16:44:10Z")

</div>

No, there is no smart trick to apply here. If you don’t know a lot more about the specific structure in your `P`, you basically still have random entries in your matrix to access. Whether or not the values in a column are set to the same value is not going to matter all that much compared to checking whether an index should be multiplied at all.

The only way I can think of that _may_ be faster than this is by only saving a single element per column. This _may_ lead to always-hot caches, but I doubt it will be enough to meaningfully impact the runtime because I’m fairly certain that the index lookups dominate here.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [August 30, 2021, 4:48pm UTC](https://discourse.julialang.org/t/taking-advantage-of-rowwise-constant-sparse-matrix/67366/8 "2021-08-30T16:48:58Z")

</div>

Thanks for taking the time anyway. Two things:

1° My _rows_ are constants, not my _columns_.  
2° Indeed, the lookup might be the problem more than the floats multiplications. With more involved multiplications, this pays out:

```julia
using SparseArrays
n,m = 1000,1000
ElType = BigFloat
x = ElType.(rand(m))
P = sparse(ElType.((rand(n,m).>0.97) .* rand(m)));

# check predicate and decompose: 
function decompose(P)
    @assert all([sum(unique(P[i,:]).!=0)<=1 for i in 1:size(P,1)])
    b = [sum(unique(P[i,:])) for i in 1:size(P,1)]
    M = sparse(BitMatrix(P .!= 0))
    return b,M
end

b,M = decompose(P)
@inline my_prod(b,M,x) = @inbounds b.*(M*x)
my_prod(b,M,x) ≈ P*x

# test : 
using BenchmarkTools
@btime P*$(x);
@btime my_prod(b,M,$x);

```

gives:

```julia
julia> @btime P*$(x);
  9.727 ms (126113 allocations: 6.67 MiB)

julia> @btime my_prod(b,M,$x);
  6.875 ms (248222 allocations: 11.14 MiB)

```

and for `DoubleFloats.Double64` numbers this gives:

```julia
julia> @btime P*$(x);
  276.300 μs (1 allocation: 15.75 KiB)

julia> @btime my_prod(b,M,$x);
  154.400 μs (2 allocations: 31.50 KiB)

```

For Float64, the indexing overhead is too high and the multiplication cost too low for this approach to be worth it. But still, I’m not just using SparseArrays multiplications to replace SparseArrays multiplications 😉