Wishart distribution

Occurrence

The Wishart distribution arises as the distribution of the sample covariance matrix for a sample from a multivariate normal distribution. It occurs frequently in likelihood-ratio tests in multivariate statistical analysis. It also arises in the spectral theory of random matrices and in multidimensional Bayesian analysis. It is also encountered in wireless communications, while analyzing the performance of Rayleigh fading MIMO wireless channels.

Probability density function

The Wishart distribution can be characterized by its probability density function as follows:

Let X be a p × p symmetric matrix of random variables that is positive semi-definite. Let V be a (fixed) symmetric positive definite matrix of size p × p.

Then, if n ≥ p, X has a Wishart distribution with n degrees of freedom if it has the probability density function

${\displaystyle f_{\mathbf{X}}(\mathbf{X})=\frac{1}{2^{np/2}{\left|\mathbf{V}\right|}^{n/2}{\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)}{\left|\mathbf{X}\right|}^{(n-p-1)/2}{e}^{-\frac{1}{2}\mathrm{tr}({\mathbf{V}}^{-1}\mathbf{X})}}$

where ${\displaystyle \left|\mathbf{X}\right|}$ is the determinant of ${\displaystyle \mathbf{X}}$ and Γ_p is the multivariate gamma function defined as

${\displaystyle {\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)={\pi }^{p(p-1)/4}\prod _{j=1}^p\mathrm{\Gamma }\left(\frac{n}{2}-\frac{j-1}{2}\right).}$

The density above is not the joint density of all the ${\displaystyle p^2}$ elements of the random matrix X (such ${\displaystyle p^2}$-dimensional density does not exist because of the symmetry constrains ${\displaystyle X_{ij}=X_{ji}}$), it is rather the joint density of ${\displaystyle p(p+1)/2}$ elements ${\displaystyle X_{ij}}$ for ${\displaystyle i\le j}$ (, page 38). Also, the density formula above applies only to positive definite matrices ${\displaystyle \mathbf{x};}$ for other matrices the density is equal to zero.

In fact the above definition can be extended to any real n > p − 1. If n ≤ p − 1, then the Wishart no longer has a density—instead it represents a singular distribution that takes values in a lower-dimension subspace of the space of p × p matrices.

Spectral density

The joint-eigenvalue density for the eigenvalues ${\displaystyle {\lambda }_1,\dots ,{\lambda }_p\ge 0}$ of a random matrix ${\displaystyle \mathbf{X}\sim W_p(\mathbf{I},n)}$ is,

${\displaystyle c_{n,p}{e}^{-\frac{1}{2}\sum _i{\lambda }_i}\prod {\lambda }_i^{(n-p-1)/2}\prod _{i<j}|{\lambda }_i-{\lambda }_j|}$

where ${\displaystyle c_{n,p}}$ is a constant. The spectral density can be marginalized to yield the density of a single eigenvalue, by evaluating a Selberg integral.

The spectral density can be integrated to give the probability that all eigenvalues of a Wishart random matrix lie within an interval.

Use in Bayesian statistics

In Bayesian statistics, in the context of the multivariate normal distribution, the Wishart distribution is the conjugate prior to the precision matrix Ω = Σ⁻¹, where Σ is the covariance matrix. The use of Normal-Wishart conjugate priors (for mean and precision) is particularly common for vector autoregression models.

Choice of parameters

The least informative, proper Wishart prior is obtained by setting n = p.

A common choice for V leverages the fact that the mean of X ~W_p(V, n) is nV. Then V is chosen so that nV equals an initial guess for X. For instance, when estimating a precision matrix Σ⁻¹ ~ W_p(V, n) a reasonable choice for V would be n⁻¹Σ₀⁻¹, where Σ₀ is some prior estimate for the covariance matrix Σ.

Properties

Log-expectation

The following formula plays a role in variational Bayes derivations for Bayes networks involving the Wishart distribution. From equation (2.63),

