Sparse-Based Estimation Performance for Partially Known Overcomplete Large-Systems

TL;DR

This paper develops a low-cost, accurate Bayesian lower bound for estimating sparse signals in large systems with known interfering subspaces, combining compressed sensing and random matrix theory.

Contribution

It introduces a novel closed-form Bayesian lower bound for sparse vector estimation that leverages CS and RMT, applicable in doubly asymptotic regimes.

Findings

Closed-form lower bounds accurately predict MSE of sparse estimators.

The bounds are computationally efficient and interpretable.

Joint estimation/rejection improves signal recovery in noisy scenarios.

Abstract

We assume the direct sum <A> o <B> for the signal subspace. As a result of post- measurement, a number of operational contexts presuppose the a priori knowledge of the LB -dimensional "interfering" subspace <B> and the goal is to estimate the LA am- plitudes corresponding to subspace <A>. Taking into account the knowledge of the orthogonal "interfering" subspace <B>\perp, the Bayesian estimation lower bound is de- rivedfortheLA-sparsevectorinthedoublyasymptoticscenario,i.e. N,LA,LB -> \infty with a finite asymptotic ratio. By jointly exploiting the Compressed Sensing (CS) and the Random Matrix Theory (RMT) frameworks, closed-form expressions for the lower bound on the estimation of the non-zero entries of a sparse vector of interest are derived and studied. The derived closed-form expressions enjoy several interesting features: (i) a simple interpretable expression, (ii) a very low computational cost especially in the doubly asymptotic scenario, (iii) an accurate prediction of the mean-square-error (MSE) of popular sparse-based estimators and (iv) the lower bound remains true for any amplitudes vector priors. Finally, several idealized scenarios are compared to the derived bound for a common output signal-to-noise-ratio (SNR) which shows the in- terest of the joint estimation/rejection methodology derived herein.

Equations

$$\mathbf{y}=\mathbf{\Psi }\mathbf{s}+\mathbf{n}=\mathbf{\Psi }\mathbf{\Phi }\mathbf{x}+\mathbf{n}$$

$$\bar{\mathbf{y}}={\mathbf{U}}_{{\mathbf{B}}^{\psi }}^T\mathbf{y}={\mathbf{U}}_{{\mathbf{B}}^{\psi }}^T\mathbf{H}\mathbf{x}+\bar{\mathbf{n}}$$

$$C_{{\mathbf{U}}_{{\mathbf{B}}^{\psi }}^T{\mathbf{A}}^{\psi }}$$

$$\frac{1}{L_A}\mathrm{Tr}\left\{{\left({\mathbf{F}}^T\mathbf{F}\right)}^{-1}\right\}$$

$$\frac{1}{N-L_B}\mathrm{Tr}\left\{{\mathbf{F}}^T\mathbf{F}\right\}$$

$$C_{[{\mathbf{A}}^{\psi }{\mathbf{B}}^{\psi }]}$$

$$=\frac{{\sigma }_{\alpha }^2}{\mathrm{SN}{\mathrm{R}}_0^{(\mathrm{na})}}\frac{\mathrm{Tr}\left\{{\left([{\mathbf{A}}^{\psi }{\mathbf{B}}^{\psi }{]}^T[{\mathbf{A}}^{\psi }{\mathbf{B}}^{\psi }]\right)}^{-1}\right\}}{L}$$

$$\cdot \frac{{\sigma }_{\alpha }^2\mathrm{Tr}\left\{\left({\mathbf{A}}^{\psi T}{\mathbf{A}}^{\psi }\right)\right\}+{\sigma }_{\beta }^2\mathrm{Tr}\left\{\left({\mathbf{B}}^{\psi T}{\mathbf{B}}^{\psi }\right)\right\}}{N},$$

