 # Generate a positive definite matrix

How do I randomly generate a positive definite matrix？ What functions can directly implement this step？

Depends on exactly what you want. A trivial way is

``````A = randn(n,n); A = A'*A; A = (A + A')/2
``````

If you want to control the eigenvalues, you can use something like

``````Q, _ = qr(randn(n, n)); D = Diagonal(eigvals); A = Q*D*Q'
``````

where `eigvals` is a vector of eigenvalues.

4 Likes

Aside: having gotten used to JuliaMono (thanks, @cormullion!), the typesetting of that `*` looks really weird and confusing…

2 Likes

Marginal note: I’m very happy if I’ve contributed even a tiny amount to your legendary productivity, Tim. (Even though I’ve seen a complaint about the asterisk. But it was Hacker News, so to be expected .)

5 Likes

Is that `A = (A + A')/2` for numerical reasons? It doesn’t seem to be necessary:

``````As = [(A = randn(n,n); A'*A) for _=1:10_000]
symm_err(A) = maximum(abs.(A-A'))
maximum(symm_err.(As)) # == 0.0
``````

It doesn’t seem to be necessary

Hmm, interesting. I know I’ve found cases where it is (where the result has ulp-level errors), so I just add it by default now. But I get the same thing you do from your test case.

You can map the Cholesky factor of a PD matrix into a vector of reals. The code below uses `TransformVariables.CorrCholeskyFactor`, which maps to a Cholesky factor of the correlation matrix (the diagonal of `U'*U` is ones):

``````using TransformVariables
n = 10
C = CorrCholeskyFactor(n)
U = C(rand(dimension(C)))
σ = abs.(randn(n))
L = U * Diagonal(σ)
A = L' * L
``````

then `A` is guaranteed to be PD, also, all PD matrices of this size will be generated for some random values.

In practice, you would want to work with the factors (eg `L`) for almost all numerical applications.