# Weighted sum of sparse matrices

**URL:** <https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564>\
**Category:** New to Julia\
**Created:** [August 25, 2017, 6:48pm UTC](https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564 "2017-08-25T18:48:38Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [August 25, 2017, 6:48pm UTC](https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564/1 "2017-08-25T18:48:38Z")

</div>

I have a simple problem:

```julia
julia> const DIM = 1000;
julia> const ORDER = 500;
julia> const matrix_reps = [sprand(DIM, DIM, 0.01) for _ in 1:ORDER;
julia> const weights = sprandn(ORDER, 0.1);

```

I want to compute weighted sum of `matrix_reps`, i.e. \sum\_i weights[i]\cdot matrix\_reps[i].  
This is what I tried so far:

```julia
julia> function w_sumFUSED(x::AbstractVector{T}, mreps::Vector{TT}) where {S<:Number, T<:Number, TT<:AbstractArray{S}}
           res = spzeros(S, size(first(mreps))...)
           for i in findn(x)
               res .+= S(x[i]).*mreps[i]
           end
           return res
       end
w_sumFUSED (generic function with 1 method)

julia> function w_sumBLAS(x::AbstractVector{T}, mreps::Vector{TT}) where {S<:Number, T<:Number, TT<:AbstractArray{S}}
           res = spzeros(S, size(first(mreps))...)
           for i in findn(x)
               Base.LinAlg.axpy!(S(x[i]), mreps[i], res)
           end
           return res
       end
w_sumBLAS (generic function with 1 method)

julia> function w_sumSUM(x::AbstractVector{T}, mreps::Vector{TT}) where {S<:Number, T<:Number, TT<:AbstractArray{S}}
           return sum(S(x[i])*mreps[i] for i in findn(x))
       end
w_sumSUM (generic function with 1 method)

```

timings:

```julia
julia> function time_everything(vect, mreps)
           w_sumSUM(vect[1:10], mreps[1:10]);
           println("SUM:")
           @time resSUM = w_sumSUM(vect, mreps);

           w_sumBLAS(vect[1:10], mreps[1:10])
           println("BLAS:")
           @time resBLAS = w_sumBLAS(vect, mreps);

           w_sumFUSED(vect[1:10], mreps[1:10])
           println("fused:")
           @time resFUSED = w_sumFUSED(vect, mreps);

           @assert resFUSED == resBLAS
           @assert resFUSED == resSUM
       end
time_everything (generic function with 1 method)

julia> time_everything(weights, matrix_reps)
SUM:
  0.145636 seconds (688 allocations: 151.844 MiB, 47.66% gc time)
BLAS:
 14.278971 seconds (41 allocations: 10.010 MiB)
fused:
 15.503572 seconds (41 allocations: 10.010 MiB)

```

Question(s):

- why is BLAS so slow, despite so little allocations? (checked against `@code_typewarn`)
- why is FUSED so slow, despite so little allocations? (checked against `@code_typewarn`)
- (general) what is the best method for performing arithmetic on sparse matices?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 25, 2017, 8:39pm UTC](https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564/2 "2017-08-25T20:39:05Z")

</div>

`axpy!` is not calling BLAS here, because BLAS doesn’t operate on sparse matrices.

Summing two sparse matrices is going to be slow if the sparsity patterns don’t match, because it is going to involve reallocating the sparse-matrix data structure and moving a lot of data around. Google how sparse matrices are stored (Julia uses compressed sparse column, but compressed sparse row is similar).

It will be faster if you precompute the sparsity pattern of the result, which is what I think `sum` may be doing?

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [August 26, 2017, 10:03am UTC](https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564/3 "2017-08-26T10:03:44Z")

</div>

Thanks for the answer,  
I didn’t expect inserting new values is so expensive, As far as I understand summing sparse matrices (when sparsity patterns don’t match) involves insertion into `nzval`s, `rowvals` and incrementing `colptr`.

Also: after changing `res = spzeros(S, size(first(mreps))...)` to dense `res = zeros(S, size(first(mreps))...)` (and thus avoiding the problem You pointed out?)  
`BLAS` and `FUSED` are faster, but still not as fast as `SUM`:

```julia
SUM:
  0.083839 seconds (943 allocations: 265.821 MiB, 6.02% gc time)
BLAS:
  0.791856 seconds (3 allocations: 7.630 MiB)
fused:
  0.759364 seconds (3 allocations: 7.630 MiB, 0.04% gc time)

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 27, 2017, 11:35am UTC](https://discourse.julialang.org/t/weighted-sum-of-sparse-matrices/5564/4 "2017-08-27T11:35:10Z")

</div>

> [@abulak](#):
>
> I didn’t expect inserting new values is so expensive, As far as I understand summing sparse matrices (when sparsity patterns don’t match) involves insertion into nzvals, rowvals and incrementing colptr.

How do you think you insert a value into the middle of an array without reallocating and moving everything in the array?

I don’t think broadcast operations on mixed combinations of sparse and dense arrays are particularly fast. It might be interesting to try a variant where you allocate a dense output array, but use explicit loops over the nonzero entries to do the sums.