Numerically Stable Evaluation of Moments of Random Gram Matrices with Applications

TL;DR

This paper introduces a numerically stable method for computing positive moments of one-side correlated random Gram matrices, improving accuracy in high-dimensional settings and enabling better eigenvalue distribution approximations.

Contribution

It provides a new stable closed-form approach for evaluating moments of correlated Gram matrices, enhancing numerical accuracy over existing methods.

Findings

The method achieves higher numerical stability in high-dimensional regimes.

It allows accurate approximation of eigenvalue distributions.

The approach improves the reliability of moment-based spectral analysis.

Abstract

This paper is focuses on the computation of the positive moments of one-side correlated random Gram matrices. Closed-form expressions for the moments can be obtained easily, but numerical evaluation thereof is prone to numerical stability, especially in high-dimensional settings. This letter provides a numerically stable method that efficiently computes the positive moments in closed-form. The developed expressions are more accurate and can lead to higher accuracy levels when fed to moment based-approaches. As an application, we show how the obtained moments can be used to approximate the marginal distribution of the eigenvalues of random Gram matrices.

Tables

Table 1. Table I : Positive moments evaluation of the Gram matrix 𝐖 𝐖 \mathbf{W} with 𝚲 𝚲 \mathbf{\Lambda} as in ( 14 ) with ξ = 0.85 𝜉 0.85 \xi=0.85 . Both settings are considered: config 1 : n t = 3 , q = 5 : subscript config 1 formulae-sequence subscript 𝑛 𝑡 3 𝑞 5 \text{config}_{1}:n_{t}=3,q=5 and config 2 : n t = 3 , q = 20 : subscript config 2 formulae-sequence subscript 𝑛 𝑡 3 𝑞 20 \text{config}_{2}:n_{t}=3,q=20 for different moment order values, p 𝑝 p .

	Formula in [3]	Empirical(${10}^6$ realizations)	Proposed
${\text{config}}_1$, $p=1$	0.5563	0.5563	0.5562
${\text{config}}_2$, $p=1$	-2.1177e+06	0.9613	0.9612
${\text{config}}_1$, $p=5$	1.1029	1.1031	1.1032
${\text{config}}_2$, $p=5$	7.7212e+03	4.7575	4.7562
${\text{config}}_1$, $p=8$	5.3799	5.4332	5.3989
${\text{config}}_2$, $p=8$	5.4568e+04	37.3178	37.47

Equations

$$f_{\lambda }\left(\lambda \right)$$

$$\times \sum _{l=1}^q\sum _{k=q-n_t+1}^q\frac{{\lambda }^{n_t+k-q-1}{e}^{-\lambda /{\beta }_l}{\beta }_l^{q-n_t-1}}{\Gamma \left(n_t-q+k\right)}{D}_{l,k},$$

$$\{\mathbf{\Psi }{\}}_{m,n}={\beta }_m^{n-1}.$$

$${\mathbf{\Psi }}^{-1}$$

$$=\frac{1}{{\prod }_{i<j}^q\left({\beta }_j-{\beta }_i\right)}{\mathbf{D}}^T,$$

$$f_{\lambda }\left(\lambda \right)=\frac{1}{n_t}\sum _{l=1}^q\sum _{k=q-n_t+1}^q\frac{{\lambda }^{n_t+k-q-1}{e}^{-\lambda /{\beta }_l}{\beta }_l^{q-n_t-1}}{\Gamma \left(n_t-q+k\right)}{\mathbf{\Psi }}_{k,l}^{-1}.$$

$$F_{\lambda }\left(\lambda \right)=\frac{1}{n_t}\sum _{l=1}^q\sum _{k=q-n_t+1}^q\frac{{\beta }_l^{k-1}\gamma \left(n_t-q+k,\lambda /{\beta }_l\right)}{\Gamma \left(n_t-q+k\right)}{\mathbf{\Psi }}_{k,l}^{-1},$$

$$\mathbb{E}\left[g(\mathbf{W})\right]$$

$$={\int }_0^{\infty }g(\lambda )f_{\lambda }(\lambda )d\lambda ,$$