$$C_{{\mathbf{A}}^{\psi }}=\frac{{\sigma }_1^2}{L_A}\mathrm{Tr}\left\{{\left({\mathbf{A}}^{\psi T}{\mathbf{A}}^{\psi }\right)}^{-1}\right\}=\frac{{\sigma }_{\alpha }^2}{\mathrm{SN}{\mathrm{R}}_1^{(\mathrm{na})}}\frac{\mathrm{Tr}\left\{\left({\mathbf{A}}^{\psi T}{\mathbf{A}}^{\psi }\right)\right\}}{N}\frac{\mathrm{Tr}\left\{{\left({\mathbf{A}}^{\psi T}{\mathbf{A}}^{\psi }\right)}^{-1}\right\}}{L_A},$$

$$S_{{\mu }_{mp}}(z)=\frac{-1}{z}+\frac{\tilde{\rho }}{z}\frac{S_{{\mu }_{mp}}(z)}{\left(1+S_{{\mu }_{mp}}(z)\right)}$$

$$\int f(x){\hat{\mu }}_{{\mathbf{F}}^T\mathbf{F}}(dx)$$

$$\frac{1}{L_A}\mathrm{Tr}\left\{{\left(\rho {\mathbf{F}}^T\mathbf{F}\right)}^k\right\}$$

Full text

††journal: Signal Processing

1 Introduction

The Compressive Sampling or Compressed Sensing (CS) is an attractive domain which gives new trends for people interested in sampling theory of sparse signals [1, 2, 3]. The CS theory states that a sparse signal, i.e., a signal that can be decomposed as few non-zero values in a given basis (Fourier, wavelets, etc.) can be sampled at a rate $T_{S}$ lower than the one predicted by the Shannon’s theory. This paradigm has been successfully exploited for solving ill-posed problems arising for instance in bio-medical analysis, RADAR detection, array processing, wireless communications and radioastronomy imaging. In the CS framework, it is well known that any matrix ${\bf H}$ of size $N\times L$ generated from an i.i.d. centered sub-Gaussian distribution with a variance of $1/N$ verifies the Restricted Isometry Property (RIP) [2] with a high probability [1]. On the other hand, the doubly asymptotic spectrum and the empirical moments of the product ${\bf H}^{T}{\bf H}$ have been extensively studied in the context of the Random Matrix Theory (RMT) [4].

In the literature, CS and RMT techniques are usually applied to the noisy linear model where there is no interfering signals. However, in a wide range of real life applications, the signal of interest is often corrupted by a partially known interfering signal and an additive noise (see [5, 6, 7, 8, 9] for instance). This context motivates this work. More specifically, the CS and the RMT frameworks will be associated to derive new analytical closed-form expressions for the Bayesian lower bound [10] on the estimation of a sparse amplitude vector [11] for the noisy linear model corrupted by a partially known interfering signal.

2 Compressed Sensing (CS) integrating an a priori knowledge

2.1 Definition of the CS model

Let ${\bf y}$ an observed vector of $N$ measurements corrupted by an additive white centered zero-mean, Gaussian circular noise vector of variance $\sigma^{2}$. The standard CS model [2, 1, 3] is defined according to

[TABLE]

where ${\bf\Psi}$ is the known measurement matrix of size $N\times K$ with $N<K$, the vector ${\bf s}={\bf\Phi}{\bf x}$ of size $K\times 1$ admits an $L$-sparse representation, denoted by ${\bf x}$, in the basis ${\bf\Phi}$ (which could be Fourier basis, Wavelets basis, canonical basis, etc.) with $L<N$ and where ${\bf H}\stackrel{{\scriptstyle{\rm def.}}}{{=}}{\bf\Psi}{\bf\Phi}$ is often called the overcomplete dictionary.

One of the main problems risen up by the theory of the Compressed Sensing relates to the minimum number of measurements $N$ needed for retrieving the $L$-sparse vector $\bf{x}$. To address this problem, the authors of [2, 1, 3] have defined the Restricted Isometry Property (RIP). A standard strategy, called universal design strategy to ensure that dictionary $\bf{H}$ satisfies the RIP condition with high probability, is to generate the i.i.d. entries of dictionary matrix $\bf{H}$ following a sub-Gaussian distribution with zero mean and variance $1/N$ [1].

