# Matrix (dense) vector (sparse) product improvement

**URL:** <https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [June 23, 2019, 6:01pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579 "2019-06-23T18:01:55Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [June 23, 2019, 6:01pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/1 "2019-06-23T18:01:55Z")

</div>

Products `W'*x_sp ` seem very slow when `W` is a dense array and `x_sp` is an sparse vector. Odly enough this does not happen if we don’t have and adjoint (`W*x_sp ` is fast but `W'*x_sp ` is very slow).

Is there an easy way to overload `*` operator in such a way that `*(W',x)` uses my own code instead of the default one? I’m not sure how to do it because `typeof(W')` is still an Array.

Would it be sensible to add a function to cover this in base so that performance is improved?

Minimal working example:

#### Case with adjoint

```julia

using SparseArrays
using BenchmarkTools

function mat_dense_vec_sparse(W_dense, x_sp)
    n_rows_W, n_cols_W = size(W_dense)
    result = zeros(eltype(W_dense), n_rows_W)
    
    @inbounds for j in 1:n_rows_W
        for i in x_sp.nzind
            result[j] += W_dense[j,i] * x_sp[i] 
        end
    end
    return result
end

percentage_sparse_tests = [0.5, 0.4, 0.3, 0.2, 0.1, 0.01, 0.001,0.0001]
results = []

for percent_sparse in percentage_sparse_tests
    W = rand(n_features, n_classes);
    x_sp = sprand(n_features, percent_sparse);
    b = zeros(n_classes);
    
    time1 = @benchmark W'*x_sp;
    time2 = @benchmark mat_dense_vec_sparse(W', x_sp);

    result1 = W'*x_sp + b
    result2 = mat_dense_vec_sparse(W', x_sp)

    time1_meantime = Int(round(mean(time1.times)))
    time2_meantime = Int(round(mean(time2.times)))
    improvement = round(time1_meantime/ time2_meantime, digits = 2)
    percentage_nonzeros = round(100*(length(nonzeros(x_sp))/n_features), digits=2)
    
    x = (percentage_sparse=percentage_nonzeros, t1=time1_meantime, t2=time2_meantime, improvement=improvement)
    push!(results, x)
    
    println("\nTrue % nonzeros:", percentage_nonzeros, "\t percent_sparse given:", 100*percent_sparse)

    print("Improvement: ", improvement,"x ")
    print("\tW'*x_sp : ", time1_meantime)
    print("\tCustom:", time2_meantime)
    print("\tisapprox: ", isapprox(result1, result2))
    println(" ")

end

```

```julia-auto
True % nonzeros:49.7 percent_sparse given:50
Improvement: 1.85x W'*x_sp : 847733	Custom:458838	isapprox: true 

True % nonzeros:40.17 percent_sparse given:40
Improvement: 2.03x W'*x_sp : 746418	Custom:368248	isapprox: true 

True % nonzeros:29.97 percent_sparse given:30
Improvement: 2.54x W'*x_sp : 690894	Custom:272434	isapprox: true 

True % nonzeros:20.43 percent_sparse given:20
Improvement: 2.91x W'*x_sp : 554531	Custom:190674	isapprox: true 

True % nonzeros:10.45 percent_sparse given:10
Improvement: 4.13x W'*x_sp : 403204	Custom:97743	isapprox: true 

True % nonzeros:0.99 percent_sparse given:1
Improvement: 105.13x W'*x_sp : 183549	Custom:1746	isapprox: true 

True % nonzeros:0.13 percent_sparse given:0.1
Improvement: 452.83x W'*x_sp : 108227	Custom:239	isapprox: true 

True % nonzeros:0.01 percent_sparse given:0.01
Improvement: 552.5x W'*x_sp : 41990	Custom:76	isapprox: true 

```

#### Case without adjoint

```julia
percentage_sparse_tests = [0.5, 0.4, 0.3, 0.2, 0.1, 0.01, 0.001,0.0001]
results = []
W_t = copy(W')

for percent_sparse in percentage_sparse_tests
    W = rand(n_features, n_classes);
    x_sp = sprand(n_features, percent_sparse);
    b = zeros(n_classes);
    
    time1 = @benchmark W_t*x_sp;
    time2 = @benchmark mat_dense_vec_sparse(W_t, x_sp);

    result1 = W_t*x_sp + b
    result2 = mat_dense_vec_sparse(W_t, x_sp)

    time1_meantime = Int(round(mean(time1.times)))
    time2_meantime = Int(round(mean(time2.times)))
    improvement = round(time1_meantime/ time2_meantime, digits = 2)
    percentage_nonzeros = round(100*(length(nonzeros(x_sp))/n_features), digits=2)
    
    x = (percentage_sparse=percentage_nonzeros, t1=time1_meantime, t2=time2_meantime, improvement=improvement)
    push!(results, x)
    
    println("\nTrue % nonzeros:", percentage_nonzeros, "\t percent_sparse given:", 100*percent_sparse)

    print("Improvement: ", improvement,"x ")
    print("\tW'*x_sp : ", time1_meantime)
    print("\tCustom:", time2_meantime)
    print("\tisapprox: ", isapprox(result1, result2))
    println(" ")

end

```

```julia-auto
True % nonzeros:49.43 percent_sparse given:50
Improvement: 0.03x W'*x_sp : 15724	Custom:461304	isapprox: true 

True % nonzeros:40.27 percent_sparse given:40
Improvement: 0.03x W'*x_sp : 12873	Custom:372671	isapprox: true 

True % nonzeros:29.87 percent_sparse given:30
Improvement: 0.03x W'*x_sp : 9583	Custom:282339	isapprox: true 

True % nonzeros:20.22 percent_sparse given:20
Improvement: 0.04x W'*x_sp : 6634	Custom:186716	isapprox: true 

True % nonzeros:9.77 percent_sparse given:10
Improvement: 0.04x W'*x_sp : 3289	Custom:92114	isapprox: true 

True % nonzeros:1.07 percent_sparse given:1
Improvement: 0.19x W'*x_sp : 399	Custom:2134	isapprox: true 

True % nonzeros:0.08 percent_sparse given:0.1
Improvement: 0.54x W'*x_sp : 77	Custom:143	isapprox: true 

True % nonzeros:0.0 percent_sparse given:0.01
Improvement: 1.1x W'*x_sp : 55	Custom:50	isapprox: true 

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [June 23, 2019, 6:10pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/2 "2019-06-23T18:10:05Z")

</div>

Isn’t the transpose of a dense matrix an adjoint?

```julia
julia> typeof(A')
Adjoint{Float64,Array{Float64,2}}

```

That should help you write the method.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [June 23, 2019, 6:14pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/3 "2019-06-23T18:14:16Z")

</div>

I tried the following:

```julia
import Base 
import Base.*

*(W::Adjoint, x::SparseVector) = mat_dense_vec_sparse(W_dense, x_sp)

```

But I get…

```julia

MethodError: *(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64}) is ambiguous. Candidates:
  *(W::Adjoint, x::SparseVector) in Main at In[276]:4
  *(adjA::Adjoint{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:103
Possible fix, define
  *(::Adjoint{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, ::SparseVector{S,Ti} where Ti<:Integer)

```

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [June 23, 2019, 7:20pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/4 "2019-06-23T19:20:16Z")

</div>

If I try to follow the “Possible fix” that is returned in `MethodError` It seems my function is not chosen. In any case the input types seem to match the definition.

```julia

function mat_dense_vec_sparse(W_dense, x_sp)
    n_rows_W, n_cols_W = size(W_dense)
    result = zeros(eltype(W_dense), n_rows_W)
    
    @inbounds for j in 1:n_rows_W
        for i in x_sp.nzind
            result[j] += W_dense[j,i] * x_sp[i] 
        end
    end
    return result
end

my_prod(W::Adjoint{T, AbstractArray{T,2}}, x::SparseVector{T,I}) where {T<:AbstractFloat, I<:Int} = mat_dense_vec_sparse(W_dense, x_sp) 

W = rand(n_features, n_classes);
x_sp = sprand(n_features, 0.1);
b = zeros(n_classes);

typeof(W'), typeof(x_sp)

my_prod(W',x_sp)

```

```julia
typeof(W'), typeof(x_sp)
(Adjoint{Float64,Array{Float64,2}}, SparseVector{Float64,Int64})

```

```julia
my_prod(W',x_sp)
MethodError: no method matching my_prod(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64})
Closest candidates are:
  my_prod(!Matched::Adjoint{T<:AbstractFloat,AbstractArray{T<:AbstractFloat,2}}, ::SparseVector{T<:AbstractFloat,I<:Int64}) where {T<:AbstractFloat, I<:Int64} at In[29]:1
  my_prod(!Matched::Adjoint{T<:AbstractFloat,AbstractArray{T<:AbstractFloat,2}}, ::SparseVector{T<:AbstractFloat,Ti} where Ti<:Integer) where {T<:AbstractFloat, I<:Int64} at In[33]:1

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [June 24, 2019, 4:47am UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/5 "2019-06-24T04:47:54Z")

</div>

I would try to do what it suggests.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [June 24, 2019, 10:25am UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/6 "2019-06-24T10:25:10Z")

</div>

Well isn’t the following definition what was suggested?

```julia
my_prod(W::Adjoint{T, AbstractArray{T,2}}, x::SparseVector{T,I}) where {T<:AbstractFloat, I<:Int} =
     mat_dense_vec_sparse(W_dense, x_sp) 

```

This seems to match the types of the elements I am trying to call the function with

```julia
typeof(W'), typeof(x_sp)
(Adjoint{Float64,Array{Float64,2}}, SparseVector{Float64,Int64})

```

In any case now I don’t get an error when I do `*(W', x_sp)` but my method is not called. When I try to do `my_prod(W',x_sp)` I still get

```julia
MethodError: no method matching my_prod(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64})
Closest candidates are:
  my_prod(!Matched::Adjoint{T<:AbstractFloat,AbstractArray{T<:AbstractFloat,2}}, ::SparseVector{T<:AbstractFloat,I<:Int64}) where {T<:AbstractFloat, I<:Int64} at In[7]:1

```

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [June 24, 2019, 10:33am UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/7 "2019-06-24T10:33:29Z")

</div>

> [@davidbp](#):
>
> Adjoint{T, AbstractArray{T,2}}

This:

```julia
Adjoint{T, AbstractArray{T,2}}

```

doesn’t mean what you (probably) think it means; you want something like `Adjoint{T, <:AbstractArray{T,2}}` for it to match, see [Types · The Julia Language](https://docs.julialang.org/en/v1/manual/types/#Parametric-Types-1)

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [June 24, 2019, 4:28pm UTC](https://discourse.julialang.org/t/matrix-dense-vector-sparse-product-improvement/25579/8 "2019-06-24T16:28:15Z")

</div>

Thanks for the remainder I was typing too much

```julia
my_prod(W::Adjoint, x::SparseVector) = mat_dense_vec_sparse(W, x) 

```

works fine and much faster than the base operation.

```julia
 # returns true
typeof(W') <: Adjoint

# returns false 
typeof(W') <: Adjoint{AbstractFloat, AbstractArray{AbstractFloat,2}} 

```

I still don’t know how to overload the \* operator in base. I can’t use the same definition as `my_prod` from above. If I try to use that method then

```julia
MethodError: *(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64}) is ambiguous. Candidates:
  *(W::Adjoint, x::SparseVector) in Main at In[48]:1
  *(adjA::Adjoint{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:103
Possible fix, define
  *(::Adjoint{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, ::SparseVector{S,Ti} where Ti<:Integer)

```