I got a strange question when I was experimenting some convex optimizations. The question is:

Suppose I randomly (say standard normal distribution) generate a $N \times N$ symmetric matrix, (for example, I generate upper triangular matrix, and fill the bottom half to make sure it is symmetric), what is the chance it is a positive definite matrix? Is there anyway to calculate the probability?

  • 1
    $\begingroup$ Try simulation ... $\endgroup$ – kjetil b halvorsen Jan 8 '18 at 19:02
  • 1
    $\begingroup$ @kjetilbhalvorsen thanks, but I am wondering what's the chance of all the eigenvalues are larger than 0. or can we even possibly to do it analytically. $\endgroup$ – Haitao Du Jan 8 '18 at 19:03
  • 6
    $\begingroup$ The answer depends on how you generate the matrix. For instance, one way generates $n$ real eigenvalues according to some distribution and then conjugates that diagonal matrix by a random orthogonal matrix. The result will be positive definite if and only if all those eigenvalues are positive. If you were to generate eigenvalues independently according to a distribution symmetric about zero, then that chance obviously is at most $2^{-n}$. To generate a PD matrix, then, choose your eigenvalues well! (For quick work, I create such matrices as covariances of multivariate Normal data.) $\endgroup$ – whuber Jan 8 '18 at 19:36
  • 11
    $\begingroup$ Not an answer to the question asked, but note that if first simulate a matrix $L$ with each entry iid normal and the same dimensions of $N$, then $N = LL^T$ is symmetric and positive definite with probability 1. $\endgroup$ – Cliff AB Jan 8 '18 at 19:38

If your matrices are drawn from standard-normal iid entries, the probability of being positive-definite is approximately $p_N\approx 3^{-N^2/4}$, so for example if $N=5$, the chance is 1/1000, and goes down quite fast after that. You can find an extended discussion of this question here.

You can somewhat intuit this answer by accepting that the eigenvalue distribution of your matrix will be approximately Wigner semicircle, which is symmetric about zero. If the eigenvalues were all independent, you'd have a $(1/2)^N$ chance of positive-definiteness by this logic. In reality you get $N^2$ behavior, both due to correlations between eigenvalues and the laws governing large deviations of eigenvalues, specifically the smallest and largest. Specifically, random eigenvalues are very much akin to charged particles, and do not like to be close to each other, hence they repel each-other (strangely enough with the same potential field as charged particles, $\propto 1/r$, where $r$ is the distance between adjacent eigenvalues). Asking them to all be positive would therefore be a very tall request.

Also, because of universality laws in random matrix theory, I strongly suspect the above probability $p_N$ will likely be the same for essentially any "reasonable" random matrix, with iid entries that have finite mean and standard deviation.

  • 5
    $\begingroup$ It is nice to know it is very low. So I will not use rejection sampling to create SPD matrix in the future. $\endgroup$ – Haitao Du Jan 8 '18 at 19:26
  • 5
    $\begingroup$ @hxd1011: if you are trying to sample SPD matrices, I suggest the method I desribed in the comments above. In addition, it may be helpful to read up on Cholesky decompositions $\endgroup$ – Cliff AB Jan 8 '18 at 19:41
  • $\begingroup$ @CliffAB thanks. I usually generate SPD matrix form covariance matrix of some data or from $A'A$ similar to what you suggested. I had the time of trying to manually put some numbers to a small matrix say $2 \times 2$ and hope it is a PD matrix. $\endgroup$ – Haitao Du Jan 8 '18 at 19:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.