# Generate random positive definite symmetric matrix

**URL:** <https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816>\
**Category:** General Usage\
**Tags:** question, statistics\
**Created:** [January 23, 2021, 1:34am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816 "2021-01-23T01:34:13Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![biona001](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/biona001/32/16497_2.png) [@biona001](https://discourse.julialang.org/u/biona001)\
**Post date:** [January 23, 2021, 1:34am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/1 "2021-01-23T01:34:13Z")

</div>

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

I could do

```julia
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.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [January 23, 2021, 1:57am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/2 "2021-01-23T01:57:10Z")

</div>

How’s this?

```julia
X=rand(n,n)
A=X'*X

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [January 23, 2021, 2:35am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/3 "2021-01-23T02:35:40Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![biona001](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/biona001/32/16497_2.png) [@biona001](https://discourse.julialang.org/u/biona001)\
**Post date:** [January 23, 2021, 2:54am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/4 "2021-01-23T02:54:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [January 23, 2021, 3:05am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/5 "2021-01-23T03:05:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [January 23, 2021, 11:05am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/6 "2021-01-23T11:05:42Z")

</div>

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

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

```

and using `Q`.

(See [What Is a Random Orthogonal Matrix? – Nick Higham](https://nhigham.com/2020/04/22/what-is-a-random-orthogonal-matrix/); 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

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

```

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

---

<div class="post-metadata">

**Author:** ![npr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/npr/32/10626_2.png) [@npr](https://discourse.julialang.org/u/npr)\
**Post date:** [January 23, 2021, 11:58am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/7 "2021-01-23T11:58:38Z")

</div>

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](https://github.com/JuliaDiff/ChainRulesTestUtils.jl/blob/478c2cba8f36dec8ec9e14183ae7b75ac7f17ff5/src/data_generation.jl))

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 23, 2021, 3:59pm UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/8 "2021-01-23T15:59:21Z")

</div>

RandomMatrices.jl

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [January 23, 2021, 4:38pm UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/9 "2021-01-23T16:38:44Z")

</div>

> [@ericphanson](#):
>
> 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

… and if you want a uniform distribibution over the positive definite matrices, that’s an [LKJ distribution](https://juliastats.org/Distributions.jl/stable/matrix/#Distributions.LKJ):

> **[Is the LKJ(1) prior uniform? "Yes" | Stephen R. Martin, PhD](https://srmart.in/is-the-lkj1-prior-uniform-yes/)**
>
> Est. reading time: 8 minutes

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [January 24, 2021, 6:21am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/10 "2021-01-24T06:21:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 24, 2021, 10:19am UTC](https://discourse.julialang.org/t/generate-random-positive-definite-symmetric-matrix/53816/11 "2021-01-24T10:19:38Z")

</div>

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

```nohighlight
eigvalrand(GaussHermite(1), n)

```

Will give the eigenvalues of an `n x n` sample.

Yes the package is need of a lot of love