# The probability for a symmetric matrix to be positive definite

Let me give a reasonable model for the question in the title. In ${\rm Sym}_n({\mathbb R})$, the positive definite matrices form a convex cone $S_n^+$. The probability I have in mind is the ratio $p_n=\theta_n/\omega_n$, where $\theta_n$ is the solid angle of $\Lambda_n$, and $\omega_n$ is the solid angle of the whole space ${\rm Sym}_n$ (the area of the unit sphere of dimension $N-1$ where $N=\frac{n(n+1)}2$). These definitions are relative to the Euclidian norm $\|M\|=\sqrt{{\rm Tr}(M^2)}$ ; this is the most natural among Euclidian norms, because it is invariant under unitary conjugation.

Because $S_2^+$ is a circular cone, I could compute $p_2=\frac{2-\sqrt2}4\sim0.146$ . Is there a known close formula for $p_n$? If not, is there a known asymptotics?

More generally, we may define open convex cones $$\Lambda_n^0\subset\Lambda_n^1\subset\cdots\subset\Lambda_n^{n-1}$$ in the following way: $M\mapsto\det M$ is a homogeneous polynomial, hyperbolic in the direction of the identity matrix $I_n$. Thus its successive derivatives in this direction are hyperbolic too. The $k$th derivative defines a "future cone" $\Lambda_n^k$, these cones being nested. For instance, $\Lambda_n^0=S_n^+$. It turns out that this derivative is, up to a constant, $\sigma_{n-k}(\vec\lambda)$, where $\sigma_j$ is the $j$th elementary symmetric polynomial and $\vec\lambda$ the spectrum of $M$. Therefore $\Lambda_n^k$ is defined by the inequalities $$\sigma_1(\vec\lambda)\ge0,\ldots,\sigma_{n-k}(\vec\lambda)\ge0.$$ For instance, $\Lambda_n^{n-1}$ is the half-space defined by ${\rm Tr}M\ge0$.

Let us define again $p_{n,k}$ the probability for $M\in{\rm Sym}_n$ to belong to $\Lambda_n^k$. Thus $p_{n,0}=p_n$ and $p_{n,n-1}=\frac12$.

What is the distribution of $(p_{n,0},\ldots,p_{n,n-1})$, asymptotically as $n\rightarrow+\infty$?

Edit. As Mikaël mentionned, it is equivalent, and easier for calculations, to consider the standard Gaussian measure (GOE) over ${\bf Sym}_n$.

• If you want to use a more standard terminology (to look for references for examples), you should say that the probability is with respect to the GOE random matrix model (=gaussian vector in $Sym(\mathbb R)$ with covariance $Id$). Hence, $p_n$ is the probability that a GOE matrix is positive. There are exact formulas for the eigenvalue distribution of a GOE. You can perhaps derive from them close formulas for $p_n$. The asymptotic behaviour of $p_n$ is also probably known (look for large deviations results for GOE). – Mikael de la Salle Jan 9 '13 at 23:24
• I don't know about exact values or precise asymptotics, but from Mikael's observation it's not too much work to see that $p_n$ goes to 0 exponentially fast. – Mark Meckes Jan 10 '13 at 15:11
• I don't understand your second question. The $p_{n,k}$ are deterministic, so what do you mean by their distribution? – Mark Meckes Jan 10 '13 at 15:13
• @Mark. Given $n$, you may plot $p_{n,k}$ as a function of $\frac{k}n$. Then, as $n\rightarrow+\infty$, the plots must tend to some graph, and I am interested in this graph. – Denis Serre Jan 10 '13 at 16:32
• "this is the most natural norm, because it is invariant under unitary conjugation": actually there are lots of norms invariant under unitary conjugation: the Frobenius norm, the Euclidean induced norm, Schatten norms, Ky Fan norms... – Federico Poloni Mar 15 '13 at 13:47