$$\begin{array}{rl}{\mu }_{\mathbf{W}}\left(p\right)& =\mathbb{E}\left[{\lambda }^p\right]\\ & =\frac{1}{n_t}\mathbb{E}\mathrm{tr}\left[{\mathbf{W}}^p\right]\\ & ={\int }_0^{\infty }{\lambda }^p{f}_{\lambda }\left(\lambda \right)\\ & =\frac{1}{n_t}\sum _{k=q-n_t+1}^q\frac{\Gamma \left(n_t+p+k-q\right)}{\Gamma \left(n_t-q+k\right)}\sum _{l=1}^q{\mathbf{\Psi }}_{k,l}^{-1}{\beta }_l^{p+k-1}.\end{array}$$

$$\sum _{l=1}^q{\mathbf{\Psi }}_{k,l}^{-1}{\beta }_l^{p+k-1},$$

$${\alpha }_{k,\tau }=\sum _{l=1}^q{\mathbf{\Psi }}_{k,l}^{-1}{\beta }_l^{\tau }.$$

$$\mathbf{\Psi }{\mathbf{\alpha }}_{\tau }={\mathbf{\beta }}_{\tau }.$$

$${\beta }_k^{\tau }=\sum _{l=1}^q{\beta }_k^l{\alpha }_{l,\tau },k=1,\cdots ,q$$

$$P\left(X\right)=\sum _{k=1}^q{\alpha }_{k,\tau }{X}^{k-1}-X^{\tau }.$$

$$P\left(X\right)=Q\left(X\right)\prod _{i=1}^q\left(X-{\beta }_i\right),$$

$$\{\begin{array}{ll}& a_{q+1}{b}_1+a_q{b}_2+\cdots a_{2q+1-\tau }{b}_{\tau -q+1}=0\\ & a_{q+1}{b}_2+a_q{b}_3+\cdots a_{2q+2-\tau }{b}_{\tau -q+1}=0\\ & ⋮\\ & a_{q+1}{b}_{\tau -q+1}=-1.\end{array}$$

$$\mathbf{\Phi }\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_{\tau -q+1}\text{par}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ -1\end{array}\right],$$

$$\mathbf{\Phi }=\left[\begin{array}{ccccc}{a}_{q+1}& a_q& \cdots & & \cdots a_{2q+1-\tau }\\ 0& a_{q+1}& a_q& \cdots & \cdots a_{2q+2-\tau }\\ ⋮& \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & a_q\\ 0& \cdots & & 0& a_{q+1}\end{array}\right].$$

$$\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_{\tau -q+1}\text{par}\end{array}\right]=\mathbf{\Phi }\setminus \left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ -1\text{par}\end{array}\right].$$

$${\alpha }_{j,\tau }=\sum _{k=1}^{\tau -q+1}{b}_k{a}_{j+2-k}.$$

$$\mathbf{\Lambda }=\left(1-\xi \right)\mathrm{diag}\left(1,\xi ,{\xi }^2,\cdots ,{\xi }^{q-1}\right),0\le \xi \le 1,$$

$${\int }_0^{\infty }\frac{-\log f\left({\lambda }^2\right)}{1+{\lambda }^2}d\lambda =\infty$$

$$\underset{\lambda \to \infty }{\lim }\frac{-\lambda \frac{\partial f\left(\lambda \right)}{\partial \lambda }}{f\left(\lambda \right)}=\infty$$

$$f_{\lambda }\left(\lambda \right)=\frac{{\lambda }^{\nu }{e}^{-\lambda /c}}{c^{\nu +1}}\sum _{i=0}^{\infty }{\delta }_i{\mathcal{L}}_i\left(\nu ,\lambda /c\right),$$

$${\mathcal{L}}_i\left(\nu ,\lambda \right)=\sum _{k=0}^i{\left(-1\right)}^k\frac{\Gamma \left(\nu +i+1\right){\lambda }^{i-k}}{k!\left(i-k\right)!\Gamma \left(\nu +i-k+1\right)},$$