${\displaystyle \mathrm{E}[\ln \left|\mathbf{X}\right|]={\psi }_p\left(\frac{n}{2}\right)+p\ln (2)+\ln |\mathbf{V}|}$

where ${\displaystyle {\psi }_p}$ is the multivariate digamma function (the derivative of the log of the multivariate gamma function).

Log-variance

The following variance computation could be of help in Bayesian statistics:

${\displaystyle \mathrm{Var}\left[\ln \left|\mathbf{X}\right|\right]=\sum _{i=1}^p{\psi }_1\left(\frac{n+1-i}{2}\right)}$

where ${\displaystyle {\psi }_1}$ is the trigamma function. This comes up when computing the Fisher information of the Wishart random variable.

Entropy

The information entropy of the distribution has the following formula:

${\displaystyle \mathrm{H}\left[\mathbf{X}\right]=-\ln \left(B(\mathbf{V},n)\right)-\frac{n-p-1}{2}\mathrm{E}\left[\ln \left|\mathbf{X}\right|\right]+\frac{np}{2}}$

where B(V, n) is the normalizing constant of the distribution:

${\displaystyle B(\mathbf{V},n)=\frac{1}{{\left|\mathbf{V}\right|}^{n/2}{2}^{np/2}{\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)}.}$

This can be expanded as follows:

${\displaystyle \begin{array}{cc}\mathrm{H}\left[\mathbf{X}\right]& =\frac{n}{2}\ln \left|\mathbf{V}\right|+\frac{np}{2}\ln 2+\ln {\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)-\frac{n-p-1}{2}\mathrm{E}\left[\ln \left|\mathbf{X}\right|\right]+\frac{np}{2}\\ & =\frac{n}{2}\ln \left|\mathbf{V}\right|+\frac{np}{2}\ln 2+\ln {\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)-\frac{n-p-1}{2}\left({\psi }_p\left(\frac{n}{2}\right)+p\ln 2+\ln \left|\mathbf{V}\right|\right)+\frac{np}{2}\\ & =\frac{n}{2}\ln \left|\mathbf{V}\right|+\frac{np}{2}\ln 2+\ln {\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)-\frac{n-p-1}{2}{\psi }_p\left(\frac{n}{2}\right)-\frac{n-p-1}{2}\left(p\ln 2+\ln \left|\mathbf{V}\right|\right)+\frac{np}{2}\\ & =\frac{p+1}{2}\ln \left|\mathbf{V}\right|+\frac{1}{2}p(p+1)\ln 2+\ln {\mathrm{\Gamma }}_p\left(\frac{n}{2}\right)-\frac{n-p-1}{2}{\psi }_p\left(\frac{n}{2}\right)+\frac{np}{2}\end{array}}$

Cross-entropy

The cross-entropy of two Wishart distributions ${\displaystyle p_0}$ with parameters ${\displaystyle n_0,V_0}$ and ${\displaystyle p_1}$ with parameters ${\displaystyle n_1,V_1}$ is

