# Confidence Intervals for Wishart Random Matrices

In this post we briefly describe Wishart-distributed random matrices and a result from (Chiani, 2017), which provides an algorithm for calculating the probability that a standard Wishart-distributed matrix's eigenvalues lie within a given interval. The functions described in this post are implemented in Python here, including Algorithm 1 from (Chiani, 2017). We also talk about quantile functions and non-standard Wishart-distributed matrices.

## The Wishart distribution

Suppose we have the $p\times p$ random matrix

$\begin{equation*} \bm{A} = \sum_{i=1}^n \bm{x}_i\bm{x}_i^T, \end{equation*}$where $n\geq p$ and each $\bm{x}_i\in\mathbb{R}^p$ is independently drawn from a zero-mean multivariate normal distribution

$\begin{equation*} \bm{x}_i \sim \mathcal{N}_p(\bm{0},\bm{\Sigma}) \end{equation*}$with positive definite covariance matrix $\bm{\Sigma}$. Then $\bm{A}$ is Wishart-distributed, and we write

$\begin{equation*} \bm{A} \sim \mathcal{W}_p(\bm{\Sigma},n), \end{equation*}$where $n$ is called the *degrees of freedom*. We refer to the special case
$\mathcal{W}_p(\bm{I}_p,n)$ as the *standard* Wishart distribution, where
$\bm{I}_p$ is the $p\times p$ identity matrix.

The Wishart distribution arises as the conjugate prior to the inverse covariance matrix of a multivariate normal distribution in Bayesian statistics, among other places.

One property of the Wishart distribution that we will make use of later is that if $\bm{A}\sim\mathcal{W}_p(\bm{\Sigma},n)$ and $\bm{C}\in\mathbb{R}^{p \times p}$ is a constant, full-rank matrix, then

$\begin{equation}\label{1} \bm{C}\bm{A}\bm{C}^T \sim \mathcal{W}_p(\bm{C}\bm{\Sigma}\bm{C}^T,n). \end{equation}$## Confidence bounds on the eigenvalues of standard Wishart matrices

Suppose we want to compute the probability that the eigenvalues of some standard Wishart-distributed matrix

$\begin{equation*} \bm{S} \sim \mathcal{W}_p(\bm{I}_p,n) \end{equation*}$are contained in a given interval $[a,b]$, with $0\leq a\leq b$. We denote this probability as

$\begin{equation}\label{2} \mathrm{Pr}[a\leq\lambda_{\min}(\bm{S}), \lambda_{\max}(\bm{S})\leq b], \end{equation}$where $\lambda_{\min}(\bm{S})$ and $\lambda_{\max}(\bm{S})$ are the minimum and maximum eigenvalues of $\bm{S}$, respectively. The probability $\eqref{2}$ is equivalent to

$\begin{equation}\label{3} \mathrm{Pr}[a\bm{I}_p \preccurlyeq \bm{S} \preccurlyeq b\bm{I}_p] \end{equation}$(see Appendix A.5.2 of (Boyd & Vandenberghe, 2004)), where the notation $\bm{A}\preccurlyeq\bm{B}$ means that $\bm{B}-\bm{A}$ is positive semidefinite.

Algorithm 1 of (Chiani, 2017) tells us how to compute the function $\psi(a,b)$ such that

$\begin{equation*} \mathrm{Pr}[a\bm{I}_p \preccurlyeq \bm{S} \preccurlyeq b\bm{I}_p] = \psi(a, b), \end{equation*}$and I've implemented this algorithm in Python here. This immediately gives us the cumulative density functions (CDFs) of the minimum and maximum eigenvalues of $\bm{S}$:

$\begin{align*} \mathrm{Pr}[a\leq\lambda_{\min}(\bm{S})] &= \psi(a, \infty) = 1 - C_{\min}(a), \\ \mathrm{Pr}[\lambda_{\max}(\bm{S})\leq b] &= \psi(0, b) = C_{\max}(b), \end{align*}$where $C_{\min}$ and $C_{\max}$ are the CDFs for the minimum and maximum eigenvalues, respectively.

## Quantile functions

We now have a method for computing the CDFs of the eigenvalues, but we do not have the inverse CDFs (i.e., the quantile functions), which compute the bounds required to achieve a given probability. For example, suppose we want to find the value $b\geq0$ such that $\lambda_{\max}(\bm{S})\leq b$ with a given probability $\rho\in[0,1]$. One way to do this is by solving the optimization problem

$\begin{equation*} \mathrm{argmin}_{b\geq0} (\rho - C_{\max}(b))^2, \end{equation*}$which is also implemented here, using one of scipy's built-in solvers.

## Non-standard Wishart matrices

We can extend the results above to the non-standard Wishart-distributed matrix

$\begin{equation*} \bm{A} \sim \mathcal{W}_p(\bm{\Sigma},n). \end{equation*}$In particular, generalizing $\eqref{3}$ we get

$\begin{equation}\label{4} \mathrm{Pr}[a\bm{\Sigma} \preccurlyeq \bm{A} \preccurlyeq b\bm{\Sigma}] = \psi(a,b). \end{equation}$To see this, let $\bm{L}$ be the Cholesky decomposition of $\bm{\Sigma}$, such that $\bm{L}\bm{L}^T=\bm{\Sigma}$. Since $\bm{\Sigma}$ is positive definite, $\bm{L}$ is full-rank and invertible. Then from $\eqref{1}$ we know that $\bm{L}^{-1}\bm{A}\bm{L}^{-T}\sim\mathcal{W}_p(\bm{L}^{-1}\bm{\Sigma}\bm{L}^{-T},n)=\mathcal{W}_p(\bm{I}_p,n)$, which means that $\mathrm{Pr}[a\bm{I}_p\preccurlyeq\bm{L}^{-1}\bm{A}\bm{L}^{-T}\preccurlyeq b\bm{I}_p]=\psi(a,b)$. Finally, since $a\bm{I}_p\preccurlyeq\bm{L}^{-1}\bm{A}\bm{L}^{-T}\preccurlyeq b\bm{I}_p$ if and only if $a\bm{\Sigma}\preccurlyeq\bm{A}\preccurlyeq b\bm{\Sigma}$ (this follows from Observation 7.1.8 of (Horn & Johnson, 2013)), we obtain $\eqref{4}$.

Notice that $\eqref{4}$ is no longer expressed in terms of eigenvalues of $\bm{A}$, but only expresses the probability that $\bm{A}$ is bounded by two matrices. This could be useful for formulating chance constraints on a positive definite matrix in convex optimization.

*Thanks to Abhishek Goudar for reading a draft of this.*

## References

- Boyd, Stephen and Lieven Vandenberghe (2004).
*Convex Optimization*. Cambridge University Press. (Freely available here.) - Chiani, M. (2017). "On the Probability That All Eigenvalues of
Gaussian, Wishart, and Double Wishart Random Matrices Lie Within an
Interval," in
*IEEE Transactions on Information Theory*, vol. 63, no. 7, pp. 4521–4531, July 2017, doi: 10.1109/TIT.2017.2694846, arxiv: 1502.04189. - Horn, Roger A. and Charles R. Johnson (2013).
*Matrix Analysis*. Cambridge University Press.