# 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:** 13\
**Page:** 1

<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, 6:49pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/1 "2024-05-27T18:49:29Z")

</div>

Hi everyone,

I was computing some rather large covariance matrices and noticed that that implementation in Statistics is not multithreaded. TLDR, I want the following to be as fast as possible (ideally, for even larger matrices):

```julia
julia> a = rand(20000,1500);
@time vcov(a);
  8.743414 seconds (1.51 k allocations: 17.235 MiB)

```

I went down the rabbit hole for some other packages that implement this, but it seems to me that neither of them are making use of threads. The only package I found that tried this was VectorizedStatistics.jl, but I found this to be slower in practice than the implementation from Base.

I feel like it can’t be so difficult to make this multithreaded, but I cannot imagine no one has done this before, unless there are good reasons that I am missing.  
I read somewhere here on discourse that one typically does not want to actually use `cov` for large matrices, but rather use other estimators (i.e. from CovarianceEstimators.jl) that implement shrinking. However as far as I can tell they all seem to not be multithreaded either.

Ultimately, I would like to use a Cholesky factorization to solve a linear system. Is this problem so ill-conditioned for large matrices that it is pointless to even try to make it more efficient? So far it was good sufficient to add some small diagonal shift to the covariance matrix.

I hope someone can give me more ideas 🙂

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [May 27, 2024, 7:08pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/2 "2024-05-27T19:08:48Z")

</div>

Perhaps you can comment on the application you have in mind. In some geospatial applications we avoid building the covariance of all locations at once, because that doesn’t fit into memory, or is too slow.

---

<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, 7:42pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/3 "2024-05-27T19:42:24Z")

</div>