$${\delta }_i=\sum _{k=0}^i\frac{{\left(-1\right)}^k}{c^{i-k}}\frac{i!}{k!\left(i-k\right)!\Gamma \left(\nu +i-k+1\right)}{\mu }_{\mathbf{W}}\left(i-k\right).$$

$$f_{\lambda ,K}\left(\lambda \right)=\frac{{\lambda }^{\nu }{e}^{-\lambda /c}}{c^{\nu +1}}\sum _{i=0}^K{\delta }_i{\mathcal{L}}_i\left(\nu ,\lambda /c\right).$$

$$\begin{array}{rl}{F}_{\lambda ,K}\left(\lambda \right)=\sum _{i=0}^K{\delta }_i\sum _{k=0}^i& {\left(-1\right)}^k\frac{\gamma \left(i+\nu -k+1,\lambda /c\right)}{k!\left(i-k\right)!}.\end{array}$$

$$f_{\lambda }(\lambda )=\sum _{l=1}^q\sum _{k=q-n_t+1}^q{\gamma }_{k,l}{\lambda }^{n_t+k-q-1}{e}^{-\lambda /{\beta }_l},$$

$$\begin{array}{rl}{f}_{\lambda }({\lambda }^2)& =\sum _{l=1}^q\sum _{k=q-s+1}^q{\gamma }_{k,l}{\lambda }^{2(n_t+k-q-1)}{e}^{-{\lambda }^2/{\beta }_l}\\ & \le e^{-{\lambda }^2/{\beta }_q}\sum _{l=1}^q\sum _{k=q-s+1}^q{\gamma }_{k,l}{\lambda }^{2(n_t+k-q-1)}.\end{array}$$

$$-\log \left(f_{\lambda }({\lambda }^2)\right)\ge {\lambda }^2/{\beta }_q-\log \left(\sum _{l=1}^q\sum _{k=q-n_t+1}^q{\gamma }_{k,l}{\lambda }^{2(n_t+k-q-1)}\right).$$

Full text

I Introduction

Gram random matrices with one-sided correlation naturally arise in the context of signal processing [1] and wireless communications [2, 3]. For instance, in signal processing, inverse moments of this kind of matrices are used to evaluate the performance of linear estimators such as the best linear unbiased estimator (BLUE) and optimize the design of some covariance matrix estimators [1, 4]. In wireless communications, Gram random matrices arise as a key element in the computation of the ergodic capacity of amplify and forward (AF) multiple input multiple output (MIMO) dual-hop systems [3].

Given a random matrix $\mathbf{H}\in\mathbb{C}^{q\times n_{t}}$ such that $\mathbf{H}=\mathbf{\Lambda}^{\frac{1}{2}}\mathbf{X}$, where $\mathbf{X}$ is a standard complex Gaussian random matrix and $\mathbf{\Lambda}$ is a Hermitian positive definite matrix, the authors in [3] derive a closed form expression of the marginal probability density function (PDF) of an unordered eigenvalue of the Gram matrix $\mathbf{W}=\mathbf{H}^{H}\mathbf{H}$ for arbitrary dimensions $n_{t}$ and $q$. The result, although being very useful, involves the inversion of a large Vandermonde matrix and as such might not be always numerically stable, especially in the following situations:

• •

The dimensions $q$ belongs to moderate to large values.

• •

The gap between the eigenvalues of $\mathbf{\Lambda}$ is small.

To solve this problem, an expression of the exact marginal PDF has been proposed in [5] when some eigenvalues are identical. However, the problem of numerical stability remains when some eigenvalues of $\boldsymbol{\Lambda}$ are different but close to each other. Motivated by these facts, we provide a more stable method to compute the positive moments of $\mathbf{W}$ without the need to invert large Vandermonde matrices. The contributions of this letter are summarized as follows:

• •

We provide a numerically stable method to evaluate the positive moments of $\mathbf{W}$.

• •

Using Laguerre polynomials and based on the calculated positive moments, we provide a numerically stable approximation of the marginal probability density function (PDF) of the eigenvalues of $\mathbf{W}$.

