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.

Aside: having gotten used to JuliaMono (thanks, @cormullion!), the typesetting of that `*`

looks really weird and confusing…

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

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.