2.2 Exploiting the "interfering" subspace knowledge

In many real life applications, we do have the knowledge of information given by the physics of the context. Those useful information help in tailoring models that precisely take into account the knowledge of particular frequencies [12] for spectral analysis purpose, spatial angles for array processing [5] or RADAR processing, and have demonstrated their power through biomedical analysis or radioastronomy imaging. So, we adopt the following "signal+interference" model ${\bf s}={\bf A}\boldsymbol{\alpha}+{\bf i}\ \mbox{with}\ {\bf i}={\bf B}\boldsymbol{\beta}$ where $[{\bf A}]_{k,\ell}=g(kT_{S}-\tau_{\ell})$ with $1\leq\ell\leq L_{A}$, $[{\bf B}]_{k,\ell^{\prime}}=g(kT_{S}-\tilde{\tau}_{\ell^{\prime}})$ with $1\leq^{\prime}\ell\leq L_{B}$ are the "steering matrices" parametrized by the regular discretization at rate $T_{S}$ of a known waveform $g(t)$ along the time space. More precisely, $\mathcal{T}=\{\tau_{\ell},1\leq\ell\leq L_{A}\}$ stands for the time-delays of the $L_{A}$ sources of interest $\boldsymbol{\alpha}$ and $\mathcal{\tilde{T}}=\{\tilde{\tau}_{\ell},1\leq\ell\leq L_{B}\}$ is associated to the $L_{B}$ interfering sources $\boldsymbol{\beta}$. In the sequel, it is assumed that $(i)$ $\langle{\bf A}\rangle$ and $\langle{\bf B}\rangle$ are two disjoint subspaces, meaning that there is no time overlapping between the sources of interest and of the interfering sources and $(ii)$ $\langle{\bf B}\rangle$ is known or previously estimated (matrix $\bf{A}$ and $\langle{\bf A}\rangle$ are unknown). For instance, the learning of $\langle{\bf B}\rangle$ is based on pre-estimation of the clutter echo time-delays in RADAR processing or by known strongly shining "calibrator stars" in radioastronomy imaging. The problem of interest is to estimate vector $\boldsymbol{\alpha}$ based on a measurement vector where the contribution of ${\bf i}$ has been removed using the knowledge of $\langle{\bf B}\rangle$. The standard "signal+interference" model described by signal ${\bf s}$ can be extended in the CS framework of model (1) following a straightforward strategy. Let ${\bf\Phi}$ be a basis matrix such as $[{\boldsymbol{\Phi}}]_{k,k^{\prime}}=g((k-k^{\prime})T_{S})$ where $1\leq k,k^{\prime}\leq K$. For a sufficiently fine partition, i.e., for $K>N>L=L_{A}+L_{B}$, we have $\langle{\bf A}\rangle\oplus\langle{\bf B}\rangle\subset\langle{\bf\Phi}\rangle$. Let ${\bf U}_{{\bf B}^{\psi}}$ be a $N\times(N-L_{B})$ orthonormal basis matrix such as $\langle{\bf U}_{{\bf B}^{\psi}}\rangle=\langle{\bf B}^{\psi}\rangle^{\perp}$. We have finally the deflated observation, defined according to

[TABLE]

where ${\bf\bar{n}}={\bf U}_{{\bf B}^{\psi}}^{T}{\bf n}$ and ${\bf x}$ is a $(K-L_{A})$-sparse such as ${\bf x}_{\mathcal{T}}=\boldsymbol{\alpha}$. The reader will find an illustration of the procedure in Fig. 1.

3 ECRB for projected measurements and a large random dictionary

3.1 Dealing with projected measurements

