Numerical inaccuracies when constructing correlation matrices

General Usage
2

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

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
1 Like