# How to randomly generate a positive semidefinite matrix?

How would I randomly generate a positive semidefinite matrix? I'm aware how to generate a random $$n\times n$$ matrix with real values between -1 and 1 with

SeedRandom
KK=RandomReal[{-1,1},{n,n}];


But how do I generate one with the extra PSD constraint?

• Maybe generate eigenvalues/eigenvectors separately and then build the matrix ? – b.gates.you.know.what Mar 11 at 13:59
• How would I code this? – Glassjawed Mar 11 at 14:00

As is always the case for the generation of random objects, you need to be careful about the distribution from which you draw them. In the case of random positive semi-definite matrices I would try to draw them from a Haar measure, meaning that they should be drawn from a distribution that is invariant under unitary/orthogonal transformations. This can be achieved in @kglr's solution by drawing the random numbers from a normal distribution instead of a hypercube:

rpsdmH[n_] := Module[{k},
While[Not[PositiveSemidefiniteMatrixQ[
k = RandomVariate[NormalDistribution[], {n, n}]]]]; k]


For $$2\times2$$ matrices we can easily check that the eigenvectors are now uniformly distributed (i.e., there are no preferred axes):

Histogram[Table[Mod[ArcTan @@ Eigenvectors[rpsdmH][], π], {10^5}],
{0, π, π/100}] This method is, however, very slow because the probability of hitting a positive semi-definite matrix decreases exponentially with n:

rpsdmH // AbsoluteTiming // First
(*    2.35343    *)


A much more efficient way is to take a random $$n\times n$$ matrix and square it, so that all eigenvalues will be nonnegative:

rpsdmH[n_] := Transpose[#].# &[RandomVariate[NormalDistribution[], {n,n}]]

rpsdmH // AbsoluteTiming // First
(*    0.040501    *)


More natural would be to generate complex-valued matrices with the same trick:

rcpsdmH[n_] := ConjugateTranspose[#].# &[RandomVariate[NormalDistribution[], {n,n,2}].{1,I}]

rcpsdmH
(*    {{11.3155 + 0. I,      1.7642 - 1.58122 I,   3.73334 - 3.09205 I},
{1.7642 + 1.58122 I,  5.28292 + 0. I,       2.12236 - 0.146192 I},
{3.73334 + 3.09205 I, 2.12236 + 0.146192 I, 3.85174 + 0. I}}    *)

Eigenvalues[%]
(*    {14.7023, 4.51453, 1.23333}    *)


This construction still leaves open the distribution of the scale of the generated matrices (expressed as the histogram of traces or determinants). You may need some "radial" scaling to achieve your goals. Often what is needed is random PSD matrices with unit trace, which you can get with

rcpsdmH1[n_] := #/Tr[#] &[rcpsdmH[n]]

rcpsdmH1
(*    {{0.130678 + 0. I,        -0.233684 - 0.105608 I, 0.115911 + 0.124215 I},
{-0.233684 + 0.105608 I, 0.598174 + 0. I,        -0.306882 - 0.115144 I},
{0.115911 - 0.124215 I,  -0.306882 + 0.115144 I, 0.271148 + 0. I}}    *)


Here's the distributions of the smallest, middle, and largest eigenvalues for such randomly generated $$3\times3$$ complex matrices:

Histogram[Transpose[Table[Sort[Eigenvalues[rcpsdmH1]], {10^5}]],
{0, 1, 1/100}] For further "radial" scaling you can apply many functions to the generated matrices. As an example, take the "matrix square-root" MatrixPower[#, 1/2], which is well-defined for PSD matrices and somewhat undoes the squaring operation used to construct them. The distributions of the eigenvalues are now more uniform, as shown in this example for $$5\times5$$ matrices:

Histogram[Transpose[Table[Sort[Eigenvalues[MatrixPower[rcpsdmH1, 1/2]]], {10^5}]],
{0, 1, 1/100}] ClearAll[rpsdm]
rpsdm[n_] := Module[{k},
While[Not[PositiveSemidefiniteMatrixQ[k = RandomReal[{-1, 1}, {n, n}]]]]; k]


Examples:

Row[Panel[MatrixForm @ #,
Row[{"PositiveSemidefiniteMatrixQ:", PositiveSemidefiniteMatrixQ @ #}]] & /@
Table[rpsdm, 3], Spacer] If I don't care very much about the distribution, but just want a symmetric positive-definite matrix (e.g. for software test or demonstration purposes), I do something like this:

m = RandomReal[NormalDistribution[], {4, 4}];

p = m.Transpose[m];

SymmetricMatrixQ[p]
(* True *)

Eigenvalues[p]
(* {9.41105, 4.52997, 0.728631, 0.112682} *)


If I want positive semi-definite, this is easily achievable too:

m = RandomReal[NormalDistribution[], {4, 3}];

p = m.Transpose[m];

SymmetricMatrixQ[p]
(* True *)

Eigenvalues[p] // Chop
(* {7.05972, 4.62517, 0.357622, 0} *)