Let ${\rm MSE}=\frac{1}{L_{A}}\mathbb{E}_{\bar{{\bf y}},\boldsymbol{\alpha}}\left[\|\boldsymbol{\hat{\alpha}}({\bf\bar{y}})-\boldsymbol{\alpha}\|^{2}\right]$ be the normalized Bayesian Mean Squared Error for an estimate $\boldsymbol{\hat{\alpha}}({\bf\bar{y}})$ of $\boldsymbol{\alpha}$. The Expected Cramér-Rao Bound (ECRB) [10], denoted by $C_{\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}}$ for the random amplitude vector $\boldsymbol{\alpha}$, of unspecified distribution $p(\boldsymbol{\alpha})$ given the observation model (2) fulfills relation ${\rm MSE}\geq C_{\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}}=\frac{\sigma^{2}}{L_{A}}\ {\rm Tr}\left\{\left(\mathbf{A}^{\psi T}\mathbf{P}_{\mathbf{B^{\psi}}}^{\perp}\bf{A}^{\psi}\right)^{-1}\right\}$ where ${\bf A}^{\psi}=\boldsymbol{\Psi}{\bf A}$. Introduce model $(\mathcal{M})$: $\bar{\bf{y}}|\boldsymbol{\alpha}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}),$ where $\boldsymbol{\mu}={\bf U}_{{\bf B}^{\psi}}^{T}{\bf A}^{\psi}\boldsymbol{\alpha}$ and $\boldsymbol{\Sigma}=\sigma^{2}{\bf I}_{N-L_{B}}$ which is the covariance matrix of noise ${\bf\bar{n}}$. After some calculus, the ECRB admits the following expression:

[TABLE]

where ${\rm SNR^{(na)}}=\frac{\mathbb{E}||\boldsymbol{\mu}||^{2}}{{\rm Tr}\{\boldsymbol{\Sigma}\}}=\frac{\sigma_{\alpha}^{2}{\rm Tr}\{\mathbf{A}^{\psi T}\mathbf{P}_{\mathbf{B^{\psi}}}^{\perp}\bf{A}^{\psi}\}}{\sigma^{2}(N-L_{B})}$ is the output and non-asymptotic SNR.

3.2 Doubly asymptotic regime

The practical interest of CRB-type expressions have been exposed in [13, 14] but we show in this work that expression (3) can be reduced to a very simple closed form expression with the advantage of remaining valid even for the low sample regime, using some powerful results extracted from the RMT [4] where it is assumed $N,L_{A},L_{B}\rightarrow\infty$ with $N/L_{A}\rightarrow\rho$ and $L_{B}/L_{A}\rightarrow c$. Towards this goal, the following Lemma is provided.

Lemma 1**.**

Let $\mathbf{F}=\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}\in\mathbb{R}^{(N-L_{B})\times L_{A}}$ whose elements $\left\{F_{ij}\right\}_{i,j=1\ldots N-L_{B},L_{A}}$ are zero mean and i.i.d. with variance $\frac{1}{N}$. Now, for $N,L_{A},L_{B}\rightarrow\infty$, and $N/L_{A}\rightarrow\rho>1$, $(N-L_{B})/L_{A}\rightarrow\tilde{\rho}=\rho-c>1$, then

[TABLE]

where $a.s.$ stands for the almost sure convergence.

Proof.

See the appendix. ∎

Under the assumptions of Lemma 1 and using (3), a very compact expression of $C_{\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}}^{\infty}$ is enunciated by the following.

Result 1**.**

Assume that $N,L_{A},L_{B}\rightarrow\infty$ and $N/L_{A}\rightarrow\rho>1$, $(N-L_{B})/L_{A}\rightarrow\tilde{\rho}>1$, then, we have $C_{\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}}\overset{a.s.}{\longrightarrow}C^{\infty}_{\mathbf{U}_{\mathbf{B}^{\psi}}^{T}\mathbf{A}^{\psi}}=\frac{\sigma_{\alpha}^{2}}{{\rm SNR}}\frac{1}{\tilde{\rho}-1}$ where ${\rm SNR}=\frac{\sigma_{\alpha}^{2}}{\sigma^{2}\rho}$ is the almost sure doubly asymptotic equivalent of ${\rm SNR^{(na)}}$.

