Distribution of the Determinant of a Complex-Valued Sample Correlation Matrix

In this post, we look at the distribution of the determinant of the sample correlation matrix of the realizations of a complex-valued Gaussian random vector. The distribution for real-valued Gaussian random vector was developed in [1], and we largely follow the developed framework. Thanks to Prashant Karunakaran for bringing this problem and the many applications to my attention in late 2017/early 2018.

Let $\boldsymbol{x}$ be a Gaussian random vector of length $p$ with mean $\boldsymbol{\mu} \in \mathbb{C}^p$ and covariance $\boldsymbol{\Sigma} \in \mathbb{C}^{p\times p}.$ Let $\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \dots, \boldsymbol{x}^{(n)}$ denote $n$ realizations, $n \geq p,$ of $\boldsymbol{x}.$ In the terminology of [1], the adjusted sample covariance matrix is given by $$\boldsymbol{S} = \frac{1}{n}\sum_{i = 1}^{n}(\boldsymbol{x}^{(i)} - \bar{\boldsymbol{x}})(\boldsymbol{x}^{(i)} - \bar{\boldsymbol{x}})^\mathrm{H},$$ where $\bar{\boldsymbol{x}}$ is the sample mean given by $$\bar{\boldsymbol{x}} = \frac{1}{n}\sum_{i = 1}^{n}\boldsymbol{x}^{(i)}.$$ Note that the adjusted sample covariance matrix is positive semi-definite.

The correlation matrix $\boldsymbol{R}$ is defined as: $$\boldsymbol{R} = \boldsymbol{D}^{-\frac{1}{2}} \boldsymbol{S} \boldsymbol{D}^{-\frac{1}{2}},$$ where $\boldsymbol{D} = \mathrm{Diag}(\boldsymbol{S})$ is a diagonal matrix with the diagonal elements of $\boldsymbol{S}$ on the main diagonal. Hence, $\boldsymbol{R}$ has unit diagonal elements and is independent of the variance of the elements of $\boldsymbol{x}.$ 

Now, for real-valued $\boldsymbol{x},$ the determinant of $\boldsymbol{R},$ denoted by $|\boldsymbol{R}|,$ is shown in [1, Theorem 2] to be a product of $p-1$ Beta-distributed scalar variables $\mathrm{Beta}(\frac{n-i}{2},\frac{i-1}{2}),$  $i=1,\dots,p-1.$ The density of the product can be given in terms of the $\mathrm{MeijerG}$ function [2] as follows [1, Theorem 2]:

$$g_\mathbb{R}(x;n,p) = \frac{\left[\Gamma(\frac{n-1}{2})\right]^{(p-1)} }{\Gamma(\frac{n-2}{2})\dots\Gamma(\frac{n-p}{2})} \mathrm{MeijerG}^{\begin{bmatrix}p-1 & 0 \\ p-1 & p-1\end{bmatrix}}\left(x\middle|\begin{matrix}\frac{n-3}{2},\dots,\frac{n-3}{2}\\ \frac{n-4}{2},\dots,\frac{n-(p+2)}{2}\end{matrix}\right).$$

Analogously, for complex-valued $\boldsymbol{x},$ $|\boldsymbol{R}|$ is a product of $p-1$ Beta-distributed scalar variables $\mathrm{Beta}(n-i,i),$  $i=1,\dots,p-1.$ The density of the product can now be given in terms of the $\mathrm{MeijerG}$ function, in a straightfoward manner, as follows.

$$g_\mathbb{C}(x;n,p) = \frac{\left[\Gamma(n-1)\right]^{(p-1)} }{\Gamma(n-1)\dots\Gamma(n-p+1)} \mathrm{MeijerG}^{\begin{bmatrix}p-1 & 0 \\ p-1 & p-1\end{bmatrix}}\left(x\middle|\begin{matrix}n-1,\dots,n-1\\ n-2,\dots,n-p\end{matrix}\right).$$

In the following, a Mathematica program for numerical simulation and the corresponding output are provided.

gC[x_, n_,   p_] := (Gamma[n])^(p - 1) /    Product[Gamma[n - i], {i, 1, p - 1}] MeijerG[{{},     Table[n - 1, {i, 1, p - 1}]}, {Table[n - i, {i, 2, p}], {}}, x]   r[x_] := Module[{d}, d = DiagonalMatrix[Diagonal[x]];    MatrixPower[d, -1/2] . x . MatrixPower[d, -1/2]] \[ScriptCapitalD] =   MatrixPropertyDistribution[    Det[r[(xr + I xi) .       ConjugateTranspose[xr + I xi]]], {xr \[Distributed]      MatrixNormalDistribution[IdentityMatrix[p], IdentityMatrix[n]],     xi \[Distributed]      MatrixNormalDistribution[IdentityMatrix[p],       IdentityMatrix[n]]}] ;   data = Re[RandomVariate[\[ScriptCapitalD], 100000]] ; \[ScriptCapitalD]1 = SmoothKernelDistribution[data] ;   Plot[{PDF[\[ScriptCapitalD]1, x], gC[x, n, p]}, {x, 0, 1},  PlotLabels -> {"Numerical", "Analytical"},  AxesLabel -> {"u", "p(u)"}]