Edit: According to Dean and Majumdar, the precise value of $c$ in my answer below is $c=\frac{\log 3}{4}$ (and $c=\frac{\log 3}{2}$ for GUE random matrices). I did not read their argument, but I have been told that it can be considered as rigourous. I heard about this result through the recent work of Gayet and Welschinger on the mean Betti number of random hypersurfaces. I am a bit surprised that this computation was not made before 2008.

Let me just expand my comment. You are talking about the uniform measure on the unit sphere of the euclidean space $Sym_n(\mathbb R)$, but for measuring subsets that are homogeneous it is equivalent to talk about the standard gaussian measure on $Sym_n(\mathbb R)$. This measure is called in random matrix theory the Gaussian Orthogonal Ensemble (GOE). In particular $p_n$ is the probability that a matrix in the GOE is positive definite. Since there are explicit formulas for the probability distribution of the eigenvalues of a GOE matrix (this is probably what Robert Bryant is proving), there migth be explicit formulas for $p_n$.

Anyway, the asymptotics are known from general large deviation results for random matrices (due to Ben Arous and Guionnet, PTRF 1997)~: $p_n$ goes to zero as $e^{-c n^2}$ for some constant $c>0$. The constant is equal to the infimum, over all probability measures $\mu$ on $\mathbb R^+$, of the quantity $$\frac{1}{2} (\int x^2 d\mu(x) - \Sigma(\mu)) - \frac 3 8 - \frac 1 4 \log 2$$ where $\Sigma(\mu)$ is Voiculescu's free entropy $\iint \log|x-y| d\mu(x) d\mu(y)$. You can probably explicitely compute $c$. It is even possible that this was known before Ben Arous and Guionnet's work, since their results are much more general.