4 Benchmarking ECRBs and estimators

This section is devoted to give a relation of order between the ECRB given by (3) with respect to two other ECRBs viewed as benchmarks and to analyze the behavior of sparse-based estimators. Let $(\mathcal{M}_{0}):{\bf y}_{0}|\boldsymbol{\alpha},\boldsymbol{\beta}\sim\mathcal{N}\left(\mathbf{A}^{\psi}\boldsymbol{\alpha}+\mathbf{B}^{\psi}\boldsymbol{\beta},\sigma_{0}^{2}{\bf I}_{N}\right)$ and $(\mathcal{M}_{1}):{\bf y}_{1}|\boldsymbol{\alpha}\sim\mathcal{N}\left(\mathbf{A}^{\psi}\boldsymbol{\alpha},\sigma_{1}^{2}{\bf I}_{N}\right)$. Model $\mathcal{M}_{0}$ is associated with the scenario where no ad-hoc strategy is developed to mitigate the corruption from the interference signals. In other words, the interference signals are wrongly interpreted as signals of interest. So, this bound does not solve the problem of interest and is given by

[TABLE]

where ${\rm SNR_{0}^{(na)}}=\frac{\sigma_{\alpha}^{2}\ {\rm Tr}\left\{\left(\mathbf{A}^{\psi T}\bf{A}^{\psi}\right)\right\}\ +\ \sigma_{\beta}^{2}\ {\rm Tr}\left\{\left(\mathbf{B}^{\psi T}\bf{B}^{\psi}\right)\right\}}{\sigma^{2}_{0}N}$ and ${\rm SIR}=\sigma_{\alpha}^{2}/\sigma_{\beta}^{2}$. The second model $\mathcal{M}_{1}$ is associated with the ideal free-interference scenario. This bound admits the following expression:

[TABLE]

where ${\rm SNR_{1}^{(na)}}=\frac{\sigma_{\alpha}^{2}\ {\rm Tr}\left\{\left(\mathbf{A}^{\psi T}\mathbf{A}^{\psi}\right)\right\}}{\sigma^{2}_{1}N}.$ Using a similar methodology as before, we have the following result given in the doubly asymptotic regime context.

Result 2**.**

Assume that $N,L_{A}\to\infty$ with $N/L_{A}\rightarrow\rho>1$, $N/L\rightarrow\bar{\rho}>1$, then, we have $C_{[\mathbf{A}^{\psi}\mathbf{B}^{\psi}]}\overset{a.s.}{\longrightarrow}C_{[\mathbf{A}^{\psi}\mathbf{B}^{\psi}]}^{\infty}=\frac{\sigma_{\alpha}^{2}}{{\rm SNR}_{0}}\left(\frac{1-{\rm SIR}^{-1}}{\rho}+\frac{{\rm SIR}^{-1}}{\bar{\rho}}\right)\frac{\bar{\rho}}{\bar{\rho}-1}$ and $C_{\mathbf{A}^{\psi}}\overset{a.s.}{\longrightarrow}C_{\mathbf{A}^{\psi}}^{\infty}=\frac{\sigma_{\alpha}^{2}}{{\rm SNR}_{1}}\frac{1}{\rho-1}$ where ${\rm SNR}_{0}=\frac{(\sigma_{\alpha}^{2}-\sigma_{\beta}^{2})}{\sigma_{0}^{2}\rho}+\frac{\sigma_{\beta}^{2}}{\sigma_{0}^{2}\bar{\rho}}$ and ${\rm SNR}_{1}=\frac{\sigma_{\alpha}^{2}}{\sigma_{1}^{2}\rho}$ are the almost sure doubly asymptotic equivalent of $\;{\rm SNR_{0}^{(na)}}$ and $\;{\rm SNR_{1}^{(na)}}$, respectively.