The remainder of this letter is organized as follows. In Section II, we provide the main steps to efficiently compute the positive moments of $\mathbf{W}$. In Section III, we propose to approximate the marginal PDF based on the computed positive moments and using Laguerre polynomials. In Section IV, we present some numerical results to validate our method and finally we conclude our work in Section V.

II A numerically stable method to compute the moments of random Gram matrices

II-A Problem statement

Let $\mathbf{H}=\boldsymbol{\Lambda}^{\frac{1}{2}}{\bf X}\in\mathbb{C}^{q\times n_{t}}$ where $\mathbf{\Lambda}\in\mathbb{C}^{q\times q}$ is a positive definite matrix with distinct eigenvalues $0<\beta_{1}<\beta_{2}<\cdots<\beta_{q}$ and ${\bf X}$ a standard complex Gaussian matrix. Assume that $n_{t}\leq q$. The marginal PDF of an unordered eigenvalue $\lambda$ of ${\bf W}={\bf H}^{H}{\bf H}$ is given by [3, Lemma 1]

[TABLE]

where $D_{l,k}=\{\mathbf{D}\}_{l,k}$ is the $\left(l,k\right)th$ cofactor of the Vandermonde $q\times q$ matrix $\mathbf{\Psi}$ whose $\left(m,n\right)$th entry is

[TABLE]

Expressing the inverse of $\boldsymbol{\Psi}$ as

[TABLE]

the PDF in (II-A) simplifies to

[TABLE]

The cumulative density function (CDF) can thus be easily derived as

[TABLE]

where $\Gamma\left(.\right)$ and $\gamma\left(.,.\right)$ are respectively the standard Gamma and the lower incomplete Gamma functions.

Knowing the marginal PDF, it is possible to compute the expected value of any functional $g$ of the eigenvalues of ${\bf W}$. Indeed, we have

[TABLE]

where $\left\{\lambda_{i}({\bf W})\right\}_{i=1}^{n_{t}}$ are the eigenvalues of ${\bf W}$. Equation (5) is very useful in practice as it can be leveraged to compute performance metrics of many wireless communication and signal processing schemes. Examples include the ergodic capacity, the SINR at the output of the MMSE receiver and the MSE of the BLUE estimator which correspond respectively to selecting $g(x)$ as $g(x)=\log_{2}(x+\sigma^{2})$, $g(x)=\frac{1}{x+\sigma^{2}}$ where $\sigma^{2}$ is the noise variance and $g(x)=x^{-1}$.

When it comes to numerically compute $\mathbb{E}\left[g({\bf W})\right]$, it is easy to see that, when the eigenvalues $\left\{\beta_{i}\right\}$ are very close causing the matrix $\boldsymbol{\Psi}$ to be ill-conditioned, some numerical stability issues might occur. In this work, we show that for some functionals $g$, namely polynomials, it is possible to evaluate $\mathbb{E}\left[g({\bf W})\right]$ in a stable way. This allows us, using moment approximation techniques, to obtain a numerically stable approximation of the marginal PDF. We believe that the same approximation method can also be extended to approximate $\mathbb{E}\left[g({\bf W})\right]$ for any functional $g$ of interest.

II-B A Numerically stable method to compute positive moments

In this section, we propose a numerically stable technique to compute the positive moments of ${\bf W}$. Let $p\in\mathbb{N}$, the $p$-th moment of matrix ${\bf W}$ is given by

[TABLE]

In many practical scenarios, numerical instability might originate from the computation of the following quantity

[TABLE]

as $\mathbf{\Psi}$ is ill-conditioned. To overcome this issue, we propose an alternative way that avoids computing the inverse of $\boldsymbol{\Psi}$. For $k\in\llbracket 1,q\rrbracket$ and $\tau\in\llbracket p+q-n_{t},p+q-1\rrbracket$, define $\alpha_{k,\tau}$ as

[TABLE]

The basic idea is based on the observation that $\boldsymbol{\alpha}_{\tau}\triangleq\left[\alpha_{1,\tau},\cdots,\alpha_{q,\tau}\right]^{T}$ is solution to the following linear system

[TABLE]

where $\boldsymbol{\beta}_{\tau}=\left[\beta_{1}^{\tau},\cdots,\beta_{q}^{\tau}\right]^{T}$. If $0\leq\tau\leq q$, then a straightforward solution to (7) is given by $\boldsymbol{\alpha}_{\tau}=\left[{\bf 0}_{\tau\times 1}^{T},1,{\bf 0}_{q-\tau\times 1}^{T}\right]^{T}$. From now on, we assume that $\tau>q$.

Writing (7) in the following equivalent way

[TABLE]

we can easily see that $\left\{\beta_{i}\right\}_{i=1}^{q}$ are roots of the following polynomial:

[TABLE]

Hence, there exists $Q(X)$ a polynomial with degree $\tau-q$ such that:

[TABLE]

Note that exact knowledge of $P(X)$ leads to the determination of the unknown coefficients $\alpha_{k,\tau}$, since they are by construction among the coefficients of $P(X)$. To fully characterize $P(X)$, we first observe that

• •

the coefficients of $P$ associated with exponents $X^{q},X^{q+1},\cdots,X^{\tau-1}$ are all zero.

• •

the coefficient associated with $X^{\tau}=-1$.

Let $\left\{a_{i}\right\}_{i=1}^{q+1}$ be the coefficients of $\prod_{i=1}^{q}\left(X-\beta_{i}\right)$ (i.e, $\prod_{i=1}^{q}\left(X-\beta_{i}\right)=\sum_{i=1}^{q+1}a_{i}X^{i-1}$), which can be exactly obtained using the Newton-Girard algorithm [6]. Let $\left\{b_{i}\right\}_{i=1}^{\tau-q+1}$ be the coefficients of $Q(X)$ so that: $Q\left(X\right)=\sum_{k=1}^{\tau-q+1}b_{k}X^{k-1}$. From the available information about the coefficients of $P$, we can show that $\left\{b_{k}\right\}_{k=1}^{\tau-q+1}$ satisfy the following set of equations

[TABLE]

where we use the convention that $a_{j}=0$ if $j\leq 0$ or $j>q+1$. The system of equation in (10) can be also expressed in the following matrix form:

[TABLE]

where $\mathbf{\Phi}$ is the upper triangular matrix given by

[TABLE]

Vector ${\bf b}=\left[b_{1},\cdots,b_{\tau-q+1}\right]^{T}$ can be thus determined by taking the inverse of matrix $\boldsymbol{\Phi}$ as

[TABLE]

From a numerical standpoint, this operation, involving inversion of an upper triangular matrix, can be solved in a stable fashion using back-substitution algorithm and is, as such, much more stable than the inversion of matrix $\boldsymbol{\Psi}$ required in the evaluation of (6). Once coefficients $\left\{b_{i}\right\}_{i=1}^{\tau-q+1}$ are obtained, $\left\{\alpha_{k,\tau}\right\}$ can be evaluated as 111This can be seen by using the fact that $P(X)=\sum_{i=1}^{q+1}\sum_{k=1}^{\tau-q+1}a_{i}b_{k}X^{i+k-2}$.

[TABLE]

To validate our procedure, we compute the positive moments of the Gram matrix $\mathbf{W}$ in the case where the correlation matrix $\mathbf{\Lambda}$ follows the following model [1]

[TABLE]

where the coefficient $\xi$ indicates the forgetting factor. This kind of matrices arise in covariance matrix estimation and more precisely in exponentially weighted sample covariance matrix (more details can be found in [1], section III-B) . Note that for moderate to large values of $q$, the eigenvalues of $\mathbf{\Lambda}$ given by $\left(1-\xi\right),\left(1-\xi\right)\xi,\cdots,\left(1-\xi\right)\xi^{q-1}$ are very close to each other, which might cause singularity issues when using the formula in (II-A). We consider two different configurations $\text{config}_{1}$ and $\text{config}_{2}$ corresponding respectively to $(n_{t}=3,q=5)$ and $(n_{t}=3,q=20)$. For both configurations, we evaluate the moments using (II-A) and the proposed method. We compare the obtained moments with the empirical ones evaluated over $10^{6}$ realizations.

The results are summarized in Table I. As a first observation, we notice that our method provides very close results to the empirical moments while the evaluation of the moments using (II-A) becomes totally inaccurate in configuration $\text{config}_{2}$ associated with a higher $q$. This clearly demonstrates the efficiency and the accuracy of our method in calculating the positive moments.

III Moment-Based Approach for Density Approximation

In this section, we show that the knowledge of all positive moments $\mu_{\mathbf{W}}\left(k\right)$, $k=1,2,\cdots$ can be leveraged to approximate the PDF $f_{\lambda}$. In general, retrieving a positive PDF from the knowledge of all its moments is known as the Stieltjes moment [7, 8] problem. We say that a PDF is called M-determinate if it can be uniquely determined by its moments. A sufficient condition for a PDF to be M-determinate is given by the Krein and the Lin conditions summarized below

Theorem 1**.**

[8]** Let $f$ be a distribution defined in the real half-line $\left(0,\infty\right)$. If the following conditions are satisfied:

The Krein condition:

[TABLE] 2. 2.

The Lin condition: $f$ is differentiable and

[TABLE]

Then, $f_{\lambda}$ is M-determinate.

Proposition 1**.**

The PDF in (3) is M-determinate.

Proof.

See Appendix for a proof. ∎

Now that we prove that $f_{\lambda}$ is M-determinate, an approximation of the marginal density, involving laguerre polynomials, can be derived as [9]:

[TABLE]

where $c=\frac{\mu_{\mathbf{W}}\left(2\right)-\mu_{\mathbf{W}}\left(1\right)^{2}}{\mu_{\mathbf{W}}\left(1\right)}$, $\nu=\frac{\mu_{\mathbf{W}}\left(1\right)}{c}-1$,

[TABLE]

is the Laguerre polynomial of order $i$ in $\lambda$ and parameter $\nu$ and

[TABLE]

Truncating the series in (17) at order $K$ yields the following approximation for the marginal PDF

[TABLE]

The CDF can thus be approximated as follows

[TABLE]

IV Selected numerical Results

In this section, we investigate the accuracy of the proposed PDF and CDF moment-based approach approximation. To this end, we compare them with their empirical counterparts and those evaluated using the results in [3].

In Figure 1, we assume that $\mathbf{\Lambda}$ follows the same model as in (14) with $\xi=0.85$, $n_{t}=3$ and $q=5$. We compare the accuracy of our approach with the corresponding empirical density and the formula provided in [3]. It can be noticed that our approximation becomes more accurate by increasing the truncation order $K$. As evidenced from Figures 1 and 2, a good approximation can be achieved starting from $K=30$ for both PDF and CDF.

In Figures 3 and 4, we increase the value of $q$ to $q=20$. In this case, the formula provided in [3] presents severe numerical instability and thus could not be plotted in this case. On the other hand, our moment-based approach achieves a very good approximation starting from $K=30$. For $K=45$, we can see that we have a perfect match with the empirical PDF and CDF.

V Conclusion

In this paper, we propose a numerically stable method that efficiently compute the positive moments of one-side correlated Gram matrices. From a practical standpoint, these moments can be used to approximate the marginal distribution and CDF of the eigenvalues of Large Gram random matrices and thus constitute an efficient alternative to conventional methods which become highly inaccurate in high dimensional settings.

Appendix (Proof of Proposition 1)

Krein condition

We start by rewrite the PDF in (3) as

[TABLE]

where $\gamma_{k,l}=\frac{\beta_{l}^{q-n_{t}-1}\mathbf{\Psi}_{k,l}^{-1}}{n_{t}\>\Gamma(n_{t}-q+k)}$. Then,

[TABLE]

Thus,

[TABLE]

and

[TABLE]

Integrating the first term of the right-hand side term provides infinity while integrating the second term results in a finite value. Thus, the integral diverges to infinity which fulfills the Krein condition.

Lin condition

Using the modified expression in (22), we have

[TABLE]

Then,

[TABLE]

Thus,

[TABLE]

This completes the proof of the proposition.