${\displaystyle \begin{array}{cc}H(p_0,p_1)& ={\mathrm{E}}_{p_0}[-\log p_1]\\ & ={\mathrm{E}}_{p_0}\left[-\log \frac{{\left|\mathbf{X}\right|}^{(n_1-p_1-1)/2}{e}^{-\mathrm{tr}({\mathbf{V}}_1^{-1}\mathbf{X})/2}}{2^{n_1{p}_1/2}{\left|{\mathbf{V}}_1\right|}^{n_1/2}{\mathrm{\Gamma }}_{p_1}\left(\frac{n_1}{2}\right)}\right]\\ & =\frac{n_1{p}_1}{2}\log 2+\frac{n_1}{2}\log \left|{\mathbf{V}}_1\right|+\log {\mathrm{\Gamma }}_{p_1}(\frac{n_1}{2})-\frac{n_1-p_1-1}{2}{\mathrm{E}}_{p_0}\left[\log \left|\mathbf{X}\right|\right]+\frac{1}{2}{\mathrm{E}}_{p_0}\left[\mathrm{tr}\left({\mathbf{V}}_1^{-1}\mathbf{X}\right)\right]\\ & =\frac{n_1{p}_1}{2}\log 2+\frac{n_1}{2}\log \left|{\mathbf{V}}_1\right|+\log {\mathrm{\Gamma }}_{p_1}(\frac{n_1}{2})-\frac{n_1-p_1-1}{2}\left({\psi }_{p_0}(\frac{n_0}{2})+p_0\log 2+\log \left|{\mathbf{V}}_0\right|\right)+\frac{1}{2}\mathrm{tr}\left({\mathbf{V}}_1^{-1}{n}_0{\mathbf{V}}_0\right)\\ & =-\frac{n_1}{2}\log \left|{\mathbf{V}}_1^{-1}{\mathbf{V}}_0\right|+\frac{p_1+1}{2}\log \left|{\mathbf{V}}_0\right|+\frac{n_0}{2}\mathrm{tr}\left({\mathbf{V}}_1^{-1}{\mathbf{V}}_0\right)+\log {\mathrm{\Gamma }}_{p_1}\left(\frac{n_1}{2}\right)-\frac{n_1-p_1-1}{2}{\psi }_{p_0}(\frac{n_0}{2})+\frac{n_1(p_1-p_0)+p_0(p_1+1)}{2}\log 2\end{array}}$

Note that when ${\displaystyle p_0=p_1}$ and ${\displaystyle n_0=n_1}$ we recover the entropy.

KL-divergence

The Kullback–Leibler divergence of ${\displaystyle p_1}$ from ${\displaystyle p_0}$ is

${\displaystyle \begin{array}{cc}{D}_{KL}(p_0\Vert p_1)& =H(p_0,p_1)-H(p_0)\\ & =-\frac{n_1}{2}\log |{\mathbf{V}}_1^{-1}{\mathbf{V}}_0|+\frac{n_0}{2}(\mathrm{tr}({\mathbf{V}}_1^{-1}{\mathbf{V}}_0)-p)+\log \frac{{\mathrm{\Gamma }}_p\left(\frac{n_1}{2}\right)}{{\mathrm{\Gamma }}_p\left(\frac{n_0}{2}\right)}+\frac{n_0-n_1}{2}{\psi }_p\left(\frac{n_0}{2}\right)\end{array}}$

Characteristic function

The characteristic function of the Wishart distribution is

${\displaystyle \mathrm{\Theta }\mapsto \mathrm{E}\left[\exp \left(i\mathrm{tr}\left(\mathbf{X}\mathbf{\Theta }\right)\right)\right]={\left|1-2i\mathbf{\Theta }\mathbf{V}\right|}^{-n/2}}$

where E[⋅] denotes expectation. (Here Θ is any matrix with the same dimensions as V, 1 indicates the identity matrix, and i is a square root of −1). Properly interpreting this formula requires a little care, because noninteger complex powers are multivalued; when n is noninteger, the correct branch must be determined via analytic continuation.

Theorem

If a p × p random matrix X has a Wishart distribution with m degrees of freedom and variance matrix V — write ${\displaystyle \mathbf{X}\sim {\mathcal{W}}_p(\mathbf{V},m)}$ — and C is a q × p matrix of rank q, then

${\displaystyle \mathbf{C}\mathbf{X}{\mathbf{C}}^T\sim {\mathcal{W}}_q\left(\mathbf{C}\mathbf{V}{\mathbf{C}}^T,m\right).}$

Corollary 1

If z is a nonzero p × 1 constant vector, then:

${\displaystyle {\sigma }_z^{-2}{\mathbf{z}}^T\mathbf{X}\mathbf{z}\sim {\chi }_m^2.}$

In this case, ${\displaystyle {\chi }_m^2}$ is the chi-squared distribution and ${\displaystyle {\sigma }_z^2={\mathbf{z}}^T\mathbf{V}\mathbf{z}}$ (note that ${\displaystyle {\sigma }_z^2}$ is a constant; it is positive because V is positive definite).