4.1 Simulation Tests

We propose now different scenarios for (i) numerically showing the limit of the doubly asymptotic approximation and (ii) for analyzing the behavior of popular sparse-based estimators when they are tailored with the account of the interference signal. We have plotted on Fig. 2 the Mean Square Error of the approximation between analytic expressions (3), (7), (8) of the ECRBs and their doubly asymptotic expressions (Results 1 and 2) with respect to the number of samples. This figure leaves no room for doubt about the speed at which the asymptotic expressions for the ECRBs merge the theoretical ones. Finally, we have computed (Fig. 3) for severals growing dimensions of the $L_{B}$-dimensional subspace $\langle\bf{B}\rangle$ and six punctual values of $L_{A}$, selected from the range $[0.01N:0.3N]$, the ECRBs. Like Fig. 2, we clearly see that both plots merge whatever the values for $L_{A}$ and $L_{B}$.

We now turn the discussion to the behavior of the tailored estimators. To be sufficiently general, we have not designed the observation model with specific waveforms or delays but directly implemented dictionary $\mathbf{H}$ as presented by section 2.1. The Basis Pursuit DeNoise [15], the CoSaMP [16] and the Orthogonal Matching Pursuit [17] sparse estimators have been computed with a deflated observation signal and a deflated dictionary for them to estimate the $L_{A}$ parameters of interest only. We have run 500 Monte-Carlo trials for each scenario with a common output SNR and confronted the Mean-Square-Error (MSE) with the ECRBs (3), (7), (8) and their doubly asymptotic equivalents given in Results 1 and 2 for both scenarios drawn through Fig. 4. When the observation signal is composed by as much equally powered sources of interest as interfering ones (Fig. 4), we observe that $C_{\mathbf{U}_{\mathbf{B}^{T}}\mathbf{A}}$ is clearly next to $C_{\mathbf{A}}$. Then, a quick analysis on both plots, reveals that none of the algorithms can perform well in a low SNR regime. The interfering signal has been almost rejected in that situation and only the OMP estimator seems to be in capacity to do from a 20 dB SNR. The CoSaMP performs as well as the interfering signal was of interest, it is so unable to reject the influence of the interfering part and the BPDN does not reach any lower bounds. Curiously thereafter, when the signal of interest is buried into many interfering and much powerful elements (Fig. 5), we notice a significant gain (more than 15 dB) for each estimators when they account for the knowledge of the interferences. In that difficult scenario, the OMP still reaches from a 20 dB SNR the lower bound $C_{\mathbf{U}_{\mathbf{B}^{T}}\mathbf{A}}$ which still remains close to the ideal bound $C_{\mathbf{A}}$ for which there is no interfering elements.

5 Conclusion

In this work, the problem of interest is to derive the Bayesian performance bound for the estimation of the $L_{A}$ non-zero amplitudes of a sparse signal of interest based on the observation of a $N\times 1$ compressed measurement vector ${\bf y}$. Vector ${\bf y}$ follows an overcomplete Bayesian linear model corrupted by a set of interfering signals spanning an a priori known $L_{B}$-dimensional subspace. Based on this standard assumption, the measurement vector ${\bf y}$ is confined in a $N-L_{B}$ subspace thanks to an orthogonal deflation technique. In addition, the proposed analysis is done in the asymptotic framework, i.e., for $N,L_{A},L_{B}\rightarrow\infty$ with finite asymptotic ratios. Our methodology allows to obtain $(i)$ an easily interpretable expression of the bound, $(ii)$ a cheap computational cost especially in the doubly asymptotic scenario, $(iii)$ an accurate prediction of the mean-square-error (MSE) of popular sparse-based estimators and $(iv)$ a lower bound for any amplitudes vector priors. Finally, several idealized scenarios are compared to the derived bound for a common output signal-to-noise-ratio (SNR) which shows the interest of our joint estimation/rejection methodology.

