Performance issues for matrix multiplication in a loop

[dnmiller]

The Matlab built-in ltitr is a fairly low-level function used heavily by a number of Matlab toolboxes. But even a hand-coded version in Matlab is roughly 10x faster than almost identical Juila:

https://gist.github.com/3962937

Is there anything that can be done about this, or is it more of an LLVM issue? Also, is there any interest in creating a benchmark set sort of like what the Python pandas library has?

[Keno]

Here's what I get:

julia> include("test.jl")
elapsed time: 0.769881010055542 seconds
ltitr took 15.39762020111084 ms

Matlab:

>> ltitr_test
Built-in took 2.0919 ms.
Non-built-in took 32.8778 ms.

Note that I had to fix your Matlab code to actually run the hand-coded version 50 times rather than once. I don't know enough about the ltitr function to say whether Matlab has an algorithmic advantage here, but maybe somebody else can comment on that.

[dnmiller]

My bad. Post-correction, I get very similar results with Julia v. hand-coded Matlab (about 48 vs. 45 ms). This could partially be a LAPACK vs. MKL issue, which would be way beyond the Julia dev's control. I'll try and link Julia into MKL to see what happens. [Edit: I ran into "could not show value of type LoadError" and was unable to resolve it when linking to MKL.]

[ViralBShah]

We need optimizations to the way we do indexing to create fewer temporaries. We also need to improve some of our gc to reuse temporaries, which @JeffBezanson is already working on.

I am reasonably sure that we can get built-in matlab performance in julia by directly calling the BLAS routines.

[andreasnoack]

Built-in Matlab performance is almost achieved if the calls to BLAS are completely avoided. This code

function ltitr2{T}(A::StridedMatrix{T}, B::StridedMatrix{T}, U::StridedMatrix{T})
    P = size(U, 1)
    N = size(U, 2)
    M = size(A, 1)
    X = zeros(M, N)
    size(U, 1)
    for i = 2:N
        for j = 1:M
            for k = 1:M
                X[j,i] += A[j,k]*X[k, i-1]
            end
            for k = 1:P
                X[j,i] += B[j,k]*U[k, i-1]
            end
        end
    end
    return X
end

is almost the same speed as the Matlab build in on my machine. However, the vectorised code it much easier to read so it would a good thing if the calls to BLAS were faster.

[ViralBShah]

Here's a version that should compute ltitr without any copies by using subArrays and calling Blas.gemv. It is quite cool how easy it was to do this, but surprisingly, this is actually slower on your example and on certain larger problems I tried as well.

function ltitr3(A, B, U)
    M = size(A, 1)
    K = size(U, 1)
    N = size(U, 2)
    X = zeros(M, N)
    AXvec = zeros(M)
    BUvec = zeros(M)
    for i = 2:N
        Blas.gemv!('N', 1.0, A, sub(X,1:M,i-1), 0.0, AXvec)
        Blas.gemv!('N', 1.0, B, sub(U,1:K,i-1), 0.0, BUvec)
        X[:, i] = AXvec + BUvec
    end
    return X
end

[ViralBShah]

@andreasnoackjensen I was planning to do just that, but didn't get around to. I verify that it is 8x faster for me.

[ViralBShah]

I would have expected the ltitr3 implementation that uses BLAS to be equally fast, but it isn't. The matrices are perhaps too small for overcoming the BLAS calling overhead. We may even want to have julia implementations of these routines and avoid calling BLAS.

[staticfloat]

@ViralBShah; has it been shown that using BLAS for larger problems (without using extra threads in BLAS) is faster than pure Julia?

[timholy]

It may depend on the routine and the problem size. On problems of modest size, the performance of axpy is simply unreal:

function my_axpy!{T}(alpha::T, x::Array{T}, y::Array{T})
    for i = 1:length(y)
        y[i] += alpha*x[i]
    end
end

function time_axpy(n, niter)
    y = randn(n)
    x = randn(n)
    alpha = randn()
    ycopy = copy(y)
    @time begin
        for i = 1:niter
            my_axpy!(alpha, x, y)
        end
    end
    y = ycopy
    @time begin
        for i = 1:niter
            Blas.axpy!(alpha, x, y)
        end
    end
end


julia> time_axpy(1000,1000)
elapsed time: 0.01027989387512207 seconds
elapsed time: 0.001661062240600586 seconds

julia> time_axpy(10000,1000)
elapsed time: 0.05250382423400879 seconds
elapsed time: 0.005312919616699219 seconds

julia> time_axpy(100000,1000)
elapsed time: 0.4226851463317871 seconds
elapsed time: 0.08148813247680664 seconds

julia> time_axpy(1000000,100)
elapsed time: 0.4845619201660156 seconds
elapsed time: 0.28545117378234863 seconds

julia> time_axpy(10000000,10)
elapsed time: 0.487393856048584 seconds
elapsed time: 0.28293299674987793 seconds

It's interesting that the advantage declines for larger problems. I wouldn't be surprised if, with bounds-checking off, the Julia version would be the same for large problem sizes.

[JeffBezanson]

The julia version wins only for tiny vectors:

julia> time_axpy(3,10000000)
elapsed time: 0.3249959945678711 seconds
elapsed time: 0.5162699222564697 seconds

[JeffBezanson]

Should we add ltitr2 to the library?

[andreasnoack]

I did some more benchmarking and it seems like OpenBLAS is particularly bad at this problem. On my Mac, the original implementation which uses BLAS spends ~48msec with OpenBlas vs ~30msec for vecLib but on Ubuntu in a VirtualBox the timing is ~105msec for OpenBLAS vs ~35msec for Netlib BLAS.

Then I wrote a version in Fortran that called BLAS and called it from Julia. On my Mac it was ~16msec for OpenBlas vs ~4msec for vecLib.

Another thing. I wondered if it was a good idea to avoid BLAS for very small multiplications like those in ltitr and wrote a special function for matrix multiplications in pure Julia but if I use that instead of the inner double loop in ltitr2 the function becomes much slower. Will it always be like that or will Julia eventually be able to optimise code "exported" to functions. It easily becomes a mess if you write every matrix multiplication out as double loops.

[ViralBShah]

I have seen performance problems on matrices with large aspect ratios across a variety of BLAS libraries over the years. There is no fixed rule, but writing some of these in julia is the right idea. We already have base/matmul.jl where we already have some special cases in julia.

Calling even a julia routine in the inner loop results in the allocation of a new matrix in every iteration of the loop. That is the real reason of the slowdown. Matlab for example, tends to reuse temporaries across loop iterations, and @JeffBezanson is working on similar improvements in the jb/gcwork branch. Inlining in julia works well in general, and "exporting" inner loops to functions should end up having no performance difference in most cases.

[ViralBShah]

Also, report the performance issue to the OpenBLAS folks.

Cc: @xianyi