Corollary 2

Consider the case where z^T = (0, ..., 0, 1, 0, ..., 0) (that is, the j-th element is one and all others zero). Then corollary 1 above shows that

${\displaystyle {\sigma }_{jj}^{-1}{w}_{jj}\sim {\chi }_m^2}$

gives the marginal distribution of each of the elements on the matrix's diagonal.

George Seber points out that the Wishart distribution is not called the “multivariate chi-squared distribution” because the marginal distribution of the off-diagonal elements is not chi-squared. Seber prefers to reserve the term multivariate for the case when all univariate marginals belong to the same family.

Estimator of the multivariate normal distribution

The Wishart distribution is the sampling distribution of the maximum-likelihood estimator (MLE) of the covariance matrix of a multivariate normal distribution. A derivation of the MLE uses the spectral theorem.

Bartlett decomposition

The Bartlett decomposition of a matrix X from a p-variate Wishart distribution with scale matrix V and n degrees of freedom is the factorization:

${\displaystyle \mathbf{X}=\mathbf{\text{L}}\mathbf{\text{A}}{\mathbf{\text{A}}}^T{\mathbf{\text{L}}}^T,}$

where L is the Cholesky factor of V, and:

${\displaystyle \mathbf{A}=\left(\begin{array}{ccccc}{c}_1& 0& 0& \cdots & 0\\ n_{21}& c_2& 0& \cdots & 0\\ n_{31}& n_{32}& c_3& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ n_{p1}& n_{p2}& n_{p3}& \cdots & c_p\end{array}\right)}$

where ${\displaystyle c_i^2\sim {\chi }_{n-i+1}^2}$ and n_ij ~ N(0, 1) independently. This provides a useful method for obtaining random samples from a Wishart distribution.

Marginal distribution of matrix elements

Let V be a 2 × 2 variance matrix characterized by correlation coefficient −1 < ρ < 1 and L its lower Cholesky factor:

${\displaystyle \mathbf{V}=\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right),\mathbf{L}=\left(\begin{array}{cc}{\sigma }_1& 0\\ \rho {\sigma }_2& \sqrt{1-{\rho }^2}{\sigma }_2\end{array}\right)}$

Multiplying through the Bartlett decomposition above, we find that a random sample from the 2 × 2 Wishart distribution is

${\displaystyle \mathbf{X}=\left(\begin{array}{cc}{\sigma }_1^2{c}_1^2& {\sigma }_1{\sigma }_2\left(\rho c_1^2+\sqrt{1-{\rho }^2}{c}_1{n}_{21}\right)\\ {\sigma }_1{\sigma }_2\left(\rho c_1^2+\sqrt{1-{\rho }^2}{c}_1{n}_{21}\right)& {\sigma }_2^2\left(\left(1-{\rho }^2\right)c_2^2+{\left(\sqrt{1-{\rho }^2}{n}_{21}+\rho c_1\right)}^2\right)\end{array}\right)}$

The diagonal elements, most evidently in the first element, follow the χ² distribution with n degrees of freedom (scaled by σ²) as expected. The off-diagonal element is less familiar but can be identified as a normal variance-mean mixture where the mixing density is a χ² distribution. The corresponding marginal probability density for the off-diagonal element is therefore the variance-gamma distribution

${\displaystyle f(x_{12})=\frac{{\left|x_{12}\right|}^{\frac{n-1}{2}}}{\mathrm{\Gamma }\left(\frac{n}{2}\right)\sqrt{2^{n-1}\pi \left(1-{\rho }^2\right){\left({\sigma }_1{\sigma }_2\right)}^{n+1}}}\cdot K_{\frac{n-1}{2}}\left(\frac{\left|x_{12}\right|}{{\sigma }_1{\sigma }_2\left(1-{\rho }^2\right)}\right)\exp \left(\frac{\rho x_{12}}{{\sigma }_1{\sigma }_2(1-{\rho }^2)}\right)}$

