[julia-users] Re: Linear combination of matrices using repmat: Matlab vs Julia

[Steven G. Johnson]

On Tuesday, November 24, 2015 at 12:01:32 PM UTC-5, Pieterjan Robbe wrote:
>
> Thanks for the fast response, that explains a lot. I don't see the 7x 
> speedup on my machine though, but probably the problem is too small to see 
> significant differences.
>

That's odd, I'm benchmarking the exact same size you are, and the 
difference is very clear:

function matrixLinComb(A::Array{Float64,3},y::Array{Float64,2})
  Z = reshape(sum(reshape(y,1,1,m).*A,3),n,n)
end

function matrixLinComb2{T1,T2}(A::Array{T1,3},y::Array{T2})
    m,n,p = size(A)
    length(y) == p || throw(ArgumentError("y should be length $p"))
    C = zeros(promote_type(T1,T2), m,n)
    for k = 1:p
        yk = y[k]
        for j = 1:n, i = 1:m
            @inbounds C[i,j] += yk * A[i,j,k]
        end
    end
    return C
end

n = 100
m = 10000

A = rand(n,n,m)
y = randn(1,m)

println("matrixLinComb: ", minimum([@elapsed(matrixLinComb(A,y)) for i = 
1:10]))
println("matrixLinComb2: ", minimum([@elapsed(matrixLinComb2(A,y)) for i = 
1:10])) 


gives

matrixLinComb: 0.633029344

matrixLinComb2: 0.084245998


Note that I'm using the minimum time here, not the mean, so it 
automatically ignores JIT times and other one-time spikes.