Sure, I am experimenting with implementing a stochastic reconfiguration method: [https://arxiv.org/abs/cond-mat/0502553](https://arxiv.org/abs/cond-mat/0502553)  
(see also [Quantum Geometric Tensor and Stochastic Reconfiguration — NetKet](https://netket.readthedocs.io/en/latest/docs/sr.html))

It is somewhat similar to stochastic gradient descend. One tries to optimize a bunch of variational parameters using the quantum geometric tensor, which is essentially a covariance matrix `S = <X_k X_l> - <X_k> <X_l>` where X is a vector that is estimated using a Markov Chain.

The parameters are optimized by solving a linear system `S δα = F`  
To optimize many parameters, one should have many measurements of X\_k so that the covariance matrix is well approximated. However, actually computing the covariance matrix is rather slow.

As you say, I can imagine that practice, this is done very differently. I am absolutely no expert on this to be honest. Perhaps one uses a lazy representation of the covariance matrix instead?  
I hope this clarifies things 🙂

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 27, 2024, 8:03pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/4 "2024-05-27T20:03:00Z")

</div>

> [@Salmon](#):
>
> I was computing some rather large covariance matrices and noticed that that implementation in Statistics is not multithreaded.

The [`cov` implementation in Statistics](https://github.com/JuliaStats/Statistics.jl/blob/9addbb80bb28309a73bf214fc974dcb2dbad3276/src/Statistics.jl#L525-L587) _is_ multithreaded, because the slowest part of the `cov` computation is a matrix multiplication, which uses BLAS … and BLAS is multithreaded (as well as SIMD-vectorized).

In particular, if you have an m \times n matrix A, then to compute the covarance matrix you first subtract the column mean to obtain X, which has \Theta(mn) cost; then compute X^\* X, which has \Theta(mn^2) cost; then do [some normalization](https://en.wikipedia.org/wiki/Bessel%27s_correction), which has \Theta(n^2) cost. So, in cases like yours where m and n are large, the X^\* X computation should completely dominate.

You could try different matrix-multiplication libraries, like MKL.jl or Octavian.jl or AppleAccelerate.jl, to see if they are faster than the default OpenBLAS library on your machine. You could also try `Float32` precision to see if that is accurate enough for your application.

> [@Salmon](#):
>
> The only package I found that tried this was [VectorizedStatistics.jl](https://juliahub.com/ui/Packages/General/VectorizedStatistics), but I found this to be slower in practice than the implementation from Base.

That’s not surprising. LoopVectorization.jl is very good, but it’s best on single linear-time (BLAS1-like or maybe BLAS2-like) loops; higher-level transformations (like in Octavian.jl) are required to make matrix multiplication (BLAS3) fast enough to compete with optimized BLAS libraries for large matrices.

---

<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:15pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/5 "2024-05-27T20:15:43Z")

</div>

> [@stevengj](#):
>
> The [`cov` implementation in Statistics](https://github.com/JuliaStats/Statistics.jl/blob/9addbb80bb28309a73bf214fc974dcb2dbad3276/src/Statistics.jl#L525-L587) _is_ multithreaded, because the slowest part of the `cov` computation is a matrix multiplication, which uses BLAS … and BLAS is multithreaded.

strange, I only see one thread usage, despite setting the number of BLAS threads to 4 with `BLAS.set_num_threads(4)`.  
As a test, if I make a call to `eigen` for some large matrix, I do see higher thread usage.  
Is there something misconfigured on my end?

```julia
julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × Intel(R) Core(TM) i5-4690 CPU @ 3.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, haswell)
  Threads: 5 on 4 virtual cores
Environment:
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_NUM_THREADS = auto

```

---

<div class="post-metadata">

**Author:** ![emil\_hedevang\_sgre](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@emil\_hedevang\_sgre](https://discourse.julialang.org/u/emil_hedevang_sgre)\
**Post date:** [May 27, 2024, 8:22pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/6 "2024-05-27T20:22:10Z")

</div>

Perhaps it would make sense to solve S\,\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 to 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).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 27, 2024, 8:29pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/7 "2024-05-27T20:29:37Z")

</div>

> [@Salmon](#):
>
> strange, I only see one thread usage

How are you determining thread usage?

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

```julia
julia> using Statistics, BenchmarkTools, LinearAlgebra

julia> a = rand(20000,1500);

julia> BLAS.set_num_threads(1);

julia> @btime cov($a) evals=1;
  968.765 ms (11 allocations: 246.06 MiB)

julia> BLAS.set_num_threads(2);

julia> @btime cov($a) evals=1;
  524.270 ms (11 allocations: 246.06 MiB)

```

---

<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).

---

<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:48pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/9 "2024-05-27T20:48:00Z")

</div>

> [@stevengj](#):
>
> You could try different matrix-multiplication libraries, like [MKL.jl](https://juliahub.com/ui/Packages/General/MKL) or [Octavian.jl](https://juliahub.com/ui/Packages/General/Octavian) or [AppleAccelerate.jl](https://juliahub.com/ui/Packages/General/AppleAccelerate), to see if they are faster than the default OpenBLAS library on your machine. You could also try `Float32` precision to see if that is accurate enough for your application.

This seems to be the answer! Upon switching to MKL.jl, I see the expected result:

```julia
julia> @time cov(a);
  9.637042 seconds (12 allocations: 1.103 GiB, 0.93% gc time)

julia> BLAS.set_num_threads(4)

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

```

For `Float32` this time is even reduced by another factor of 2 (yes, I think the precision loss should not matter here).  
Thanks a bunch!  
I suppose the only remaining mystery is why BLAS does not want to use threads for me, though MKL seems to be faster here anyway.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 27, 2024, 9:35pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/10 "2024-05-27T21:35:41Z")

</div>

Could this be linked to the weird interactions between Julia threads and BLAS threads?

[https://docs.julialang.org/en/v1/manual/performance-tips/#man-multithreading-linear-algebra](https://docs.julialang.org/en/v1/manual/performance-tips/#man-multithreading-linear-algebra)

---

<div class="post-metadata">

**Author:** ![dqeeq](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@dqeeq](https://discourse.julialang.org/u/dqeeq)\
**Post date:** [May 27, 2024, 11:17pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/11 "2024-05-27T23:17:23Z")

</div>

I get performance similar to [this](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/8) whether I use `julia -t 1` or `julia -t 4`. Is that what you mean?

In particular, I get the following times with a single julia thread:

```julia
11.757666 seconds (12 allocations: 1.103 GiB, 0.04% gc time)
11.723894 seconds (12 allocations: 1.103 GiB, 0.12% gc time)

```

and the following times with 4 julia threads:

```julia
11.804927 seconds (12 allocations: 1.103 GiB, 0.05% gc time)
11.698108 seconds (12 allocations: 1.103 GiB, 0.13% gc time)

```

---

<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 28, 2024, 5:03am UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/12 "2024-05-28T05:03:25Z")

</div>

> [@gdalle](#):
>
> Could this be linked to the weird interactions between Julia threads and BLAS threads?

Hm, the page doesnt go too much into detail regarding this, but according to the guide everything should be fine when the number of BLAS threads is set correctly. So possibly being a bug with my BLAS artifact?

> [@dqeeq](#):
>
> I get performance similar to [this](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/8) whether I use `julia -t 1` or `julia -t 4`. Is that what you mean?

yup, that’s also what I get without MKL. Could you share what hardware youre running this on? Maybe thats our common denominator

---

<div class="post-metadata">

**Author:** ![dqeeq](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@dqeeq](https://discourse.julialang.org/u/dqeeq)\
**Post date:** [May 28, 2024, 10:26pm UTC](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/13 "2024-05-28T22:26:57Z")

</div>

I don’t think this is hardware-related