where K_ν(z) is the modified Bessel function of the second kind. Similar results may be found for higher dimensions. In general, if ${\displaystyle X}$ follows a Wishart distribution with parameters, ${\displaystyle \mathrm{\Sigma },n}$, then for ${\displaystyle i\ne j}$, the off-diagonal elements

${\displaystyle X_{ij}\sim \text{VG}(n,{\mathrm{\Sigma }}_{ij},({\mathrm{\Sigma }}_{ii}{\mathrm{\Sigma }}_{jj}-{\mathrm{\Sigma }}_{ij}^2{)}^{1/2},0)}$.

It is also possible to write down the moment-generating function even in the noncentral case (essentially the nth power of Craig (1936) equation 10) although the probability density becomes an infinite sum of Bessel functions.

The range of the shape parameter

It can be shown that the Wishart distribution can be defined if and only if the shape parameter n belongs to the set

${\displaystyle {\mathrm{\Lambda }}_p:=\{0,\dots ,p-1\}\cup \left(p-1,\mathrm{\infty }\right).}$

This set is named after Simon Gindikin, who introduced it in the 1970s in the context of gamma distributions on homogeneous cones. However, for the new parameters in the discrete spectrum of the Gindikin ensemble, namely,

${\displaystyle {\mathrm{\Lambda }}_p^{\ast }:=\{0,\dots ,p-1\},}$

the corresponding Wishart distribution has no Lebesgue density.

Wishart–Laguerre ensembles and β-extensions

In random matrix theory, the Wishart family is often studied through its Laguerre ensembles. For the real case (orthogonal symmetry, β=1), the joint density of the eigenvalues ${\displaystyle {\lambda }_1,\dots ,{\lambda }_p\ge 0}$ of ${\displaystyle \mathbf{X}\sim W_p(\mathbf{I},n)}$ is

${\displaystyle c_{n,p}\exp (-\frac{1}{2}\sum _{i=1}^p{\lambda }_i)\prod _{i=1}^p{\lambda }_i^{\frac{n-p-1}{2}}\prod _{1\le i<j\le p}|{\lambda }_i-{\lambda }_j|,}$

which is the Laguerre orthogonal ensemble (LOE). The complex and quaternion analogues are the Laguerre unitary (LUE, ${\displaystyle \beta =2}$) and Laguerre symplectic (LSE, ${\displaystyle \beta =4}$) ensembles, respectively.

β-Laguerre ensemble (general β > 0 {\displaystyle \beta >0} )

A further generalization, the β-Laguerre ensemble, allows the Dyson index ${\displaystyle \beta >0}$ to vary continuously. Its joint eigenvalue density has the Coulomb gas form

${\displaystyle p({\lambda }_1,\dots ,{\lambda }_p)\propto (\prod _{i=1}^p{\lambda }_i^{\alpha }{e}^{-\frac{\beta }{2}{\lambda }_i})(\prod _{i<j}|{\lambda }_i-{\lambda }_j{|}^{\beta }),\alpha >-1,}$

which reduces to LOE/LUE/LSE for ${\displaystyle \beta =1,2,4}$. For the classical Gaussian Wishart case (scale ${\displaystyle \mathbf{I}}$), one has the identification

${\displaystyle \alpha =\frac{\beta }{2}(n-p+1)-1.}$

A concrete probabilistic construction for any ${\displaystyle \beta >0}$ is provided by the Dumitriu–Edelman bidiagonal model. One samples a random bidiagonal matrix ${\displaystyle \mathbf{B}}$ with independent chi variables and sets ${\displaystyle \mathbf{L}=\mathbf{B}{\mathbf{B}}^T}$; the eigenvalues of ${\displaystyle \mathbf{L}}$ follow the β-Laguerre law with parameter ${\displaystyle \alpha }$.

Sampling

