Matrix (dense) vector (sparse) product improvement

davidbp · June 23, 2019, 6:01pm

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

PetrKryslUCSD · June 23, 2019, 6:10pm

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

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

That should help you write the method.

davidbp · June 23, 2019, 6:14pm

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[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)

davidbp · June 23, 2019, 7:20pm

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)
(Adjoint{Float64,Array{Float64,2}}, SparseVector{Float64,Int64})

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

PetrKryslUCSD · June 24, 2019, 4:47am

I would try to do what it suggests.

davidbp · June 24, 2019, 10:25am

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)
(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

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

fredrikekre · June 24, 2019, 10:33am

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

davidbp · June 24, 2019, 4:28pm

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

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)