# Numerical inaccuracies when constructing correlation matrices

**URL:** <https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [August 2, 2022, 4:18pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177 "2022-08-02T16:18:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![daniel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel/32/6480_2.png) [@daniel](https://discourse.julialang.org/u/daniel)\
**Post date:** [August 2, 2022, 4:18pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/1 "2022-08-02T16:18:16Z")

</div>

I’m trying to sample correlation matrices using the vine method explained here:  
[How to generate a large full-rank random correlation matrix with some strong correlations present? - Cross Validated](https://stats.stackexchange.com/a/125020) (which is based on the JKL 2009 paper [https://doi.org/10.1016/j.jmva.2009.04.008](https://doi.org/10.1016/j.jmva.2009.04.008)).

```julia
using Random
using LinearAlgebra
Random.seed!(123)

function random_correlation_matrix(d::Int)
    P = fill(0.0, d, d)
    S = Matrix{Float64}(I, d, d)

    for k in 1:(d - 1), i in (k + 1):d
        P[k, i] = (rand() - 0.5) * 2 # scale to [-1, 1]

        p = P[k, i]
        for l in (k - 1):-1:1
            p = p * sqrt((1 - P[l, i] ^ 2) * (1 - P[l, k] ^ 2)) + P[l, i] * P[l, k]
        end
        
        S[k, i] = p
        S[i, k] = p
    end

    permutation = shuffle(1:d)
    S = S[permutation, permutation]

    return S
end

S1 = random_correlation_matrix(10);
S2 = random_correlation_matrix(1000);

```

This works fine for lower dimensions `d`, but fails for higher dimensions:

```julia
isposdef(S1) # true
isposdef(S2) # false

```

This seems to be due to accumulating floating-point rounding errors, as the smallest eigenvalues are all close to 0 but negative:

```julia
minimum(eigvals(S2)) # -8.46e-14

```

Is there a method to heal an _almost_ positive semidefinite matrix like this (or tweaks to make the algorithm more robust)?

I’ve found a `cor_nearPD` function in the [Bigsimr.jl](https://schisslergroup.github.io/Bigsimr.jl/dev/nearest_correlation_matrix/) package that does exactly what I want, but the problem seems to be simple enough that I’d like to avoid adding another dependency.

I also tried doing an eigendecomposition and simply replace the negative eigenvalues in D, but that didn’t work out because I’m lacking the background in linear algebra to know what I’m doing.

---

<div class="post-metadata">

**Author:** ![stephenll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephenll/32/28751_2.png) [@stephenll](https://discourse.julialang.org/u/stephenll)\
**Post date:** [August 2, 2022, 5:56pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/2 "2022-08-02T17:56:31Z")

</div>

Here’s code I use. Not the greatest code, needs some updating for clarity.

```julia
using LinearAlgebra
using 
function quickpd!(x::AbstractMatrix, posdtol = 1e-8)

    n = size(x,1)

    diagIdx = 1:n+1:n^2

   # if !issymmetric(x)
        # TODO: add error
    #end

    vals, vecs = eigen(x)

    eips = posdtol * maximum(abs.(vals))

    vals[vals .< eips] .= eips

    o_diag = diag(x)

    X = vecs * (vals .* transpose(vecs))

    D = sqrt.(max.(eips, o_diag) ./ diag(X))
    
    x .= D .* X .* repeat(transpose(D), inner=(n,1))

    x[diagIdx] .= 1.0

    return x
end

```

---

<div class="post-metadata">

**Author:** ![joshualeond](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joshualeond/32/230_2.png) [@joshualeond](https://discourse.julialang.org/u/joshualeond)\
**Post date:** [August 2, 2022, 6:35pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/3 "2022-08-02T18:35:55Z")

</div>

Is the _vine method_ different than random sampling from an [LKJ](https://juliastats.org/Distributions.jl/v0.25/matrix/#Distributions.LKJ) distribution in `Distributions.jl`? I know you said you don’t want to introduce a dependency but just curious:

```julia
using Distributions, Random, LinearAlgebra
Random.seed!(123)

S1 = rand(LKJ(10, 1))
S2 = rand(LKJ(1000, 1))

```

---

<div class="post-metadata">

**Author:** ![daniel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel/32/6480_2.png) [@daniel](https://discourse.julialang.org/u/daniel)\
**Post date:** [August 2, 2022, 6:50pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/4 "2022-08-02T18:50:43Z")

</div>

Cool, thanks! I tried it out and it actually works

```julia
isposdef(quickpd!(copy(S2))) # false, but adding Symmetric:
isposdef(Symmetric(quickpd!(copy(S2)))) # true

```

Could you maybe explain what’s going on in these lines (or give a link) so I can brush up my linear algebra?

```julia
    D = sqrt.(max.(eips, o_diag) ./ diag(X))
    x .= D .* X .* repeat(transpose(D), inner=(n,1))

```

---

<div class="post-metadata">

**Author:** ![daniel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel/32/6480_2.png) [@daniel](https://discourse.julialang.org/u/daniel)\
**Post date:** [August 2, 2022, 7:01pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/5 "2022-08-02T19:01:32Z")

</div>

Yes, from my (limited) understanding using the LKJ distribution is indeed identical to this function (not sure about the eta=1 though).  
The authors of the paper linked above are where the initials for the LKJ distribution come from.

My usecase is a bit more generic though, I want the partial correlation coefficients to come from a mixture of Beta distributions, so that I can simulate correlation matrices where some of the correlation coefficients are large and some low (I just left this part out of the minimal example because it wasn’t relevant).

---

<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:** [August 2, 2022, 7:27pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/6 "2022-08-02T19:27:43Z")

</div>

> [@daniel](#):
>
> Is there a method to heal an _almost_ positive semidefinite matrix like this (or tweaks to make the algorithm more robust)?

The basic issue is that trying to make a matrix _exactly_ semidefinite is inherently nonrobust to roundoff errors. Different ways of computing the eigenvalues will give slightly different values due to differences in floating-point effects, so even if you make the eigenvalues \ge 0 for one computation you might get slightly negative eigenvalues from another computation/algorithm with the same matrix.

The question is, what are you trying to _do_ with the matrix `S2` that requires it to be strictly semidefinite? It is generally better to try to fix your _subsequent_ computation to be tolerant of tiny negative eigenvalues in such cases. For example, `sqrt(Symmetric(A))` includes a backwards-stable correction to return a real matrix square root for nearly semidefinite matrices (see [julia#35057](https://github.com/JuliaLang/julia/pull/35057)).

Alternatively, if you don’t care too much about accuracy, you can just do something like `S2 + 1e-8*norm(S2)*I`, of course (sacrificing half of the significant digits).

You should also consider the possibility that you have a bug in your code, or are mis-using the algorithm, or that the source you are using is incorrect. Your linked article claims that the algorithm generates “full-rank” correlation matrices, but your `S2` matrix is extremely rank deficient: it is a 1000x1000 matrix with numerical rank (number of non-negligible eigenvalues) of only **64**. Is this really what you want?

---

<div class="post-metadata">

**Author:** ![daniel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel/32/6480_2.png) [@daniel](https://discourse.julialang.org/u/daniel)\
**Post date:** [August 2, 2022, 9:59pm UTC](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/7 "2022-08-02T21:59:56Z")

</div>

Thank you! It does not matter to me that the matrix is strictly semidefinite (that’s impossible anyway), I just want it to behave nice enough that `isposdef(A) = true`, so I can do stuff like sampling from a multivariate normal without getting a `PosDefException`. The ultimate goal is to use it for a Gaussian copula.

I believe the full rank thing just comes from the stack exchange post. It is of course entirely possible that there is a bug in the code, I’m misusing the algorithm or that the person in that stack exchange post did something wrong, it might be even all of that!

But for now I’m happy with getting correlation matrices where some elements are high and others low, and the LKJ paper seemed like a good fit because it allows me to set the partial correlations to whatever I want.