# Matrix (dense) vector (sparse) product improvement

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

``````
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
``````
``````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

``````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
``````
``````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
``````

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

``````julia> typeof(A')
``````

That should help you write the method.

I tried the following:

``````import Base
import Base.*

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

But I get…

``````
MethodError: *(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64}) is ambiguous. Candidates:
*(W::Adjoint, x::SparseVector) in Main at In: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)
``````

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.

``````
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)
``````
``````typeof(W'), typeof(x_sp)
``````
``````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: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:1
``````

I would try to do what it suggests.

Well isn’t the following definition what was suggested?

``````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

``````typeof(W'), typeof(x_sp)
``````

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

``````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:1
``````

This:

``````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

1 Like

Thanks for the remainder I was typing too much

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

works fine and much faster than the base operation.

`````` # returns true
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
``````MethodError: *(::Adjoint{Float64,Array{Float64,2}}, ::SparseVector{Float64,Int64}) is ambiguous. Candidates: