Four ways to extract a column

[David Roodman]

This is more a comment than a question, as they say. I discovered some years ago that different syntaxes for extracting a column from a matrix have different run times. In the test below, for matrices with less than about 10,000 rows, it's faster to precompute a colvector one2N=1::N, and then do M[one2N, 3] instead of M[.,3] or M[,3] or M[|.,3\3,.|]. I found this surprising. But if one wants to extract several columns, then it seems best to do something like M[., (3,5)].

I found this to be true on a Lenovo T430 laptop with an i7 processor in Stata/IC, and on an old HP workstation with 12 cores and 12-core Stata MP.

Here's some test code.

cap mata mata drop test()
mata
mata set matastrict on
mata set mataoptimize on
mata set matalnum off

// N = rows of matrix to extract from; iter = number of repetitions; p = rowvector of cols to extract
void test(real scalar N, real scalar iter, real rowvector p) {
    real colvector one2N
    real matrix M, S, t
    real scalar i

    M = runiform(N, 5)
    one2N = 1::N

    timer_clear()
    timer_on(1)
    for (i=iter;i;i--) t = M[,p]
    timer_off(1)
    timer_on(2)
    for (i=iter;i;i--) t = M[.,p]
    timer_off(2)
    timer_on(3)
    for (i=iter;i;i--) t = M[one2N,p]
    timer_off(3)
    if (cols(p)==1) {
        S = .,p\.,p
        timer_on(4)
        for (i=iter;i;i--) t = M[|S|]
        timer_off(4)
    }
    timer()
}
test(10,1000000 ,3)
test(100,1000000,3)
test(1000,100000,3)
test(10000,10000,3)
test(1000000,100,3)

mata
test(10,1000000 ,(3,5))
test(100,1000000,(3,5))
test(1000,100000,(3,5))
test(10000,10000,(3,5))
test(1000000,100,(3,5))
end

[Sergio Correia]

Wow, that's really interesting.

I also found some interesting results, documented mostly for myself in bullet 5 here: http://scorreia.com/blog/2016/10/06/mata-tips.html . Some are expected (loop unrolling, common cases first, but others surprised me).

For instance, evaluating x[i] is much faster than evaluating x[i,.] , even if x only has one column. Thus, I often have code like return(k > 1 ? x[i, .] : x[i])

Perhaps these should be optimized away by the compiler?