6 Appendix : Proof of Lemma 1

We first need to introduce the normalized trace of the resolvent for the random form $\rho\mathbf{F}^{T}\bf{F}$ by the complex function $h(z)=\frac{1}{L_{A}}{\rm Tr}\left\{\left(\rho\mathbf{F}^{T}\mathbf{F}-z\mathbf{I}\right)^{-1}\right\}.$ Owing to the assumptions stated by Lemma 1 and due to [18], we know that the empirical distribution of $\rho\mathbf{F}^{T}\bf{F}$ defined by $\hat{\mu}_{\rho\mathbf{F}^{T}\bf{F}}=\frac{1}{L_{A}}\sum_{i=1}^{L_{A}}\delta_{\lambda_{i}(\rho\mathbf{F}^{T}\bf{F})}$ with $\delta_{\lambda_{i}}$ the Dirac measure for eigenvalue $\lambda_{i}$, that is the unique probability measure satisfying $f(x)\delta_{\lambda}=f(\lambda)$ for any continuous function $f\in\mathbb{R}$; converges almost surely in distribution towards a deterministic distribution function $\mu_{mp}$, i.e. $\hat{\mu}_{\rho\mathbf{F}^{T}\bf{F}}\overset{a.s.}{\longrightarrow}\mu_{mp}$. The deterministic distribution $\mu_{mp}$ is supported on the compact interval $[\lambda_{-},\lambda_{+}]$ possibly with a point mass at 0 [4] with generalized density $\frac{d\mu_{mp}(x)}{dx}=\frac{\sqrt{\left(\lambda_{+}-x\right)\left(x-\lambda_{-}\right)}}{2\pi\tilde{\rho}x}$ with $\lambda_{\pm}=(1\pm\sqrt{\tilde{\rho}})^{2}$ and referred to as the Marchenko-Pastur density [18]. For any distribution function $\mu$, the Stieltjes Transform of $\mu$ denoted by $S_{\mu}(z):\mathbb{C}\;\textbackslash$ Supp$(\mu)\longrightarrow\mathbb{C}$ is defined as $S_{\mu}(z)=\int\frac{\mu(d\lambda)}{\lambda-z}.$ It has been already shown [19] that the Stieltjes Transform $S_{\mu_{mp}}(z)$ of $\frac{d\mu_{mp}(x)}{dx}$ respects the following quadratic relation

[TABLE]

and proved to converge [19], towards $h(z)$. Notice now that the Stieltjes Transform $S_{\mu_{mp}}(z)$ weighted by $\frac{L_{A}}{N}\rightarrow\rho^{-1}$ and evaluated at $z=0$ corresponds to $\frac{1}{L_{A}}{\rm Tr}\left\{\left(\mathbf{F}^{T}\bf{F}\right)^{-1}\right\}$. With all these materials, (4) of lemma 1 is obtained by expressing the weighted Stieltjes Transform $S_{\mu_{mp}}(0)$ through the above expression. To prove (5) we use the spectral theorem applied to empirical distribution $\hat{\mu}_{\rho\bf{F}^{T}\bf{F}}$, to obtain

[TABLE]

Since this empirical distribution converges when weighted by $\rho$ towards the Marchenko-Pastur distribution, i.e. $\hat{\mu}_{\rho\mathbf{F}^{T}\bf{F}}\overset{a.s.}{\longrightarrow}\mu_{mp}$, we have consequently for $f$ a polynomial function the assertion

[TABLE]

which reduces obviously to $\tilde{\rho}$ for $k=1$. We have consequently $\frac{1}{L_{A}}{\rm Tr}\left\{\left(\mathbf{F}^{T}\bf{F}\right)^{k}\right\}\overset{a.s.}{\longrightarrow}\rho^{-1}\tilde{\rho}$ and after obvious manipulations (5) is reached.