Consider the problem of taking a linear combination of m (n x n)-matrices
stored in a (n x n x m)-array A. The weights are stored in a length-m vector
w.
In Matlab, we can accomplish this by
n = 100;
m = 10000;
A = rand(n,n,m);
y = randn(1,m);
times = zeros(1,10);
for cntr = 1:10
tic;
Z = sum(repmat(reshape(y,[1 1 m]),[n n 1]).*A,3);
times(cntr) = toc;
end
sprintf('time = %f',mean(times))
This gives time = 0.3259 on my machine using R2014b. (Note: I am aware of
the ridiculous amount of memory allocated by this solution, however, it
turns out to be 25% faster than using bsxfun, for instance)
A Julia-implementation that accomplishes the same thing is
function matrixLinComb(A::Array{Float64,3},y::Array{Float64,2})
Z = reshape(sum(reshape(y,1,1,m).*A,3),n,n)
end
n = 100
m = 10000
A = rand(n,n,m)
y = randn(1,m)
times = zeros(10)
for cntr = 1:10
times[cntr] = @elapsed matrixLinComb(A,y)
end
println("time = $(mean(times))")
This gives time = 0.5130 using v0.4.0 . I've tried out different other
implementations (broadcast, reshaping A and *multiply etc.), though this
turned out the be the fastest method. Does anyone has an idea why Matlab's
repmat function is superior to the Julia-implementation (or, as
my colleagues call it - Julia doesn't work)? Thanks!
