Generate random positive definite symmetric matrix

General Usage
,
1

Can Julia or some existing package simulate random covariance (symmetric + positive definite) matrices?

I could do

using LinearAlgebra
function random_covariance_matrix(n::Int)
    x = rand(n, n)
    xsym = 0.5(x + x') # make x symmetric
    return xsym + n*I # diagonally dominant matrices are positive definite
end
A = random_covariance_matrix(3)
isposdef(A)
issymmetric(A)

or something similar, but I’m not sure if diagonal dominance (or something similar) is somewhat more special than a typical covariance matrix.

1 Like
2

How’s this?

X=rand(n,n)
A=X'*X
11 Likes
3

What does “random” mean here? As in, what distribution do you want?

2 Likes
4

I’m considering a multivariate normal. I think @ctkelley’s solution is better than mine, thank you.

5

Maybe you want X = randn(n, n) instead of rand.

3 Likes
6

If you want the distribution to be invariant under orthogonal rotations (just like the uniform distribution over real numbers is invariant under translations), then you want the Haar distribution, which you can sample from via

Q, R = qr(randn(n, n))
Q = Q*Diagonal(sign.(diag(R)))

and using Q.

(See What Is a Random Orthogonal Matrix? – Nick Higham; the above is my translation of the matlab code there, but I didn’t test it!)

Edit: whoops, I completely mixed up the question, my answer is for generating orthogonal matrices, not PSD. They’re connected in my head because you can generate a random PSD matrix by taking such a Q and using it as the matrix of eigenvectors of your PSD matrix, eg

Q * Diagonal(x*rand(n)) * transpose(Q)

where x is how large you want the largest eigenvalue to possibly be.

5 Likes
7

if you aren’t concerned about the values, only that it’s symmetric + positive definite, you can use A = let X = randn(n, n); X * X' + I; end (which is what we use for tests in ChainRules)

2 Likes
8

RandomMatrices.jl

4 Likes
9

… and if you want a uniform distribibution over the positive definite matrices, that’s an LKJ distribution:

2 Likes
10

I don’t think it’s obvious from the docs of RandomMatrices.jl how to solve the problem.
EDIT: Unless you already know random matrix theory.

2 Likes
11

GOE, that is GaussHermite(1) is random symmetric with Gaussian entries. It’s a lot faster, O(n^2) to sample the eigenvalues to answer the question of positive definite ness

eigvalrand(GaussHermite(1), n)

Will give the eigenvalues of an n x n sample.

Yes the package is need of a lot of love

4 Likes