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