For your second question, I am pretty sure that the limiting graph of $t \in [0,1( \mapsto p_{n,E(tn)}$ is $0$ ($E(x)$ is the integer part of $x$). But this is probably not what you really want to ask.

• You're right. If the limiting graph is identically zero (something I anticipated), then I am interesting in a non-trivial asymptotics, when we rescale this graph appropriately. – Denis Serre Jan 10 '13 at 18:20
• My guess (no proof) is that the non-tivial asymptotics will occur for the n-uple $(\log(p_{n,i})/n^2)_{i=0}^n$. This should also follow from the Ben Arous-Guionnet large deviation results. – Mikael de la Salle Jan 10 '13 at 19:08
• One strange corollary of this: Imagine exposing your matrix minor by minor (so that after step $k$ the $k \times k$ upper left submatrix is exposed). By Sylvester's criterion, $M$ is positive definite iff the determinants of each exposed submatrix are positive. So what this is saying is that the probability the $n^{th}$ determinant is positive, conditioned on the previous determinants being positive, decays exponentially in $n$. I find this counterintuitive, especially given that individual entries have symmetric distribution. – Kevin P. Costello Jan 10 '13 at 20:48

This is really a long comment: I don't know the answers, but I expect that these numbers have been calculated by the geometric probabilists (and physicists) quite some time ago.

If you use the fact that the group $\mathrm{SO}(n)$ acts on the symmetric matrices in the obvious way and parametrize the 'hemisphere' of symmetric matrices of positive trace and norm $1$ by the symmetric matrices of zero trace via the $\mathrm{SO}(n)$-equivariant mapping $$P = \frac{(I+z)}{(1+|z|^2)^{1/2}},$$ then you can cover the symmetric matrices of trace zero by acting by $\mathrm{SO}(n)$ on the space $D_n\simeq\mathbb{R}^{n-1}$ of trace-free diagonal matrices by setting $$z= A\ r\ A^T$$ where $r = \mathrm{diag}(r_1,\ldots,r_n)$ with $r_1+\cdots+r_n=0$. Then you can compute that the volume form on the hemisphere of symmetric matrices of positive trace and norm $1$ pulls back under this mapping of $\mathrm{SO}(n)\times D_n\to{\mathrm{Sym}}_n({\mathbb{R}})$ to be a constant times the volume form on $\mathrm{SO}(n)$ times $$\rho = \frac{\prod_{i<j}(r_i{-}r_j)}{(1+r_1^2+\cdots+r_n^2)^{n(n+1)/2}}\ dr_1\wedge dr_2\wedge\cdots dr_{n-1}\ .$$ (Remember that $r_n=-(r_1+\cdots+r_{n-1})$.) Thus, the probability you want is one-half the integral of $\rho$ over the whole of $D_n$ divided into the integral of $D_n$ over the $(n{-}1)$-simplex defined by the $n$ inequalities $$r_i+1 > 0.$$ This looks like the sort of integral that one encounters in the study of random matrices in the GOE theory, so I'd check there to see how these are evaluated.

The other cones that you mention will define other regions in $D_n$, and it's possible that there are methods for computing these integrals as well.

I came across this question while preparing a talk on this paper. There we gave the explicit formula for $p_n$ (in the GOE) that the other answers anticipated could be found by random matrix methods.

The formula is given as follows:

First, let $n\geq 1$ be any integer, and define $n':=2 \lceil n/2 \rceil$. When real symmetric matrices are chosen according to the $n$-dimensional Gaussian Orthogonal Ensemble, the probability of positive definiteness is given by: \begin{equation}p_n=\frac{{\mathrm{ Pf}}(A)}{2^{n(n+3)/4}\prod_{m=1}^n\Gamma(\frac m 2)},\label{theorem3} \end{equation} where $A$ is the $n'\times n'$ skew-symmetric matrix whose $(i,j)$-entry $a_{ij}$ is given for $i<j$ by: \begin{align}\label{pfaffianformula} a_{ij} = \begin{cases} \,2^{i+j-2}\Gamma(\frac{i+j}{2})\left(\beta_{\frac{1}{2}}(\frac{i}{2},\frac{j}{2})- \beta_{\frac{1}{2}}(\frac{j}{2},\frac{i}{2})\right) & \mbox{if }~ i<j\leq n, \\[.075in] \,2^{i-1}\Gamma(\frac{i}{2}) & \mbox{if }~ i<j=n+1.\\[.075in] \end{cases} \end{align} $($Note that the second case arises only when $n$ is odd.$)$

Here $\mathrm{Pf}$ is the Pfaffian, $\Gamma$ denotes the usual gamma function $\Gamma(s)=\int_{0}^\infty x^{s-1}e^{-x} \mathrm{d} x$, and $\beta_{t}$ is the usual incomplete beta function $\beta_{t}(i,j) = \int_0^{t} x^{i-1}(1-x)^{j-1} \mathrm{d} x.$ This formula follows from de Bruijn's identity for evaluating determinantal integrals.

As for the asymptotics, the leading term of $\log p_n$ is indeed $-\frac{\log(3)}{4} n^2$ as noted in Mikael's answer. If you are interested in the leading order asymptotic for $p_n$ it can be found from the results of Borot, Eynard, Majumdar and Nadal. Specialising to the GOE they give: \begin{align} p_n \sim \exp\left(-\frac{\log3}{4} n^2 - \frac{1}{2}\log\left(1+\frac{2}{\sqrt{3}}\right)n - \frac{1}{24}\log n + \alpha_1 \right)(1+o(1)).\label{asympexp} \end{align} The constant $\alpha_1$ is given by: \begin{align} \alpha_1= \frac{\zeta'(-1)}{2} - \frac{\log(2)}{12} -\frac{\log(3)}{16}+ \frac{1}{4}\log\left(1+\frac{2}{\sqrt{3}}\right)\simeq-0.0172, \end{align} where $\zeta(z)$ is the Riemann zeta function.

Edit: The exact formula lends itself easily to numerics. I have added a plot of $\log(p_n)$ for $1\leq n \leq 400.$ Fitting the data for $201\leq n \leq 400$ to $a n^2 + b n + c\log n +d$ we find $$−0.27465n^2 - 0.38382n -0.041504 \log n -0.018061.$$ A numerical approximation for the asymptotic formula is $$−0.27465n^2 - 0.38382n -0.041667 \log n -0.017223,$$ so we see good agreement. 