The following figure shows that the numerical and analytical results match perfectly for the example case $n=4, k=6.$

[1] Pham-Gia, T. and Choulakian, V. (2014) Distribution of the Sample Correlation Matrix and      Applications. Open Journal of Statistics, 4, 330-344. doi: 10.4236/ojs.2014.45033.

[2] Weisstein, Eric W. "Meijer G-Function." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/MeijerG-Function.html

Distribution of A Simple Function of Gamma Variables

In this post, we look at a nifty result presented in [Krishnamoorthy2019] where the probability density function (pdf) of the function $\dfrac{r_1}{c+r_2}$ two independent Gamma distributed random variables $r_1 \sim \mathcal{G}(k_1,\theta_1)$ and $r_2 \sim \mathcal{G}(k_2, \theta_2)$ is derived.

The derivation is an exercise (for the scalar case) in computing the pdf of functions of variables by constructing the joint pdf and marginalizing based on the Jacobian. A similar approach can also be used for matrix-variate distributions (which will probably be a good topic for another post.)

Theorem

Let $c > 0,$ and let $r_1 \sim \mathcal{G}(k_1,\theta_1)$ and $r_2 \sim \mathcal{G}(k_2, \theta_2)$ be independent random variables, then, the pdf of $r = \dfrac{r_1}{c+r_2}$, denoted by $p_r(r;k_1,\theta_1,k_2,\theta_2,c),$ is given by

$\begin{equation}K_\mathrm{r} r^{k_1-1} \exp\left(-\frac{rc}{\theta_1}\right) U\left(k_2,k_1+k_2+1;c\left(\frac{r}{\theta_1}+\frac{1}{\theta_2}\right)\right),\end{equation}$ 

where $$U(a,b;z) = \frac{1}{\Gamma(a)} \int_{0}^{+\infty} \exp(-zx) x^{a-1} (1+x)^{b-a-1} \mathrm{d}x$$ is the hypergeometric U-function [Olver2010, Chapter 13, Kummer function], and $K_\mathrm{r}$ is a constant ensuring that the integral over the pdf equals one.

Proof

As $r_1$ and $r_2$ are independent, their joint pdf, denoted by $p_{r_1,r_2}(r_1,r_2;k_1,\theta_1,k_2,\theta_2),$ is given by

$\begin{equation}\frac{1}{\Gamma(k_1) \theta_1^{k_1} \Gamma(k_2) \theta_2^{k_2}} r_1^{k_1-1} r_2^{k_2-1} \exp\left(-\frac{r_1}{\theta_1}-\frac{r_2}{\theta_2}\right).\end{equation}$

Applying transformation $r = \frac{r_1}{c+r_2}$ with Jacobian $\frac{\mathrm{d} r_1}{\mathrm{d} r} = c+r_2,$ we obtain the transformed pdf, $p_{r,r_2}(r,r_2;k_1,\theta_1,k_2,\theta_2,c)$, as

$\begin{align}K_\mathrm{r}' r^{k_1-1} (c+r_2)^{k_1} r_2^{k_2-1}\exp\left(-\frac{rc}{\theta_1}-\left(\frac{r}{\theta_1} +\frac{1}{\theta_2}\right)r_2\right),\end{align}$

where $K_\mathrm{r}'$ is a constant ensuring that the integral over the pdf equals one. Next, 

$p_r(r;k_1,\theta_1,k_2,\theta_2,c)$ is obtained by marginalization as

$\begin{equation}p_r(r;k_1,\theta_1,k_2,\theta_2,c) = \int_{0}^{+\infty} \hspace{-0.25cm} p_{r,r_2}(r,r_2;k_1,\theta_1,k_2,\theta_2,c) \mathrm{d} r_2,\end{equation}$

where the integration is conducted using the integral representation of the hypergeometric U-function from [Olver2010, Chapter 13, Kummer function] to obtain the expression in the theorem statement.

References

[Krishnamoorthy2019] A. Krishnamoorthy and R. Schober, “Precoder design for two-user uplink MIMO-NOMA with simultaneous triangularization,” in Proc. IEEE Global Commun. Conf., Dec. 2019, pp. 1–6. https://doi.org/10.1109/GLOBECOM38437.2019.9014161

[Olver2010] F. W. Olver, D. Lozier, R. Boisvert and C. Clark, "NIST digital library of mathematical functions", Release, vol. 1, pp. 14, 2010.