Most efficient implementation of covariance matrix

[Salmon]

@stevengj

How are you determining thread usage?

If I simply look at the performance, the impact of multiple threads is clearly apparent

Not for me, unfortunately:

julia> a = rand(4800,10000);
julia> BLAS.set_num_threads(4)

julia> @time cov(a);
 11.101160 seconds (12 allocations: 1.103 GiB)

julia> BLAS.set_num_threads(1)

julia> @time cov(a);
 10.656343 seconds (12 allocations: 1.103 GiB, 0.59% gc time)

I can only guess there is some issue with the BLAS call.

@emil_hedevang_sgre

Perhaps it would make sense to solve S,\delta\alpha=FSδα=FS,\delta\alpha=F using an interative method. In this case you don’t need to materialize the covariance matrix, you just need the result of the linear transformation that maps \delta\alphaδα\delta\alpha to S,\delta\alphaSδαS,\delta\alpha. Up to some fiddling, this map basically reduces to two matrix-vector multiplications. Have a look at LinearMaps.jl. (Please be warned that I haven’t thought through the details).

This is an interesting idea, thanks for the link. I believe this will not be faster here, because I expect that evaluating the matrix vector product once will likely have to actually compute each matrix element (though its probably a good idea to check this, perhaps there is a shortcut I am overlooking).