• Classical ${\displaystyle \beta =1,2,4}$. Bartlett's decomposition gives ${\displaystyle \mathbf{X}=\mathbf{L}\mathbf{A}{\mathbf{A}}^T{\mathbf{L}}^T}$ with ${\displaystyle \mathbf{A}}$ lower-triangular having chi/normal entries, yielding LOE/LUE/LSE exactly.
• General ${\displaystyle \beta >0}$. In the Dumitriu–Edelman model, for matrix size ${\displaystyle p}$ and parameter ${\displaystyle \alpha >-1}$, one samples independently

${\displaystyle B_{ii}\sim \frac{1}{\sqrt{\beta }}{\chi }_{\beta (\alpha +p-i+1)},B_{i,i-1}\sim \frac{1}{\sqrt{\beta }}{\chi }_{\beta (p-i+1)}(i=2,\dots ,p),}$

and sets ${\displaystyle \mathbf{L}=\mathbf{B}{\mathbf{B}}^T}$. Then the eigenvalues of ${\displaystyle \mathbf{L}}$ have the β-Laguerre joint density above.

Rectangular data matrices

If one wishes to realize these spectra as the singular values of a rectangular ${\displaystyle M\times N}$ matrix (with ${\displaystyle M\ge N=p}$), draw independent Haar-distributed ${\displaystyle \mathbf{U}\in O(M)}$ (or ${\displaystyle U(M)}$) and ${\displaystyle \mathbf{V}\in O(N)}$, and set

${\displaystyle X=U\left[\begin{array}{c}\mathrm{diag}(\sqrt{{\lambda }_1},\dots ,\sqrt{{\lambda }_p})\\ 0\end{array}\right]V^{\ast },}$

so that ${\displaystyle X^{\ast }X}$ has eigenvalues ${\displaystyle \{{\lambda }_i\}}$. For i.i.d. Gaussian entries (true Wishart), ${\displaystyle \beta }$ is fixed by the field (real, complex, quaternion) and ${\displaystyle \alpha }$ reduces to ${\displaystyle \frac{\beta }{2}(n-p+1)-1}$.

Hard-edge behavior and universality

At the “hard edge” (near ${\displaystyle \lambda =0}$), β-Laguerre ensembles exhibit Bessel-kernel correlations and level repulsion of order ${\displaystyle \beta }$. In particular, with ${\displaystyle p\to \mathrm{\infty }}$ and fixed ${\displaystyle \alpha }$, the distribution of the smallest eigenvalue converges to a universal hard-edge (Bessel) law that depends only on ${\displaystyle \alpha }$ and ${\displaystyle \beta }$.

Macroscopic limit (Marchenko–Pastur law)

Under proportional growth ${\displaystyle p/n\to \gamma \in (0,\mathrm{\infty })}$, the empirical spectral distribution of ${\displaystyle n^{-1}\mathbf{G}{\mathbf{G}}^T}$ (with i.i.d. entries of variance ${\displaystyle {\sigma }^2}$) converges almost surely to the Marchenko–Pastur law, supported on ${\displaystyle {\sigma }^2(1\pm \sqrt{\gamma }{)}^2}$. When ${\displaystyle \gamma =1}$, the density diverges as ${\displaystyle x^{-1/2}}$ at the origin (a hard-edge singularity).

Relationships to other distributions

• The Wishart distribution is related to the inverse-Wishart distribution, denoted by ${\displaystyle W_p^{-1}}$, as follows: If X ~ W_p(V, n) and if we do the change of variables C = X⁻¹, then ${\displaystyle \mathbf{C}\sim W_p^{-1}({\mathbf{V}}^{-1},n)}$. This relationship may be derived by noting that the absolute value of the Jacobian determinant of this change of variables is |C|^p+1, see for example equation (15.15) in.
• In Bayesian statistics, the Wishart distribution is a conjugate prior for the precision parameter of the multivariate normal distribution, when the mean parameter is known.
• A generalization is the multivariate gamma distribution.
• A different type of generalization is the normal-Wishart distribution, essentially the product of a multivariate normal distribution with a Wishart distribution.