# Most efficient implementation of covariance matrix

**URL:** https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797
**Category:** Performance
**Tags:** statistics, multithreading, nerdsnipe, cov, covariance-matrix
**Created:** [May 27, 2024, 6:49pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797 "2024-05-27T18:49:29Z")
**Posts on this page:** 1
**Showing post:** 8

<div class="post-metadata">

### Author: ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)
#### Post date: [May 27, 2024, 8:36pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/8 "2024-05-27T20:36:54Z")

</div>

> [@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
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](https://github.com/JuliaLinearAlgebra/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).

---

_[View the full topic](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797)._