An Asymptotically MSE-Optimal Estimator
based on Gaussian Mixture Models

Michael Koller, Benedikt Fesl, Nurettin Turan, and
Wolfgang Utschick, Fellow, IEEE
This work was supported by the Deutsche Forschungsgemeinschaft under
grant UT 36/21.
Preliminary results have been published at ICASSP'22 <cit.>.
The authors are with Professur für Methoden der Signalverarbeitung, Technische Universität München,
80333 München, Germany, e-mail<EMAIL_ADDRESS>
Koller et al.

This paper investigates a channel estimator based on gmm in the context of linear inverse problems with additive Gaussian noise.
We fit a gmm to given channel samples to obtain an analytic pdf which approximates the true channel pdf.
Then, a cme corresponding to this approximating pdf is computed in closed form and used as an approximation of the optimal cme based on the true channel pdf.
This optimal cme cannot be calculated analytically because the true channel pdf is generally unknown.
We present mild conditions which allow us to prove the convergence of the gmm-based cme to the optimal cme as the number of gmm components is increased.
Additionally, we investigate the estimator's computational complexity and present simplifications based on common model-based insights.
Further, we study the estimator's behavior in numerical experiments including mimo and wideband systems.

asymptotic convergence,
conditional mean channel estimation,
Gaussian mixture models,
machine learning,
spatial channel model

§ INTRODUCTION

Channel estimation plays a critical role in future mobile communications systems, e.g., <cit.>.
The mse minimizing channel estimator is known as cme.
Computing the cme in closed form requires analytic knowledge of the channel pdf.
Even if the pdf was known, calculating the cme might not be possible analytically or not be tractable practically.
Increasingly, advanced channel models (e.g., <cit.>) or simulators (e.g., <cit.>) are used to generate large amounts of realistic channel samples.
In a real application, channel samples can, for example, be collected at the base station to be used in addition to or instead of the simulated data.
Importantly, such data represent the whole scenario (or environment) in which the base station is placed.
It is thus interesting to investigate data-based algorithms to design channel estimators which are applicable to a whole scenario.
Many estimators have been proposed in this context
(cf., e.g., <Ref> for details).
In particular, estimators based on cs and on machine learning have recently been proposed (see, e.g., <cit.>).
To our knowledge, the algorithms' optimality has not been studied if an arbitrary channel pdf is assumed.

In this paper, we study a gmm-based channel estimator.
The estimator itself has already been investigated in the case where the channel pdf is given by a gmm.
One of our contributions is to provide a strong motivation to employ the gmm-based estimator even if the channel pdf is not a gmm.
To this end, we show (in a proof and in numerical simulations) that even if the channel is not gmm distributed,
the gmm-based estimator converges to the optimal cme as the number of gmm components increases.

In detail, the following approach is taken in this paper.
First, channel samples are used to fit a gmm.
Since gmm can approximate any continuous pdf <cit.>,
the fitted gmm is a pdf which approximates the unknown true channel pdf.
Second, we analytically compute a cme for channels distributed according to the gmm pdf.
Since the gmm pdf approximates the true channel pdf,
we ask whether the gmm-based cme approximates the true cme.

A related work is <cit.>, where the authors assume that the channel is gmm-distributed and study the available closed-form cme for example in the asymptotic high snr regime to derive pilot signals.
In this work, we do not assume that the channel is gmm-distributed.
Instead, a gmm is used as an approximation of the true channel pdf and we analyze whether the corresponding cme is an approximation of the true cme.
This can be viewed as a study of the gmm estimator in the high number of gmm components regime.
A main contribution of our paper is to prove that as the number of components increases, the gmm cme converges to the optimal cme if the observation matrix is invertible (cf. <Ref>).
For noninvertible observation matrices, we make a weaker statement.
Moreover, we analyze the gmm estimator's computational complexity and show how the complexity can be reduced in different estimation scenarios.

We study the gmm estimator in numerical simulations where we consider both mimo and wideband channel estimation scenarios with both invertible and noninvertible observation matrices.
The considered channel data come from a 3GPP channel model <cit.> and from the QuaDRiGa channel simulator <cit.> so that they are not gmm-distributed by construction.
The generated data represent a scenario where for example a base station covers a certain sector with users whose positions are drawn uniformly at random.
The obtained gmm estimator is then suited for channel estimation in the whole scenario.
Already for a finite number of gmm components,
the gmm estimator shows a performance close to the optimal cme in the numerical simulations.
We emphasize the gmm estimator's broad applicability by comparing it to state-of-the-art algorithms from the literature.

The paper is structured as follows.
<Ref> introduces the signal model discussed throughout the paper as well as particular instances thereof which are used in numerical simulations.
<Ref> reviews gmm and channel estimation literature which employs them.
The main part is <Ref> where we investigate the gmm-based cme and study its convergence to the optimal cme as well as its computational complexity.
<Ref> present state-of-the-art channel estimation algorithms, channel models, and numerical simulations.

The supremum norm of a continuous function \( f: \R^N \to \R \) is given by \( \| f \|_\infty = \sup_{\mbx\in\R^N} |f(\mbx)| \),
and \( \| \mbx \| \) is the Euclidean norm of \( \mbx \in \C^N \).
A real- or complex-valued normal distribution with mean vector \( \mbmu \) and covariance matrix \( \mbC \) is denoted by \( \calN(\mbmu, \mbC) \) or \( \calN_{\C}(\mbmu, \mbC) \), respectively.
The vectorization (stacking columns) of a matrix \( \mbX \in \C^{m\times N} \) is written as \( \vect(\mbX) \in \C^{mN} \), and \( \mbA \otimes \mbB \in \C^{m_1m_2\times N_1N_2} \) is the Kronecker product of \( \mbA \in \C^{m_1\times N_1} \) and \( \mbB \in \C^{m_2\times N_2} \).

§ SIGNAL MODELS

We consider the generic signal model

\begin{equation}\label{eq:signal_model}
\mby
= \mbA \mbh + \mbn, \quad\quad \mbn \sim \calN_{\C}(\mbzero, \mbSigma)
\end{equation}

where \( \mbh \in \C^N \) is the channel, \( \mbA \in \C^{m\times N} \) is the observation matrix, and \( \mbn \in \C^m \) is additive white Gaussian noise.
The technical interpretation (e.g., number of antennas or pilots) of the dimensions \( m \) and \( N \) depends on the context.
Examples can be found in the following subsections.
The observation \( \mby \in \C^m \), the matrix \( \mbA \), the noise mean vector \( \mbzero \in \C^m \), and the noise covariance matrix \( \mbSigma \in \C^{m\times m} \) are given.
The goal of channel estimation is to recover \( \mbh \) from (<ref>).

In this paper, we study a channel estimation algorithm which is designed using the given signal model as well as a data set of channel samples.
While the main part addresses the generic signal model (<ref>), we consider the following three instances of it in numerical simulations.

§.§ Single-Input Multiple-Output Signal Model

The simo signal model is for instance appropriate if a single-antenna mobile device transmits pilot signals to a base station with \( N \) antennas which receives

\begin{equation}\label{eq:signal_model_simo}
\mby = \mbh + \mbn \in \C^N.
\end{equation}

This model is interesting for us because the observation matrix is the identity matrix and therefore invertible.
Further, the performance of the proposed channel estimator can be studied without having to take into account the difficulty of choosing a suitable observation matrix.

§.§ Multiple-Input Multiple-Output Signal Model

If a mobile user with \( \Ntx \) antennas transmits \( N_p \) pilots to a base station with \( \Nrx \) antennas,
the receive signal \( \mbY \in \C^{\Nrx\times N_p} \) can be written as

\begin{equation}\label{eq:signal_model_mimo}
\mbY = \mbH \mbP + \mbN
\end{equation}

where \( \mbH \in \C^{\Nrx\times \Ntx} \) is the channel, \( \mbP \in \C^{\Ntx\times N_p} \) is the pilot, and \( \mbN \in \C^{\Nrx \times N_p} \) is the noise matrix.
With the definitions \( \mbh = \vect(\mbH) \), \( \mby = \vect(\mbY) \), \( \mbn = \vect(\mbN) \), and \( \mbA = \mbP^\tp \otimes \mbI_{\Nrx} \), the mimo signal model (<ref>) is an instance of (<ref>).

§.§ Wideband Signal Model

If we consider a siso transmission in the spatial domain over a frequency-selective fading channel, \( \mbH\in \C^{N_c\times N_t} \) represents the time-frequency response of the channel for \( N_c \) subcarriers and \( N_t \) time slots.
When only \( N_p \) positions of the time-frequency response are occupied with pilot symbols,
then there is a selection matrix \( \mbA \in \{0,1\}^{N_p\times N_c N_t} \) which represents the pilot positions.
This leads to the observations as described in (<ref>) with \( \mbh =\vect(\mbH)\in\C^{N_c N_t} \).
Regarding the structure of the pilot positions, three different arrangements are commonly considered: block-, comb-, and lattice-type, cf. <cit.>.

§ GAUSSIAN MIXTURE MODELS IN THE LITERATURE

In this section, we briefly explain gmm and summarize channel estimation literature which makes use of gmm.

§.§ Gaussian Mixture Models

A gmm with \( K \) components is a pdf of the form <cit.>

\begin{equation}\label{eq:gmm_of_h}
f_{\mbh}^{(K)}(\mbh) = \sum_{k=1}^K p(k) \calN_{\C}(\mbh; \mbmu_k, \mbC_k)
\end{equation}

consisting of a weighted sum of \( K \) Gaussian pdf.
The probabilities \( p(k) \) are called mixing coefficients, and \( \mbmu_k \in \C^N \) and \( \mbC_k \in \C^{N\times N} \) denote the mean vector and covariance matrix of the \( k \)th gmm component, respectively.
As explained in <cit.>, gmm allow to calculate the responsibilities \( p(k \mid \mbh) \) by evaluating Gaussian likelihoods:

\begin{equation}
p(k \mid \mbh) = \frac{p(k) \calN_{\C}(\mbh; \mbmu_k, \mbC_k)}{\sum_{i=1}^K p(i) \calN_{\C}(\mbh; \mbmu_i, \mbC_i) }.
\end{equation}

This is an important property for our considerations.

Given data samples, an em algorithm can be used to fit a \( K \)-components gmm <cit.>.
The data-based fitting process determines the mixing coefficients, the mean vectors, and the covariance matrices.
A detailed introduction to gmm and the corresponding well-known em algorithm can, e.g., be found in <cit.>.

§.§ Gaussian Mixture Models in Channel Estimation Literature

In <cit.>, channels are modeled as sparse vectors whose non-zero coefficients are gmm-distributed.
A combination of em and amp is then introduced for channel estimation.
The algorithm simultaneously estimates the gmm parameters.
Building on this work, <cit.> models the beam domain channels via gmm in the context of uplink channel estimation with pilot contamination.
The authors of <cit.> then investigate the approach further including a new initialization technique for the algorithm.
In <cit.>, the beam domain channel is also assumed to be gmm-distributed and a modification of learned amp is proposed for sparse channel estimation.

The authors of <cit.> model temporal channel variations as gmm, e.g., in order to predict channel states.
The authors of <cit.> propose to improve channel estimation techniques by using gmm as a better characterization of the noise in communications environments than it is given by the additive white Gaussian noise model.
In <cit.>, a gmm prior is used for the unknown data symbols in semi-blind channel estimation.

gmm are also employed for channel clustering tasks.
For example, <cit.> use gmm for channel multipath clustering.
In <cit.>, a power weighted gmm is proposed to increase the clustering performance.
Another variation of gmm, called rotationally invariant gmm, can be found in <cit.>.

In <cit.>, the true channel pdf is assumed to be equal to a gmm and the authors then investigate the corresponding cme to optimize the pilot matrix.
To this end, the asymptotic high-snr regime of the cme is studied.
Further, an information-theoretic criterion for pilot optimization is introduced because the mse of the estimator has no closed-form expression and is thus not suitable as optimization criterion <cit.>.

In this paper, we study the same gmm-based estimator as the authors of <cit.>.
However, we do not assume that the true channel pdf is equal to a gmm.
Instead, we take the gmm as an approximation of the true channel pdf and we ask whether the corresponding gmm-based estimator is then an approximation of the true cme (based on the true channel pdf).
In this sense, we study the gmm-based estimator's behavior in the high number of components regime and our work complements <cit.> by motivating the application of the estimator in a wider class of channel models.

§ MAIN PART

The mse-optimal channel estimate for the model (<ref>) is given by the conditional expectation \( \expec[\mbh \mid \mby] \), cf., e.g., <cit.>.
However, the true channel pdf is generally not known and, therefore, \( \expec[\mbh \mid \mby] \) can generally not be calculated analytically.
Even if the true channel pdf was known, the cme \( \expec[\mbh \mid \mby] \) might still not have an analytic expression.
In this section, we investigate a gmm-based cme with closed-form expression and prove that it converges to the optimal cme as the number of gmm components is increased.
Further, we discuss its computational complexity and how the complexity can be reduced.

§.§ Channel Estimator

gmm are known to be able to approximate any continuous pdf arbitrarily well <cit.>.
In particular, if \( f_{\mbh} \) denotes the pdf of the channel which is assumed to be continuous, then there exists a sequence \( ( f_{\mbh}^{(K)} )_{K=1}^\infty \) of gmm which converges uniformly to \( f_{\mbh} \).
To define a gmm-based estimator, let \( f_{\mbn} \) and \( f_{\mby} \) be the pdf of the noise and the observation, respectively, and let us first observe the following:

\begin{equation}\label{eq:fh_given_y}
f_{\mbh\mid\mby}(\mbh\mid\mby) = \frac{f_{\mby\mid\mbh}(\mby\mid\mbh)f_{\mbh}(\mbh)}{f_{\mby}(\mby)} = \frac{f_{\mbn}(\mby - \mbA\mbh)f_{\mbh}(\mbh)}{f_{\mby}(\mby)}.
\end{equation}

With this the optimal cme can be expressed as

\begin{equation}\label{eq:conditional_mean}
\hhat(\mby) = \expec[\mbh \mid \mby]
= \int \mbh \frac{f_{\mbn}(\mby - \mbA\mbh)f_{\mbh}(\mbh)}{f_{\mby}(\mby)} d \mbh.
\end{equation}

For every \( K \in \N \), we now consider the model

\begin{equation}
\mby^{(K)} = \mbA \mbh^{(K)} + \mbn
\end{equation}

where \( \mbh^{(K)} \) is distributed according to the gmm \( f_{\mbh}^{(K)} \)
which has the form (<ref>).
Since we have a sequence \( (f_{\mbh}^{(K)} )_{K=1}^\infty \) of gmm, the parameters \(p(k)\), \(\mbmu_k\), and \(\mbC_k\) would also depend on the sequence index \(K\) but we omit it for readability.
Let \( f_{\mby}^{(K)} \) be the pdf of \( \mby^{(K)} \).
We now define a gmm-based estimator

\begin{equation}\label{eq:conditional_mean_K}
\hhat^{(K)}(\mby) := \expec^{(K)}[\mbh^{(K)} \mid \mby]
:=\int \mbh \frac{f_{\mbn}(\mby - \mbA\mbh)f_{\mbh}^{(K)}(\mbh)}{f_{\mby}^{(K)}(\mby)} d \mbh
\end{equation}

by replacing \( f_{\mbh} \) and \( f_{\mby} \) in (<ref>) with \( f_{\mbh}^{(K)} \) and \( f_{\mby}^{(K)} \), respectively,
because similar to (<ref>) we have

\begin{equation}\label{eq:fhk_given_y}
= \frac{f_{\mbn}(\mby - \mbA\mbh)f_{\mbh}^{(K)}(\mbh)}{f_{\mby}^{(K)}(\mby)}.
\end{equation}

The law of total expectation allows us to write

\begin{equation}\label{eq:total_expectation}
\hhat^{(K)}(\mby) = \sum_{k=1}^K p(k\mid \mby) \expec^{(K)}[\mbh^{(K)} \mid \mby, k]
\end{equation}

in order to introduce the gmm mixing variable.
By definition of gmm, conditioning on one of the mixing variables yields a Gaussian random vector.
That is, \( \mbh^{(K)} \mid k \sim \calN_{\C}(\mbmu_k, \mbC_k) \) is the \( k \)th Gaussian in the gmm \( f_{\mbh}^{(K)} \), see also (<ref>).
Since \( \mbh^{(K)} \mid k \) is Gaussian, also \( \mbA \mbh^{(K)} \mid k \) and therefore \( \mby^{(K)} \mid k \) are Gaussian.
The conditional mean vector and conditional covariance matrix of \( \mby^{(K)} \mid k \) are \( \expec[\mby^{(K)} \mid k] = \mbA \meanhk \) and

\begin{equation}
\expec[(\mby^{(K)} - \mbA \meanhk)(\mby^{(K)} - \mbA \meanhk)^\herm \mid k]
= \mbA \covhk \mbA^\herm + \mbSigma,
\end{equation}

The well-known lmmse formula can now be used to compute

\begin{equation}\label{eq:lmmse_formula}
\expec^{(K)}[\mbh^{(K)} \mid \mby, k] =
\covhk \mbA^\herm (\mbA \covhk \mbA^\herm + \mbSigma)^{-1} (\mby - \mbA \meanhk) + \meanhk
\end{equation}

which can be plugged into (<ref>).

In order to calculate \( p(k \mid \mby) \) in (<ref>), we compute the pdf

\begin{equation}\label{eq:gmm_y}
f_{\mby}^{(K)}(\mby) = \sum_{k=1}^K p(k) \calN_{\C}(\mby; \mbA \meanhk, \mbA \covhk \mbA^\herm + \mbSigma),
\end{equation}

which is a gmm.
gmm allow to calculate the responsibilities by evaluating Gaussian likelihoods (cf. <Ref>):

\begin{equation}\label{eq:responsibilities}
p(k \mid \mby) = \frac{p(k) \calN_{\C}(\mby; \mbA \meanhk, \mbA \covhk \mbA^\herm + \mbSigma)}{\sum_{i=1}^K p(i) \calN_{\C}(\mby; \mbA \meanhi, \mbA \covhi \mbA^\herm + \mbSigma) }.
\end{equation}

Plugging this into (<ref>) shows that as soon as the mixing coefficients \( p(k) \) as well as the means \( \meanhk \) and covariances \( \covhk \) are given, the estimator \( \hest \) can be computed in closed form by combining (<ref>), (<ref>), and (<ref>), which results in (<ref>).
As discussed next in <Ref>, various quantities of the gmm estimator can be precomputed at this point to save computational complexity:
the products involving the known observation matrix \( \mbA \) and means \( \mbmu_k \) and covariance matrices \(\mbC_k\) as well as the lmmse filters including the computationally costly matrix inverse in (<ref>).
To obtain the gmm parameters, the channel pdf \( f_{\mbh} \) needs to be approximated by fitting a \( K \)-components gmm to given channel samples, cf. <Ref>.

\begin{equation}\label{eq:estimator_full}
\hhat^{(K)}(\mby) = \sum_{k=1}^K
\frac{p(k) \calN_{\C}(\mby; \mbA \meanhk, \mbA \covhk \mbA^\herm + \mbSigma)}{\sum_{i=1}^K p(i) \calN_{\C}(\mby; \mbA \meanhi, \mbA \covhi \mbA^\herm + \mbSigma) }
\left(\covhk \mbA^\herm (\mbA \covhk \mbA^\herm + \mbSigma)^{-1} (\mby - \mbA \meanhk) + \meanhk\right)
\end{equation}

A possible application scenario of the gmm estimator would be to use channel samples collected at, for example, the base station of a cellular radio system to construct a site-specific gmm channel estimator.
In an initial (offline) training phase, the channel samples are used to fit a \( K \)-components gmm.
Afterwards, (online) channel estimates are computed via (<ref>).
The formula (<ref>) can also be found, e.g., in <cit.>.
One of the key differences to other work is that we provide a strong motivation to use (<ref>) even if the channel pdf \( f_{\mbh} \) is not a gmm.
<Ref> summarizes both the offline gmm training and the online channel estimation phases.
The necessary number \( M \) of training data depends, e.g., on the number \( K \) of gmm components.
We discuss this in more detail in <Ref>.

GMM Estimator

[1]

training data \( \{ \mbh_m \}_{m=1}^M \), number of components \( K \)
\( (\{ p(k) \}_{k=1}^K, \{ \mbmu_k \}_{k=1}^K, \{ \mbC_k \}_{k=1}^K) \leftarrow \mathrm{EM}(\{ \mbh_m \}_{m=1}^M, K) \)
an EM-algorithm computes all parameters of \( f_{\mbh}^{(K)} \)
.5height .4pt.5
observation \( \mby \), matrix \( \mbA \), noise matrix \( \mbSigma \)
\( \hhat^{(K)} \leftarrow \mbzero \)

\( k = 1 \) to \( K \)
\( p(k \mid \mby) \leftarrow \frac{p(k) \calN_{\C}(\mby; \mbA \meanhk, \mbA \covhk \mbA^\herm + \mbSigma)}{\sum_{i=1}^K p(i) \calN_{\C}(\mby; \mbA \meanhi, \mbA \covhi \mbA^\herm + \mbSigma) } \)
\( \tilde{\mbh} \leftarrow \covhk \mbA^\herm (\mbA \covhk \mbA^\herm + \mbSigma)^{-1} (\mby - \mbA \meanhk) + \meanhk \)
\( \hhat^{(K)} \leftarrow \hhat^{(K)} + p(k \mid \mby) \tilde{\mbh} \)

many quantities in the loop can be precomputed
\( \hhat^{(K)} \)
estimated channel, see (<ref>)

§.§ Computational Complexity

To compute \( \hest(\mby) \) in (<ref>), \( K \) responsibilities \( p(k \mid \mby) \) (<ref>) and \( K \) lmmse formulas (<ref>) need to be evaluated.
Since both the matrix \( \mbA \) and the gmm covariance matrices \( \covhk \) do not change between observations,
the inverse in (<ref>) can be precomputed offline for various snr.
Thus, evaluating (<ref>) online is dominated by matrix-vector multiplications and has a complexity of \( \calO(mN) \).
The responsibilities are calculated by evaluating Gaussian densities, as can be seen from (<ref>).
A Gaussian density with mean \( \mbmu \in \C^{m} \) and covariance matrix \( \mbC \in \C^{m\times m} \) can be written as

\begin{equation}\label{eq:gaussian_density}
\calN_{\C}(\mbx; \mbmu, \mbC) = \frac{\exp(-(\mbx - \mbmu)^\herm \mbC^{-1} (\mbx - \mbmu))}{\pi^{m} \det(\mbC)}.
\end{equation}

Again, since the gmm covariance matrices and mean vectors do not change between observations, the inverse and the determinant of the densities in (<ref>) can be precomputed offline.
Thus, the online evaluation is again dominated by matrix-vector multiplications and has a complexity of \( \calO(m^2) \).
The resulting overall complexity of computing \( \hest \) is \( \calO(K m N) \).

In some cases, as demonstrated in the following subsections, the computational complexity can be reduced by constraining the gmm covariance matrices \( \covhk \) such that corresponding matrix-vector multiplications are accelerated.
In other cases, the number of gmm parameters might be reduced by introducing covariance matrix constraints which can enhance the convergence of the em algorithm, improve the resulting estimation performance, and reduce the required amount of channel samples.
Particular choices for constraints can come from scenario-specific insights.
We demonstrate the feasibility of the following two constraint examples in <Ref>.

§.§.§ Circulant covariance matrices

A first example is a scenario, where the base station employs a ula and where the channel covariance matrix therefore is Toeplitz structured.
For large numbers of antennas, a Toeplitz matrix is well approximated by a circulant matrix <cit.>.
Any circulant matrix \( \mbC \in \C^{N\times N} \) has an eigendecomposition of the form \( \mbC = \mbF^\herm \diag(\mbc) \mbF \)
where \( \mbF \in \C^{N\times N} \) is the dft matrix and where \( \mbc \in \C^N \).
Consequently, thanks to fast Fourier transforms, matrix-vector multiplications involving circulant matrices can be performed in \( \calO(N \log(N)) \) time.
For a large number of antennas, we therefore have a motivation to use circulant covariance matrices \( \covhk = \mbF^\herm \diag(\mbc_k) \mbF \) in the gmm.

This is particularly interesting for a signal model where \( \mbA = \mbI \) and \( \mbSigma = \sigma^2 \mbI = \sigma^2 \mbF \mbF^\herm \).
In this case, the lmmse formula (<ref>) simplifies to

\begin{equation}\label{eq:nlogn_lmmse_formula}
\expec[\mbh \mid \mby, k] = \mbF^\herm \diag(\mbd_k) \mbF (\mby - \meanhk) + \meanhk
\end{equation}

where the \( i \)th entry of the vector \( \mbd_k \) is given by \( [\mbd_k]_i = \frac{[\mbc_k]_i}{[\mbc_k]_i + \sigma^2} \),
such that (<ref>) can be calculated in \( \calO(N \log(N)) \) time.
With a circulant \( \mbC = \mbF^\herm \diag(\mbc) \mbF \), (<ref>) reads as

\begin{multline}\label{eq:nlogn_gaussian_density}
\calN_{\C}(\mbx; \mbmu, \mbF^\herm \diag(\mbc) \mbF) = \\ \frac{ \exp(-(\mbF(\mbx - \mbmu))^\herm \diag(\mbc)^{-1} \mbF (\mbx - \mbmu))}{\pi^n \prod_{i=1}^N [\mbc]_i}.
\end{multline}

A first observation is that this can also be evaluated in \( \calO(N \log(N)) \) time such that computing channel estimates \( \hest \) has a complexity of \( \calO(K N \log(N)) \) if circulant covariance matrices are used in the gmm.
A second observation is that the Gaussian density (<ref>) has a significantly reduced number of parameters: \( N + N \) (mean vector \( \mbmu \) and covariance vector \( \mbc \)) in contrast to \( N + \frac{N(N+1)}{2} \) (mean vector \( \mbmu \) and covariance matrix \( \mbC \)) in (<ref>).
This simplifies the em algorithm iterations of the gmm fitting process and reduces the number of required training channel samples.
For this latter reason, even if \( \mbA \neq \mbI \), one might be interested in employing a gmm with circulant covariance matrices.
The relationship between the number of channel samples and the em algorithm's performance is demonstrated in <Ref>.
In an implementation, instead of constraining the covariance matrices to be circulant,
all channel samples can be Fourier transformed as \( \tilde{\mbh} = \mbF \mbh \) and then the gmm's covariance matrices can be constrained to be diagonal matrices due to the relation \( \mbC = \mbF^\herm \diag(\mbc) \mbF \).

§.§.§ Kronecker covariance matrices

Another example where complexity can be reduced is the mimo signal model from <Ref>.
A well-known assumption for spatial correlation scenarios is that the scattering in the vicinity of the transmitter and of the receiver are independent of each other, cf. <cit.>.
In this case, every channel covariance matrix \( \mbC \) can be decomposed into the Kronecker product of a transmit and receive side spatial covariance matrix: \( \mbC = \mbC_{\text{tx}} \otimes \mbC_{\text{rx}} \).
Here, we have a motivation to use a gmm with Kronecker product covariance matrices \( \covhk = \mbC_{\text{tx},k} \otimes \mbC_{\text{rx},k} \).

To this end, instead of fitting a single gmm using the vectorized channel data of dimension \( N = \Ntx \Nrx \),
one can fit two independent transmit and receive side gmm of dimensions \( \Ntx \) and \( \Nrx \), respectively.
This not only results in lower offline training complexity and in the ability to parallelize, but also in a smaller number of training channel samples needed because the respective gmm have much fewer parameters.
The training channel samples for these low-dimensional gmm are obtained by taking the rows (columns) of the available channel matrices independently in order to fit the transmit (receive) side gmm.
In order to then obtain the full-size covariance matrices \( \mbC_k \),
all combinatorial Kronecker products of transmit and receive side covariance matrices \( \mbC_{\text{tx},i} \) and \( \mbC_{\text{rx},j} \) are computed.
The details are described in the numerical simulations section.

In this example, plugging the Kronecker decomposition \( \covhk = \mbC_{\text{tx},k} \otimes \mbC_{\text{rx},k} \) into the lmmse formula (<ref>) does not lead to an expression that simplifies to a Kronecker product.
This is because the inverse in (<ref>) can generally not be written in terms of a Kronecker product and, thus, full matrix-vector products are necessary.
However, <cit.> explains how (<ref>) can be approximated by means of a Kronecker product in the described setting, which might be interesting if computational complexity of (<ref>) is an issue.
Nonetheless, even if the online computational complexity is not affected, Kronecker gmm covariance matrices can still be beneficial, for instance, if the number of available channel samples is small.
We demonstrate this case in the numerical simulations section.

§.§ Convergence of the Estimator

This subsection uses a universal approximation result of <cit.> to show that if \( f_{\mbh} \) is continuous, then the gmm-based estimator \( \hest \) in (<ref>) can approximate the optimal cme \( \hhat(\mby) \) in (<ref>) arbitrarily well as the number \( K \) of gmm components increases.
The intuition is that if a sequence of pdf \( f_{\mbh}^{(K)} \), which is used in (<ref>), converges to the channel pdf \( f_{\mbh} \),
we can conjecture that also \( f_{\mby}^{(K)} \) from (<ref>) converges to \( f_{\mby} \) and that then \( \hhat^{(K)}(\mby) \) in (<ref>) converges to the cme \( \hhat(\mby) \) in (<ref>).

Recall that the pdf of a complex random vector can be expressed by means of a joint pdf of its real and imaginary parts.
Therefore, this subsection considers real-valued quantities only and the results generalize to the complex-valued setting by considering stacked real and imaginary parts.

To state the main result, we adopt some definitions from <cit.>.
Let \( \calC = \{ f: \R^N \to \R: f \geq 0, \int f(\mbx) d\mbx = 1, f \text{ is contin.} \} \) denote the set of all continuous pdf.
Further, let \( g \) denote the standard Gaussian density and define the class of \( K \)-component location-scale finite Gaussian mixtures as

\begin{equation}
\calM_K = \left\{ h: h(\mbx) = \sum_{k=1}^K c_k \frac{1}{\sigma_k^N} g\left( \frac{\mbx - \mbmu_k}{\sigma_k} \right) \right\}
\end{equation}

with \( \mbmu_k \in \R^N, \sigma_k > 0, c_k \geq 0 \) for all \( k \in \{1, \dots, K\} \) and \( \sum_{k=1}^K c_k = 1 \).
Then, any continuous pdf can be approximated arbitrarily well by means of gmm, as <cit.> states:

\(
\calC_0 = \{ f \in \calC: \forall \varepsilon > 0, \exists \text{ a compact } \calK \subset \R^N \text{ such that } \sup_{\mbx \in \R^N \setminus \calK} | f(\mbx) | < \varepsilon \}
\)
denote the set of all continuous pdf which vanish at infinity.
For any \( f \in \calC_0 \), there exists a sequence \( ( f^{(K)} )_{K=1}^\infty \) with \( f^{(K)} \in \calM_K \) with

\begin{equation}
\lim_{K\to\infty} \|f - f^{(K)}\|_\infty = 0.
\end{equation}

Note that since a pdf is integrable, it always vanishes at infinity such that this is not a constraint for our considerations.
As mentioned, we now work with real quantities \( \mby = \mbA \mbh + \mbn \) where \( \mbn \) is a real Gaussian random vector with mean zero and covariance matrix \( \mbSigma \in \R^{N\times N} \) whose pdf we denote by \( f_{\mbn} \in \calC_0 \).
The pdf of \( \mbh \) is \( f_{\mbh} \) and the pdf of \( \mby \) is \( f_{\mby} \).
The following theorem is proved in Appendix <ref>.

With the notation defined above, let \( \mbA \in \R^{N\times N} \) be invertible and let \( f_{\mbh} \in \calC_0 \) be arbitrary.
Let \( ( f_{\mbh}^{(K)} )_{K=1}^\infty \) be a sequence of pdf in \( \calC_0 \) which converges uniformly to \( f_{\mbh} \).
Then, the estimator

\begin{equation}
\hhat^{(K)}(\mby) = \expec^{(K)}[\mbh^{(K)}\mid\mby] = \int \mbh \frac{f_{\Brv n}(\mby - \mbA\mbh)f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} d \mbh
\end{equation}

approximates the cme

\begin{equation}
\hhat(\mby) = \expec[\mbh\mid\mby] = \int \mbh \frac{f_{\Brv n}(\mby - \mbA\mbh)f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} d \mbh
\end{equation}

in the sense that for any radius \( r > 0 \),

\begin{equation}\label{eq:convergence_estimator}
\lim_{K\to\infty} \| \hhat(\mby) - \hhat^{(K)}(\mby) \| = 0
\end{equation}

holds uniformly for all \( \mby \) in the ball \( \calB_r = \{ \mby \in \R^{N} : \| \mby \| \leq r \} \).
Thus, in particular, by finding a suitable \( r > 0 \), (<ref>) can be seen to hold for any given \( \mby \in \R^N \).

§.§ Discussion of Theorem 2

Both estimators \( \hhat^{(K)} \) and \( \hhat \) are functions which map the current observation \( \mby \) onto corresponding channel estimates \( \hhat^{(K)}(\mby) \) and \( \hhat(\mby) \).
<Ref> proves the pointwise convergence of the function sequence \( ( \hhat^{(K)} )_{K=1}^{\infty} \) to the function \( \hhat \).
In detail, for any observation \( \mby \), the sequence \( ( \hhat^{(K)}(\mby) )_{K=1}^{\infty} \) of channel estimates converges in Euclidean norm to the optimal channel estimate \( \hhat(\mby) \).
To our knowledge, the pointwise convergence of the estimators has not been investigated yet.
The work in <cit.> can be considered to be most related to our result.
Therein, the author assumes that the random vectors \( (\mbh^{(K)}, \mby^{(K)}) \)
converge in distribution to the random vectors \( (\mbh, \mby) \) and the question
is whether the random vectors \( \expec[\mbh^{(K)} \mid \mby^{(K)}] \)
converge in distribution to the random vector \( \expec[\mbh \mid \mby] \).
Here, the condition is still considered to be a random vector
whereas we assume to condition on the current realization, which is given by the observation at the base station.
Thus, <cit.> studies a sequence of random vectors and we study a sequence of functions.
Further, <cit.> studies the convergence in distribution and we study the pointwise convergence.
Therefore, the result in <Ref> is not a consequence of <cit.>.

Next, we discuss the implications of the fact that (<ref>) holds uniformly for all \( \mby \) in a ball \( \calB_r \).
Let \( \calB_r \) be given.
If we want the error \( \| \hhat(\mby) - \hhat^{(K)}(\mby) \| \) to be smaller than a given threshold \( \varepsilon_{\text{thr}} > 0 \), then according to <Ref>, we can find a \( K_r \in \N \) such that \( \| \hhat(\mby) - \hhat^{(K)}(\mby) \| \leq \varepsilon_{\text{thr}} \) holds for all \( \mby \in \calB_r \) and \( K \geq K_r \).
This does not mean that the error is always larger than \( \varepsilon_{\text{thr}} \) for \( \mby \notin \calB_r \).
However, it can be the case, that for certain \( \mby \notin \calB_r \) the number \( K \) of components needs to be strictly larger than \( K_r \) in order for \( \| \hhat(\mby) - \hhat^{(K)}(\mby) \| \) to fall below the threshold.

A requirement of <Ref> is a sequence of pdf which converges uniformly to \( f_{\mbh} \).
By <Ref>, such a sequence always exists if we consider gmm.
However, as argued in <cit.>, there exist other mixtures with universal approximation properties as well.
It is an interesting question whether (<ref>) can be computed in closed form for other mixture models and to see if they for example need fewer components for a satisfying approximation and channel estimation.

<Ref> requires \( \mbA \) to be invertible.
Unfortunately, the proof of <Ref> cannot be conducted as presented if \( \mbA \) is not invertible, which is for example the case when we consider a wide matrix with more columns than rows.
There are multiple challenges involved.
First, the proof makes use of <Ref> which shows that \( \int \| \mbh \| f_{\mbn}(\mby - \mbA \mbh) d\mbh \) is finite for any \( \mby \).
This integral is generally not finite if \( \mbA \) is not invertible (see Appendix <ref>).
Second, for invertible \( \mbA \), we could directly show that the sequence of pdf corresponding to \( \mbA \mbh^{(K)} \) converges uniformly to the pdf of \( \mbA \mbh \).
This will likely not hold for noninvertible \( \mbA \) (see Appendix <ref>).

While a strong statement about the convergence of the estimators does not seem possible with the presented means if \( \mbA \) is not invertible,
we can make the following observation.
Since \( f_{\mbh}^{(K)} \) converges uniformly to \( f_{\mbh} \),
it converges in particular pointwise.
By Scheffe's lemma (e.g., <cit.>), it follows that the random vectors \( \mbh^{(K)} \) converge to \( \mbh \) in distribution.
We also have \( (\mbh^{(K)}, \mbn) \to (\mbh, \mbn) \) in distribution.
If we define the continuous mapping

\begin{equation}
s: \mathbb{R}^{N+m} \to \mathbb{R}^{N+m}, (\mbh, \mbn) \mapsto (\mbh, \mbA \mbh + \mbn)
\end{equation}

then, the continuous mapping theorem (e.g., <cit.>) implies the convergence of \( s(\mbh^{(K)}, \mbn) = (\mbh^{(K)}, \mby^{(K)}) \) to
\( s(\mbh, \mbn) = (\mbh, \mby) \) in distribution.
Given the convergence in distribution of a sequence \( (\mbh^{(K)}, \mby^{(K)}) \) to \( (\mbh, \mby) \),
the author of <cit.> investigates conditions which ensure the convergence of the corresponding conditional expectations \( \expec[\mbh^{(K)} \mid \mby^{(K)}] \) to \( \expec[\mbh \mid \mby] \).
The main result depends without limitation on the distribution of \( f_{\mbh} \) which is not assumed to be known in our setting.
Reciting the main result is beyond the scope of the current paper and we refer the interested reader to <cit.>.

§ CHANNEL MODELS AND RELATED CHANNEL ESTIMATORS

Before we turn to numerical simulations, we introduce the considered channel models and discuss other channel estimation algorithms which we use for comparison.
At this point, it should be noted that the proposed approach does not rely on estimating a (link-based) covariance matrix based on pilot symbols.
Therefore, we do not show comparisons with such approaches as this lies in a different field of applications.
To generate channel samples (for training and testing purposes), we define a scenario like for example a base station which covers a certain \( 120^\circ \) sector.
Afterwards, we choose user positions uniformly at random within the scenario and retrieve their corresponding channels.
The so-obtained set of channel samples can then be used to find estimators for the whole scenario.
These estimators (both the gmm estimator as well as all estimators introduced in the following) are trained/computed/defined once and then tested on the whole scenario without further modification.

§.§ Channel Models

§.§.§ 3GPP

We work with a spatial channel model <cit.> where channels are modeled conditionally Gaussian: \( \mbh \mid \mbdelta \sim \calN(\mbzero, \mbC_{\mbdelta}) \).
The random vector \( \mbdelta \) collects the angles of arrival/departure and path gains of the main propagation clusters between a mobile terminal and the base station.
The main angles are drawn independently and uniformly from the interval \( [0, 2\pi] \) and the path gains are independent zero-mean Gaussians.
The base station employs a ula for both the transmitter and the receiver such that the transmit- and receive-side spatial channel covariance matrix are given by

\begin{equation}
\mbC_{\mbdelta}^{\text{\{rx,tx\}}} = \int_{-\pi}^\pi g^{\text{\{rx,tx\}}}(\theta; \mbdelta) \mba^{\text{\{rx,tx\}}}(\theta) \mba^{\text{\{rx,tx\}}}(\theta)^\herm d \theta.
\end{equation}

\begin{equation}
\mba^{\text{\{rx,tx\}}}(\theta) = [1, e^{j\pi\sin(\theta)}, \dots, e^{j\pi(N_{\text{\{rx,tx\}}}-1)\sin(\theta)}]^\tp
\end{equation}

is the array steering vector for an angle of arrival/departure \( \theta \) and \( g \) is a power density consisting of a sum of weighted Laplace densities whose standard deviations describe the angle spread of the propagation clusters <cit.>.
The full channel covariance matrix is constructed as \( \mbC_{\mbdelta} = \mbC_{\mbdelta}^{\text{tx}} \otimes \mbC_{\mbdelta}^{\text{rx}} \) due to the assumption of independent scattering in the vicinity of transmitter and receiver, see, e.g., <cit.>.
In the simo case, \( \mbC_{\mbdelta} \) degenerates to the receive-side covariance matrix \( \mbC_{\mbdelta}^{\text{rx}} \).
For every channel sample, we generate random angles and path gains, combined in \( \mbdelta \), and then draw the sample as \( \mbh \sim \calN(\mbzero, \mbC_\mbdelta) \).

§.§.§ QuaDRiGa

Version 2.4 of the QuaDRiGa channel simulator <cit.> is used to generate channel samples.
We simulate an urban macrocell scenario at a center frequency of 2.53 GHz.
The base station's height is 25 meters and it covers a \( 120^\circ \) sector.
The minimum and maximum distances between the mobile terminals and the base station are 35 meters and 500 meters, respectively.
In 80% of the cases, the mobile terminals are located indoors at different floor levels, whereas the mobile terminals' height is 1.5 meters in the case of outdoor locations.

QuaDRiGa models the channel of the \( c \)-th carrier and \( t \)-th time symbol as
\( \mbH_{c,t} = \sum_{\ell=1}^{L} \mbG_{\ell} e^{-2\pi j f_c \tau_{\ell,t}} \)
where \( \ell \) is the path number, and the number of multi-path components \( L \) depends on whether there is los, nlos, or o2i propagation: \( L_\text{LOS} = 37 \), \( L_\text{NLOS} = 61 \) or \( L_\text{O2I} = 37 \), cf. <cit.>.
The frequency of the \( c \)-th carrier is denoted by \( f_c \) and the \( \ell\)-th path delay of the \( t \)-th time symbol by \( \tau_{\ell, t} \).
The coefficients matrix \( \mbG_{\ell} \) consists of one complex entry for each antenna pair, which comprises the attenuation of a path, the antenna radiation pattern weighting, and the polarization <cit.>.
As described in the QuaDRiGa manual <cit.>, the generated channels are post-processed to remove the path gain.

For the simulations in <Ref> and <Ref>, we generate single-carrier simo and mimo channels, respectively.
The base station is equipped with a ula with \( \Nrx \) “3GPP-3D” antennas and the mobile terminals employ \( \Ntx \) “omni-directional” antennas.
For the simulations in <Ref>, we consider a SISO system in the spatial domain with \( N_c \) carriers over a bandwidth of 360 kHz and for a time slot with 1 ms duration that is divided into \( N_t \) time symbols.
Each user moves with a certain velocity \( v \) in a random direction.

§.§ State-of-the-Art Channel Estimators

A simple baseline algorithm is the ls channel estimator which computes

\begin{equation}\label{eq:ls}
\hhat_{\text{LS}} = \mbA^\dagger \mby
\end{equation}

using the Moore-Penrose pseudoinverse \( \mbA^\dagger \).
For \( \mbA = \mbI \), there is nothing to compute, and for \( \mbA \neq \mbI \), we have a complexity of \( \calO(mN) \) because the pseudoinverse can be precomputed.
Another immediate estimator consists of first estimating a sample covariance matrix \( \mbC = \frac{1}{M} \sum_{m=1}^M \mbh_m \mbh_m^\herm \) using \( M = 10^5 \) training channel samples drawn uniformly from the whole scenario and then computing lmmse channel estimates:

\begin{equation}\label{eq:sample_cov}
\hhat_{\text{sample cov.}} = \mbC \mbA^\herm (\mbA \mbC \mbA^\herm + \mbSigma)^{-1} \mby.
\end{equation}

Since \( \mbC \) is computed in the offline phase, this estimator also has an online complexity of \( \calO(m N) \).
When we work with the 3GPP channel model from <Ref>, then the true covariance matrix \( \mbC_{\mbdelta} \) for every channel sample is available and we can compute a genie lmmse channel estimate:

\begin{equation}\label{eq:genie_lmmse}
\hhat_{\text{gen. LMMSE}} = \expec[\mbh \mid \mby, \mbdelta] = \mbC_{\mbdelta} \mbA^\herm (\mbA \mbC_{\mbdelta} \mbA^\herm + \mbSigma)^{-1} \mby
\end{equation}

which presents a lower bound for all estimators.
Note that this is not the optimal cme considered in (<ref>) because of the additional genie knowledge of \( \mbdelta \).
The inverse in (<ref>) needs to be computed for every observation because for every observation there is a corresponding \( \mbC_{\mbdelta} \).
Thus, the complexity is \( \calO(m^3 + mN) \).

Many modern channel estimation algorithms focus on cs approaches,
see, e.g., the surveys <cit.>.
In what follows, we therefore consider two cs algorithms <cit.>.
Recent non-cs algorithms often focus on machine learning methods.
For this reason, we also compare with such methods in the numerical simulations <cit.>.

cs approaches assume the channel to be (approximately) sparse: \( \mbh \approx \mbD \mbs \).
Here, \( \mbD \in \C^{N\times L} \) is a dictionary and \( \mbs \in \C^L \) is a sparse vector.
A typical choice for \( \mbD \) is an oversampled dft matrix (e.g., <cit.>).
cs algorithms then assume that \( \mby = \mbA \mbD \mbs + \mbn \) holds and they recover an estimate \( \hat{\mbs} \) of \( \mbs \) and estimate the channel as \( \hhat = \mbD \hat{\mbs} \).
A well-known cs algorithm is omp <cit.>.
omp needs to know the sparsity order.
Since order estimation is a difficult problem, we avoid it via a genie-aided approach: omp gets access to the true channel to choose the optimal sparsity order.
This yields a performance bound for omp.
As explained in <cit.>, every iteration has a complexity of \( \calO(N \log(N)) \)
and the number of iterations is equal to the genie-determined sparsity order.
Another algorithm, which we use for comparison, is amp <cit.>, which does not need to know the sparsity order.
The computational complexity of amp is not analyzed in <cit.>.
Since it is an iterative algorithm, it depends without limitation on the number of iterations.
We set the number of iterations per channel estimate to 100 in the simulations.

A cnn-based channel estimator was derived in <cit.> for the simo signal model (cf. <Ref>).
In <cit.>, the cnn estimator has been generalized to the mimo signal model (cf. <Ref>).
In the simo case, we use the cnn estimator as described in <cit.>.
The activation function is the rectified linear unit and we use the input transform based on the \( 2N \times 2N \) Fourier matrix, cf. <cit.>.
In the mimo case, we use the cnn estimator as described in <cit.>
where again the activation function is the rectified linear unit.
In all cases, the cnn is trained on samples corresponding to the channel model on which it is tested later.
The computational complexity is \( \calO(N \log(N)) \) <cit.>.

The concept of a cae was introduced in <cit.> and adapted for wideband channel estimation (cf. <Ref>) in <cit.>.
The cae is an autoencoder where the encoder is replaced by a concrete selection layer that selects the \( N_p \ll N_c N_t \) most informative features of the \( N_c N_t \)-dimensional input.
This corresponds to designing the pilot matrix by selecting the \( N_p \) pilot positions.
The decoder can then be used to perform channel estimation.
During training, noisy channels are given as input, such that the decoder of the cae performs denoising and reconstruction of the full-dimensional channels.
Hence, a new cae needs to be trained for every different snr.
In our simulations, in contrast to <cit.>, no further denoising networks are applied after the cae.
The complexity of cae is not analyzed in <cit.>.

The authors in <cit.> propose a deep cnn approach for 2D wideband channel estimation.
The estimator, called ChannelNet, consists of a combination of an image super-resolution and an image restoration network.
Thus, the networks perform interpolation and denoising of the low-dimensional observations with respect to the high-dimensional channel matrix.
The super-resolution network consists of three 2D convolution layers, whereas the image restoration network consists of 20 2D convolution layers.
Here, too, we train snr-specific networks.
The complexity of ChannelNet is not analyzed in <cit.>.

§ NUMERICAL SIMULATIONS

In all simulations, a normalized mse (nMSE) is used as performance measure.
Specifically, we generate \( T = 10^4 \) \( N \)-dimensional test channel samples \( \{ \mbh_t \}_{t=1}^T \), obtain corresponding channel estimates \( \hhat_t \), and
define \( \text{nMSE} = \frac{1}{NT} \sum_{t=1}^T \| \mbh_t - \hhat_t \|^2. \)
The noise covariance matrix is \( \mbSigma = \sigma^2 \mbI \).
The test samples are normalized such that \( \expec[\|\mbh\|^2] = N \) holds which
allows us to define an snr as \( \frac{1}{\sigma^2} \).
For training purposes, we generate \( M = 10^5 \) channel samples unless stated otherwise.
The number of training samples is always chosen large enough such that increasing \( M \) does not lead to a performance improvement during the testing phase.
The generated channel samples stem from one of the scenarios described in <Ref>.

As explained in <Ref> (see also <Ref> there), to obtain the gmm-based estimator,
we fit one gmm using the available training data via an em algorithm.
Afterwards, inverses which appear in (<ref>) are precomputed for every snr.
In contrast to this approach, the introduced neural network-based estimators need to be newly trained for every snr at which we evaluate them.
This includes searching for suitable hyperparameters for every snr.

§.§ SIMO

xlabel=\( N \), number of antennas,
ymode = log,
ymin= 1e-2,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 3,
legend entries=
legend style=at=(0.0,1.0), anchor=south west,
legend cell align = left,

table[x=antennas, y=amp, col sep=comma]

table[x=antennas, y=cnn_fft2x_relu_non_hier_False, col sep=comma]

table[x=antennas, y=genie lmmse, col sep=comma]

table[x=antennas, y=global cov, col sep=comma]

table[x=antennas, y=gmm full all, col sep=comma]

table[x=antennas, y=gmm diag all, col sep=comma]

table[x=antennas, y=ls, col sep=comma]

table[x=antennas, y=Genie_OMP, col sep=comma]

SIMO signal model (<Ref>) and 3GPP channel model (<Ref>) with one propagation cluster at 10 dB snr. The performance of the circulant gmm estimator (“circ. GMM”, <ref>, <Ref>) is shown too. In both cases, \( K = 128 \) components are used.

xlabel=SNR [dB],
ymode = log,
ymin= 1e-5,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 2,
legend entries=
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x=SNR, y=amp, col sep=comma]

table[x=SNR, y=cnn_fft2x_relu_non_hier_False, col sep=comma]

table[x=SNR, y=genie lmmse, col sep=comma]

table[x=SNR, y=global cov, col sep=comma]

table[x=SNR, y=gmm full all, col sep=comma]

table[x=SNR, y=gmm diag all, col sep=comma]

table[x=SNR, y=ls, col sep=comma]

table[x=SNR, y=Genie_OMP, col sep=comma]

SIMO signal model (<Ref>) and 3GPP channel model (<Ref>) with one propagation cluster and \( N = 128 \) antennas. The performance of the circulant gmm estimator (“circ. GMM”, <ref>, <Ref>) is shown too. In both cases, \( K = 128 \) components are used.

xlabel=SNR [dB],
ymode = log,
ymin= 5e-5,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 2,
legend entries=
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x=SNR, y=amp, col sep=comma]

table[x=SNR, y=cnn_fft2x_relu_non_hier_False, col sep=comma]

table[x=SNR, y=genie lmmse, col sep=comma]

table[x=SNR, y=global cov, col sep=comma]

table[x=SNR, y=gmm full all, col sep=comma]

table[x=SNR, y=gmm diag all, col sep=comma]

table[x=SNR, y=ls, col sep=comma]

table[x=SNR, y=Genie_OMP, col sep=comma]

SIMO signal model (<Ref>) and 3GPP channel model (<Ref>) with three propagation clusters and \( N = 128 \) antennas. The performance of the circulant gmm estimator (“circ. GMM”, <ref>, <Ref>) is shown too. In both cases, \( K = 128 \) components are used.

xlabel=SNR [dB],
ymode = log,
ymin= 1e-4,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 2,
legend entries=
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x=SNR, y=amp, col sep=comma]

table[x=SNR, y=cnn_fft2x_relu_non_hier_False, col sep=comma]

table[x=SNR, y=global cov, col sep=comma]

table[x=SNR, y=gmm full all, col sep=comma]

table[x=SNR, y=gmm diag all, col sep=comma]

table[x=SNR, y=ls, col sep=comma]

table[x=SNR, y=Genie_OMP, col sep=comma]

SIMO signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( N = 128 \) antennas. The performance of the circulant gmm estimator (“circ. GMM”, <ref>, <Ref>) is shown too. In both cases, \( K = 128 \) components are used.

<Ref> show channel estimation results for the simo signal model from <Ref>.
The cs algorithms omp and amp are used with oversampled dft dictionaries that have \( L = 4 N \) and \( L = 2 N \) columns, respectively,
because these parameters yielded the best results.
Unless stated otherwise, the gmm fitting process uses \( 19 \cdot 10^4 \) training data.

In <Ref>, we consider the 3GPP channel model from <Ref> with one propagation cluster.
The snr is 10 dB and the number \( N \) of antennas is varied.
It is interesting to see that the gmm-based estimator performs almost as well as the genie lmmse estimator.
As the number \( N \) of antennas increases, the cnn estimator starts to outperform the genie omp estimator.
The reason for this is that the assumptions under which the cnn estimator was derived in <cit.> are better fulfilled for a larger number of antennas.
Generally, the relative performance between all estimators hardly differs with different numbers of antennas.
This is an observation we have made in all our experiments.
For this reason, in what follows, the number of antennas is fixed at \( N = 128 \).

In <Ref>, we consider again the 3GPP channel model from <Ref> with one propagation cluster.
For almost all snr values, the gmm-based estimator performs almost as well as the genie lmmse estimator.
In the mid-snr range, the two cs algorithms are approximately equally good.
In <Ref>, we have three propagation clusters.
A first observation is the strong performance of the cnn estimator in the mid-snr range.
Note that we can generally not expect any estimator to reach the genie lmmse curve because it has more channel knowledge (the true covariance matrix for every sample).
In the higher snr-range, the gmm-based estimator is the only algorithm still outperforming ls estimation.

In <Ref>, we concentrate on the QuaDRiGa channel model described in <Ref> where the channel covariance matrices and therefore the genie lmmse curve are no longer available.
Here, the two cs algorithms behave not as similarly as they did in the previous experiments.
Additionally, their performance is not as convincing.
A reason might be that the channels now are not sparse enough.
The cnn estimator shows again a good performance and overall the gmm-based estimator can compete with it or is better.

In addition to the gmm estimator (<ref>),
<Ref> display the performance of the reduced-complexity gmm estimator which uses circulant covariance matrices as described in <Ref>.
As expected, the estimator's performance suffers but it is still comparable to the other algorithms.
<Ref> is particularly interesting where there is not much difference between the full- and low-complexity gmm estimators.

<Ref> shows the behavior of the gmm-based estimator for different numbers of components, \( K \).
We consider an snr of 10 dB and the same situation as in <Ref>:
The 3GPP channel model (cf. <Ref>) with three propagation clusters and \( N = 128 \) antennas.
In addition to \( K \), also the number of training data used to fit the gmm is varied.
Since the number of parameters of a gmm increases when \( K \) is increased,
more training data is necessary for a good fit.
This effect is clearly visible in <Ref>.
Overall, as long as the number of training data is high enough (\( M \geq 300 \) in the figure),
increasing \( K \) leads to an mse improvement,
which is in accordance with <Ref>.

Note that we cannot expect the gmm estimator to converge to the genie lmmse estimator (<ref>) (which is displayed in <Ref>).
The genie lmmse estimator has more knowledge (namely the true channel covariance matrix \( \mbC_{\mbdelta} \)) and is therefore not the cme, \( \hhat = \expec[\mbh \mid \mby] \), which we want to approximate in <Ref>.
The cme cannot be computed in closed form in the considered scenario which is the main motivation to study the gmm estimator in the first place.

<Ref> also shows the gmm estimator's behavior for different numbers of components, \( K \),
but now for the QuaDRiGa channel model (cf. <Ref>).
For all displayed snr, a saturation can be observed as \( K \) is increased,
and already a moderate number of components can lead to a satisfactory estimation performance.
Altogether, a smaller \( K \) tends to be sufficient for higher snr.
Generally, a suitable number of components needs to be determined based on the training data size as well as on the desired estimator complexity.

xlabel=\( K \), number of components,
ymode = log,
yticklabels=\( 4 \cdot 10^{-2} \), \( 7 \cdot 10^{-2} \), \( 10^{-1} \),
ylabel=Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 3,
legend entries=
\( M = 10 \),
\( M = 25 \),
\( M = 50 \),
\( M = 100 \),
\( M = 300 \),
\( M = 500 \),
legend cell align = left,
legend style=at=(0.0,1.0), anchor=south west,

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain10]

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain25]

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain50]

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain100]

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain300]

table[x=components, y=snr_10, col sep=comma, discard if not=ntrain500]

SIMO signal model (<Ref>) and 3GPP channel model (<Ref>) with three propagation clusters and \( N = 128 \) antennas. The SNR is 10 dB. The GMM estimator is trained using \( M \cdot 10^3 \) samples.

xlabel=\( K \), number of components,
ymode = log,
ylabel = Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 2,
legend entries=
\( \text{SNR} = -10 \text{ dB} \),
\( \text{SNR} = 0 \text{ dB} \),
\( \text{SNR} = 10 \text{ dB} \),
\( \text{SNR} = 20 \text{ dB} \),
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x=K, y=nMSE, col sep=comma, discard if not=SNR-10]

table[x=K, y=nMSE, col sep=comma, discard if not=SNR0]

table[x=K, y=nMSE, col sep=comma, discard if not=SNR10]

table[x=K, y=nMSE, col sep=comma, discard if not=SNR20]

SIMO signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( N = 128 \) antennas. The training data size is \( 19 \cdot 10^4 \).

§.§ MIMO

xlabel=SNR \( [\operatorname{dB}] \),
ymode = log,
ymin= 8*1e-3,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 1,
legend entries=
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x= SNR, y=gmm_full, col sep=comma]

table[x= SNR, y=gmm_kron, col sep=comma]

table[x= SNR, y=Genie_Aided, col sep=comma]

table[x= SNR, y=Genie_OMP, col sep=comma]

table[x= SNR, y=cnn_fft_relu, col sep=comma]

table[x=SNR, y=ls, col sep=comma]

MIMO signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( (\Nrx, \Ntx) = (32, 4) \). The GMM estimator (“GMM”, <ref>) uses \( K = 32 \) components  and the Kronecker GMM estimator (“Kron. GMM”, <ref>, <Ref>) uses \( K = 4 \times 8 \) components.

xlabel= \( M \), number of training samples,
ymode = log,
ymin= 1e-2,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 1,
legend entries=
legend style=at=(1.0,1.0), anchor=north east,
legend cell align = left,

table[x=num, y=gmm_full, col sep=comma]

table[x=num, y=gmm_kron, col sep=comma]

MIMO signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( (\Nrx, \Ntx) = (32, 4) \). The SNR is 10 dB. The Kronecker product GMM estimator (“Kron. GMM”, <ref>, <Ref>) has fewer parameters than the normal GMM estimator. Both have \( K = 32 \) components.

For the mimo simulations whose signal model is described in <Ref>,
we use a scaled dft pilot matrix \( \mbP \).
Here, we are mainly interested in comparing the gmm estimator in (<ref>) to a gmm estimator which uses Kronecker product covariance matrices as described in
We generate \( M = 10^5 \) training channel samples \( \{ \mbH_i \}_{i=1}^M \) and use them in two different ways.
Either, we use the (vectorized) channel matrices \( \mbH_i \) directly to fit a single \( K = 32 \) components gmm.
Or, we view all columns and all rows of the data as separate data sets and fit two gmm:
One transmit side gmm with \( \Ktx = 4 \) components using the rows,
one receive side gmm with \( \Krx = 8 \) components using the columns.
Afterwards, we combine the two gmm to a single full-size gmm with \( K = \Ktx \Krx = 32 \) components for the whole data set of channel matrices.
To this end, the \( K \) full-size means and covariances are directly calculated by combinatorial computation of the Kronecker products of the transmit- and receive-side gmm components.
The corresponding mixing coefficients can be computed by fixing the means and covariances and performing a single E-step (cf. <cit.>) in the em algorithm to obtain them.

<Ref> displays simulation results.
A first observation is that the two gmm estimators as well as the cnn estimator perform very similarly, with minor differences in the lower and higher snr-regimes.
Further, these three estimators outperform the sample covariance matrix-based estimator and the genie-aided omp algorithm, which uses a Kronecker product of two two-times oversampled dft matrices as dictionary.
A second observation is that there is almost no difference between the gmm estimator with or without Kronecker product covariance matrices.
This is insofar surprising as the normal gmm consists of \( K = 32 \) covariance matrices of dimension \( N \times N \) with \( N = 32 \cdot 4 = 128 \)
which means that it has \( K \frac{N(N+1)}{2} = 264192 \) covariance parameters,
whereas in contrast, the Kronecker gmm has only \( \Krx \frac{\Nrx(\Nrx+1)}{2} + \Ktx \frac{\Ntx(\Ntx+1)}{2} = 4224 + 40 = 4264 \) covariance parameters.
Since the Kronecker gmm has significantly fewer parameters,
it should require a smaller number \( M \) of training data.
This is confirmed in <Ref> where the two estimators are compared at an snr of 10 dB for varying \( M \).
The curves intersect between \( M = 10^4 \) and \( M = 10^5 \).

§.§ Wideband

In this section, we show numerical results for the wideband signal model described in Section <ref> and with the QuaDRiGa simulation setup described in Section <ref>.
We chose a typical 5G frame structure as defined in <cit.> with \( N_c = 24 \) carriers over a bandwidth of 360 kHz with 15 kHz carrier spacing and with \( N_t = 14 \) time symbols over a time slot with 1 ms duration.
The number of pilot symbols is \( N_p = 50 \), which means that 50 out of the \( 24 \times 14 = 336 \) resource elements are occupied with pilot tones.
We compare the gmm-based estimator with the lmmse estimator based on the sample covariance matrix (see (<ref>)), with the cae approach, and with the ChannelNet estimator (see <Ref>).
The training data for each approach consist of \( M = 10^5 \) channel realizations and corresponding observations from the pilot positions.

In <Ref>, we depict mse results over the snr for a scenario where every user moves at \( v=3 \) km/h speed, in which the block-type pilot arrangement has shown the best results.
One can observe superior performance of the cae over the sample covariance lmmse estimator and the ChannelNet, which may be due to the fact that the cae optimizes the pilot pattern.
However, the gmm approach is able to outperform all baseline algorithms, where the performance gap increases with increasing snr.
Further, the impact of more components (from \( K=8 \) to \( K=128 \)) is visible and results in better performance for all snr values.

<Ref> shows the same setup but now each user's velocity is randomly chosen between 0 - 300 km/h (in both training and testing sets), which makes the estimation more challenging and the lattice-type pilot arrangement superior.
First, the performance gap to the sample covariance lmmse estimator increases, which is a result of the more diverse setting that cannot be captured well by a single covariance matrix.
Further, the number of gmm components seems to play a more important role.
The gmm with \( K=128 \) components is still able to compete with the ChannelNet and cae approaches.

Finally, <Ref> shows the mse behavior for different numbers of components, \( K \), and for different amounts of training data used to fit the gmm.
Similar to the results in <Ref>, also the wideband case shows an improving performance for increasing numbers of components.
This at least provides numerical evaluation of the convergence of the gmm estimator for noninvertible observation matrices.

xlabel=SNR $[\operatorname{dB}]$,
ymode = log,
ymin= 1e-5,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 1,
legend entries=
(block); \( K = 8 \),
(block); \( K = 128 \),
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x= SNR, y=MMSE_block, col sep=comma]

table[x= SNR, y=channel_net, col sep=comma]

table[x= SNR, y=CAE, col sep=comma]

table[x= SNR, y=gmm_from_y_block, col sep=comma]

table[x= SNR, y=gmm_from_y_block, col sep=comma]

Wideband signal model (<Ref>) and QuaDRiGa channel model (<Ref>)
with \( N_p = 50 \) pilot tones for \( N_c = 24 \) carriers and \( N_t = 14 \) time symbols for \( v = 3 \) km/h.

xlabel=SNR $[\operatorname{dB}]$,
ymode = log,
ymin= 1e-4,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 1,
legend entries =
(lattice); \( K = 8 \),
(lattice); \( K = 128 \),
legend style=at=(0.0,0.0), anchor=south west,
legend cell align = left,

table[x= SNR, y=MMSE_lattice, col sep=comma]

table[x= SNR, y=channel_net_lattice, col sep=comma]

table[x= SNR, y=CAE, col sep=comma]

table[x= SNR, y=gmm_from_y_lattice, col sep=comma]

table[x= SNR, y=gmm_from_y_lattice, col sep=comma]

Wideband signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( N_p=50 \) pilot tones for \( N_c=24 \) carriers and \( N_t=14 \) time symbols for \( v \in [0,300] \) km/h.

xlabel=\( K \), number of components,
ylabel= Normalized MSE,
ylabel shift = 0.0cm,
grid = both,
legend columns = 1,
legend entries=
\( M = 10 \),
\( M = 50 \),
\( M = 100 \),
\( M = 200 \),
legend style=at=(1.0,1.0), anchor=north east,
legend cell align = left,

table[x=comp, y=gmm_from_y_block, col sep=comma]

table[x=comp, y=gmm_from_y_block, col sep=comma]

table[x=comp, y=gmm_from_y_block, col sep=comma]

table[x=comp, y=gmm_from_y_block, col sep=comma]

Wideband signal model (<Ref>) and QuaDRiGa channel model (<Ref>) with \( N_p = 50 \) pilot tones for \( N_c = 24 \) carriers and \( N_t = 14 \) time symbols for \( v = 3 \) km/h and block-type pilot pattern at an SNR of 20 dB. The GMM estimator is trained using \( M \cdot 10^3 \) samples.

§ CONCLUSION AND OUTLOOK

We studied the behavior of a gmm channel estimator when the number of gmm components is increased.
The gmm estimator is tailored for a particular
communications environment.
This is the case because the underlying gmm
uses training data which stem, for example,
from the coverage area of a base station.
Thereafter, the gmm estimator can be employed for channel estimation in this environment.

In case of an invertible observation matrix, we proved the convergence of the gmm estimator to the optimal cme.
Notably, the proof only assumes that a sequence of pdf exists which converges uniformly to the true channel pdf.
An example of such a sequence is given by gmm but no properties unique to gmm are used in the proof.
In particular, any sequence of cme which is based on a uniformly convergent sequence of pdf also converges to the optimal cme.
This invites the study of new channel estimators based on other universal approximators.

While the theoretical results are of an asymptotic nature,
we analyzed a number of practically relevant settings (simo, mimo, and wideband) in numerical simulations.
There, already with a moderate number of gmm components, the proposed estimator outperforms various state-of-the-art approaches in all depicted scenarios.
In particular, increasing the number of gmm components lead to a performance improvement even for noninvertible observation matrices.
Additionally, we demonstrated how scenario-specific insights can be used to reduce the gmm estimator's complexity.

§ PROOF OF THEOREM 2

The proof of <Ref> makes use of and is presented after the following lemma.

For an arbitrary \( \mby \in \R^N \), it holds

\begin{multline*}
\int \| \mbh \| f_{\Brv n}(\mby - \mbA\mbh) d \mbh
\leq \\ \sqrt{\det(\mbA^{-1}\mbA^{-\tp})} \sqrt{\| \mbA^{-1} \mby \|^2 + \mathrm{trace}(\mbA^{-1}\mbSigma\mbA^{-\tp})}.
\end{multline*}

Recall that \( f_{\Brv n} \) denotes a Gaussian pdf with mean zero and covariance matrix \( \mbSigma \).
We have

\begin{align}
&f_{\mbn}(\mby - \mbA\mbh) = \frac{\exp(-\frac{1}{2}(\mby - \mbA\mbh)^\tp \mbSigma^{-1} (\mby - \mbA\mbh))}{\sqrt{(2\pi)^N \det(\mbSigma)}}
\\&= \frac{\exp(-\frac{1}{2}(\mbA^{-1}\mby - \mbh)^\tp \mbA^\tp \mbSigma^{-1} \mbA (\mbA^{-1}\mby - \mbh))}{\sqrt{(2\pi)^N \det(\mbSigma)}}
\\&= \frac{\exp(-\frac{1}{2}(\mbh - \mbA^{-1}\mby)^\tp (\mbA^{-1} \mbSigma \mbA^{-\tp})^{-1} (\mbh - \mbA^{-1}\mby))}{\sqrt{(2\pi)^N \det(\mbSigma)}}.
\end{align}

Therefore, \( (\sqrt{\det(\mbA^{-1} \mbA^{-\tp})})^{-1} f_{\mbn}(\mby - \mbA\mbh) =: \tilde{f}(\mbh) \) is a Gaussian pdf with mean vector \( \tilde{\mbmu} := \mbA^{-1} \mby \) and covariance matrix \( \tilde{\mbSigma} := \mbA^{-1} \mbSigma \mbA^{-\tp} \).
We are interested in computing

\begin{align}
\nonumber
&\int \| \mbh \| f_{\Brv n}(\mby - \mbA\mbh) d \mbh
\\&= \sqrt{\det(\mbA^{-1} \mbA^{-\tp})} \int \| \mbh \| \frac{f_{\mbn}(\mby - \mbA\mbh)}{\sqrt{\det(\mbA^{-1} \mbA^{-\tp})}} d \mbh
\\&= \sqrt{\det(\mbA^{-1} \mbA^{-\tp})} \int \| \mbh \| \tilde{f}(\mbh) d \mbh
\label{eq:to_be_bounded}
\end{align}

The integral in (<ref>) computes the expected value of the norm of a random vector with the pdf \( \tilde{f} \).
Let \( \Brv w \) be a standard Gaussian random vector (with mean \( \mbzero \in \R^N \) and covariance matrix \( \mbI \in \R^{N\times N} \)) and let \( \tilde{\mbSigma}^{\frac{1}{2}} \) be a square root of the covariance matrix \( \tilde{\mbSigma} = \tilde{\mbSigma}^{\frac{\tp}{2}} \tilde{\mbSigma}^{\frac{1}{2}} \).
Then, \( \tilde{\mbmu} + \tilde{\mbSigma}^{\frac{1}{2}} \Brv w \) is a random vector with the pdf \( \tilde{f} \).
Thus, we can express the integral in (<ref>) as:

\begin{equation}
\int \| \mbh \| \tilde{f}(\mbh) d \mbh
= \expec_{\Brv w}[\| \tilde{\mbmu} + \tilde{\mbSigma}^{\frac{1}{2}} \Brv w \|]
\label{eq:expectation_wrt_w}
\end{equation}

where we compute the expected value with respect to the standard Gaussian random vector \( \Brv w \).
Jensen's inequality yields:

\begin{align}
&\left( \expec_{\Brv w}[\| \tilde{\mbmu} + \tilde{\mbSigma}^{\frac{1}{2}} \Brv w \|] \right)^2
\leq \expec_{\Brv w}[\| \tilde{\mbmu} + \tilde{\mbSigma}^{\frac{1}{2}} \Brv w \|^2]
\\&= \| \tilde{\mbmu} \|^2 + 2 \tilde{\mbmu}^\tp \tilde{\mbSigma}^{\frac{1}{2}} \expec_{\Brv w}[\Brv w] + \expec_{\Brv w}[\mathrm{trace}(\Brv w^\tp \tilde{\mbSigma}^{\frac{\tp}{2}} \tilde{\mbSigma}^{\frac{1}{2}} \Brv w)]
\\&= \| \tilde{\mbmu} \|^2 + \mathrm{trace}(\tilde{\mbSigma} \expec_{\Brv w}[\Brv w \Brv w^\tp])
= \| \tilde{\mbmu} \|^2 + \mathrm{trace}(\tilde{\mbSigma})
\\&= \| \mbA^{-1} \mby \|^2 + \mathrm{trace}(\mbA^{-1}\mbSigma\mbA^{-\tp}).
\end{align}

Now, we take the square root on both sides and plug the result into (<ref>).
Then, we use this in (<ref>) to conclude.

First, we show that the convergence of \( f_{\Brv h}^{(K)} \) to \( f_{\Brv h} \) implies the convergence of \( f_{\Brv y}^{(K)} \) to \( f_{\Brv y} \).
The pdf of \( \Brv x := \mbA \Brv h \) is

\begin{equation}
f_{\Brv x}(\mbx) = \frac{1}{|\det(\mbA)|} f_{\Brv h}(\mbA^{-1} \mbx)
\end{equation}

because \( \mbA \) is invertible.
Since the random vector \( \Brv y = \Brv x + \Brv n \) is a sum of two stochastically independent random vectors, its pdf can be computed via convolution:

\begin{equation}\label{eq:convolution}
f_{\Brv y}(\mby) = \int f_{\Brv n}(\mbs) f_{\Brv x}(\mby - \mbs) d\mbs.
\end{equation}

Similarly, \( f_{\Brv y}^{(K)} \) is obtained by replacing \( f_{\Brv x} \) with \( f_{\Brv x}^{(K)} \) in (<ref>).
For later reference, note that because \( f_{\Brv n} \) is positive (\( f_{\Brv n}(\mbs) > 0 \) for all \( \mbs \in \R^N \)) and \(f_{\Brv x}\) as well as \(f_{\Brv x}^{(K)}\) are continuous pdf, the convolution results \( f_{\Brv y}^{(K)} \) and \( f_{\Brv y} \) are positive, too.
We have

\begin{align}
&| f_{\Brv y}(\mby) - f_{\Brv y}^{(K)}(\mby) |
\\&= \left| \int f_{\Brv n}(\mbs) \left(f_{\Brv x}(\mby - \mbs) - f_{\Brv x}^{(K)}(\mby - \mbs)\right) d\mbs \right|
\\&\leq \int \left|f_{\Brv n}(\mbs) \frac{f_{\Brv h}(\mbA^{-1}(\mby - \mbs)) - f_{\Brv h}^{(K)}(\mbA^{-1}(\mby - \mbs))}{|\det(\mbA)|} \right| d\mbs
\\&\leq \frac{\| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty}{|\det(\mbA)|} \int | f_{\Brv n}(\mbs) | d\mbs = \frac{\| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty}{|\det(\mbA)|}.
\end{align}

The last integral is equal to one because \( f_{\Brv n} \) is a pdf.
Since \( \lim_{K\to\infty} \| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty = 0 \) holds by assumption, we have

\begin{equation}\label{eq:convergence_fyK}
\lim_{K\to\infty} \| f_{\Brv y} - f_{\Brv y}^{(K)} \|_\infty = 0.
\end{equation}

To show (<ref>), let \( \mby \in \calB_r \) be arbitrary.
With (<ref>) and (<ref>) in mind, we find the following upper bound:

\begin{multline}
\| \hhat - \hhat^{(K)} \|
\leq \int \| \mbh \| | f_{\Brv n}(\mby - \mbA\mbh) | \left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right| d \mbh
\\\leq \sup_{\mbh\in\R^N} \left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right| \int \| \mbh \| f_{\Brv n}(\mby - \mbA\mbh) d \mbh.
\end{multline}

The last integral is independent of \( K \) and by <Ref>, it is finite for any \( \mby \in \R^N \).
It is in particular bounded by

\begin{equation*}
\sqrt{\det(\mbA^{-1}\mbA^{-\tp})} \sqrt{\| \mbA^{-1} \|^2 r^2 + \mathrm{trace}(\mbA^{-1}\mbSigma\mbA^{-\tp})} < \infty
\end{equation*}

for all \( \mby \in \calB_r \).
Hence, as soon as

\begin{equation}\label{eq:convergence_frac_xy}
\lim_{K\to\infty} \sup_{\mbh\in\R^N} \left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right| = 0, \quad \forall \mby \in \calB_r
\end{equation}

is shown, (<ref>) is confirmed.
To prove (<ref>), we write

\begin{equation}\label{eq:frac_common_denominator}
\left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right|
\left| \frac{f_{\Brv h}(\mbh) f_{\Brv y}^{(K)}(\mby) - f_{\Brv y}(\mby) f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}(\mby) f_{\Brv y}^{(K)}(\mby)} \right|
\end{equation}

for an arbitrary \( \mbh \in \R^N \).
Now, we add \( 0 = f_{\Brv y}^{(K)}(\mby) f_{\Brv h}^{(K)}(\mbh) - f_{\Brv y}^{(K)}(\mby) f_{\Brv h}^{(K)}(\mbh) \) in the numerator on the right-hand side and apply the triangle inequality to get

\begin{align}
&\left| \frac{f_{\Brv h}(\mbh) f_{\Brv y}^{(K)}(\mby) - f_{\Brv y}(\mby) f_{\Brv h}^{(K)}(\mbh) }{f_{\Brv y}(\mby) f_{\Brv y}^{(K)}(\mby)} \right|
\\&\leq\nonumber
\frac{\left| \left(f_{\Brv h}(\mbh) - f_{\Brv h}^{(K)}(\mbh)\right) f_{\Brv y}^{(K)}(\mby) \right|}{f_{\Brv y}(\mby) f_{\Brv y}^{(K)}(\mby)}
\\&+
\frac{\left| \left(f_{\Brv y}^{(K)}(\mby) - f_{\Brv y}(\mby)\right) f_{\Brv h}^{(K)}(\mbh) \right|}{f_{\Brv y}(\mby) f_{\Brv y}^{(K)}(\mby)}
\\&\leq \label{eq:frac_upper_bound}
\frac{\| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty \| f_{\Brv y}^{(K)} \|_\infty
+ \| f_{\Brv y}^{(K)} - f_{\Brv y} \|_\infty \| f_{\Brv h}^{(K)} \|_\infty}{f_{\Brv y}(\mby) f_{\Brv y}^{(K)}(\mby)}.
\end{align}

By the compactness of \( \calB_r \) and continuity of \( f_{\Brv y} \), there exists a \( \mby_{\text{min}} \in \calB_r \) at which \( f_{\Brv y} \) attains a minimum value \( f_{\Brv y}(\mby_{\text{min}}) > 0 \) over \( \calB_r \).
Due to the uniform convergence (<ref>), there exists an index \( N_1 \in \N \) such that \( | f_{\Brv y}(\mby) - f_{\Brv y}^{(K)}(\mby)| \leq \frac{1}{2} f_{\Brv y}(\mby_{\text{min}}) \) holds for all \( K \geq N_1 \) and for all \( \mby \in \calB_r \).
The reverse triangle inequality then shows that \( f_{\Brv y}^{(K)}(\mby) \geq f_{\Brv y}(\mby) - | f_{\Brv y}(\mby) - f_{\Brv y}^{(K)}(\mby) | \geq f_{\Brv y}(\mby_{\text{min}}) - \frac{1}{2} f_{\Brv y}(\mby_{\text{min}}) = \frac{1}{2} f_{\Brv y}(\mby_{\text{min}}) \) is true.
Hence, with \( M_1 := \frac{1}{2} f_{\Brv y}(\mby_{\text{min}}) > 0 \), the inequality

\begin{equation}\label{eq:N1}
f_{\Brv y}^{(K)}(\mby) \geq M_1
\end{equation}

holds for all \( \mby \in \calB_r \) and for all \( K \geq N_1 \).
Further, since \( \| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty \to 0 \) and \( \| f_{\Brv h} \|_\infty < \infty \), there exist \( M_2 > 0 \) and \( N_2 \in \N \) such that

\begin{equation}\label{eq:N2}
\| f_{\Brv h}^{(K)} \|_\infty \leq M_2 \quad\text{for all}\quad K \geq N_2.
\end{equation}

Analogously, there exist \( M_3 > 0 \) and \( N_3 \in \N \) such that

\begin{equation}\label{eq:N3}
\| f_{\Brv y}^{(K)} \|_\infty \leq M_3 \quad\text{for all}\quad K \geq N_3.
\end{equation}

Let \( \varepsilon > 0 \) be arbitrary.
Due to \( \| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty \to 0 \), there exists an index \( N_4 \geq \max\{N_1, N_3\} \) such that

\begin{equation}\label{eq:N4}
\| f_{\Brv h} - f_{\Brv h}^{(K)} \|_\infty \leq \frac{f_{\Brv y}(\mby_{\text{min}}) M_1}{2 M_3} \varepsilon \quad\text{for all}\quad K \geq N_4.
\end{equation}

Similarly, there exists an index \( N_5 \geq \max\{N_1, N_2\} \) with

\begin{equation}\label{eq:N5}
\| f_{\Brv y} - f_{\Brv y}^{(K)} \|_\infty \leq \frac{f_{\Brv y}(\mby_{\text{min}}) M_1}{2 M_2} \varepsilon \quad\text{for all}\quad K \geq N_5.
\end{equation}

We can use the last five inequalities to bound (<ref>).
To this end, \( f_{\Brv y}(\mby_{\text{min}}) \) and (<ref>) provide bounds on the terms in the denominator, (<ref>) and (<ref>) bound the first summand in the numerator, and (<ref>) and (<ref>) bound the second summand in the numerator.
In total, this yields an upper bound on (<ref>):

\begin{align}
\nonumber
\left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right|
&\leq \frac{f_{\Brv y}(\mby_{\text{min}}) M_1}{2 M_3} \varepsilon \cdot \frac{M_3}{f_{\Brv y}(\mby_{\text{min}}) M_1} \\&+ \frac{f_{\Brv y}(\mby_{\text{min}}) M_1}{2 M_2} \varepsilon \cdot \frac{M_2}{f_{\Brv y}(\mby_{\text{min}}) M_1}.
\end{align}

for all \( K \geq \max\{N_4, N_5\} \) and for all \( \mby \in \calB_r \).
We conclude

\begin{equation}
\sup_{\mbh\in\R^N} \left| \frac{f_{\Brv h}(\mbh)}{f_{\Brv y}(\mby)} - \frac{f_{\Brv h}^{(K)}(\mbh)}{f_{\Brv y}^{(K)}(\mby)} \right|
\leq \varepsilon
\end{equation}

for all \( K \geq \max\{N_4, N_5\} \) and all \( \mby \in \calB_r \), and because \( \varepsilon \) was arbitrary, (<ref>) is confirmed, which finishes the proof.

§ INTEGRAL FOR NONINVERTIBLE MATRICES

To see why \( \int \| \mbh \| f_{\mbn}(\mby - \mbA \mbh) d\mbh \) might not be finite, consider the matrix \( \mbA = [\mbI, \mbzero] \in \R^{m\times N} \) where \( \mbI \in \R^{m\times m} \) is the identity matrix and the remaining matrix elements are zero.
Let us write \( \mbh = [\mbh_m^\tp, \mbh_{N-m}^\tp]^\tp \in \R^m \times \R^{N-m} \).
Define the set \( \calC = \{ \mbh \in \R^N \mid \| \mbh_m \| \geq 1, \| \mbh_{N-m} \| \geq 1 \} \) where the norm of both sub-vectors \( \mbh_m \) and \( \mbh_{N-m} \) is at least one such that we always have \( \| \mbh \| \geq 1 \) on \( \calC \).
We can now compute:

\begin{align}
&\int_{\R^N} \| \mbh \| f_{\mbn}(\mby - \mbA \mbh) d\mbh
\geq \int_{\calC} f_{\mbn}(\mby - \mbA \mbh) d\mbh
\\&= \int_{\| \mbh_{N-m} \| \geq 1} \int_{\| \mbh_m \| \geq 1} f_{\mbn}(\mby - \mbh_m) d\mbh_m d\mbh_{N-m}.
\end{align}

Since \( f_n \) is an \( m \)-dimensional Gaussian pdf, the inner integral is equal to some constant \( c \) with \( 0 < c < 1 \) and it follows that \( \int \| \mbh \| f_{\mbn}(\mby - \mbA \mbh) d\mbh \) is not finite.

§ ON THE UNIFORM CONVERGENCE

Let us express the pdf of \( \mbA \mbh \).
If \( \mbA \in \R^{m\times N} \) is a wide matrix with full rank \( m \), we can assume that the first \( m \) columns are linearly independent (otherwise we introduce a permutation matrix).
This allows us to partition \( \mbA = [\mbA_i, \mbA_n] \) into an invertible part \( \mbA_i \in \R^{m\times m} \)
and a noninvertible part \( \mbA_n \in \R^{m\times N-m} \).
With a corresponding partitioning of \( \mbh = [\mbh_i^\tp, \mbh_n^\tp]^\tp \in \R^m \times \R^{N-m} \),
we have \( \mbx = \mbA \mbh = \mbA_i \mbh_i + \mbA_n \mbh_n \) and
we can define an invertible mapping \( t: \R^N \to \R^N \):

\begin{align}
t: (\mbh_i, \mbh_n) &\mapsto (\mbA_i \mbh_i + \mbA_n \mbh_n, \mbh_n) = (\mbx, \mbx')
\\t^{-1}: (\mbx, \mbx') &\mapsto (\mbA_i^{-1}(\mbx - \mbA_n \mbx'), \mbx') = (\mbh_i, \mbh_n).
\end{align}

We can now compute the joint density \( f_{\mbx,\mbx'} \) with the usual transformation formula:

\begin{equation}
f_{\mbx, \mbx'}(\mbx, \mbx') =
\frac{f_{\mbh}(t^{-1}(\mbx, \mbx'))}{\left|\det\left(\frac{\partial t}{\partial \mbh}(t^{-1}(\mbx, \mbx'))\right)\right|}.
\end{equation}

Together with \( \left|\det\left(\frac{\partial t}{\partial \mbh}(t^{-1}(\mbx, \mbx'))\right)\right| = | \det(\mbA_i) | \),
we can express the pdf of \( \mbx = \mbA \mbh \) via marginalization:

\begin{equation}
= \int_{\R^{N-m}} f_{\mbx,\mbx'}(\mbx,\mbx') d \mbx'
= \int_{\R^{N-m}} \frac{f_{\mbh}(t^{-1}(\mbx, \mbx'))}{| \det(\mbA_i) |} d \mbx'.
\end{equation}

Analogously, one obtains \( f_{\mbx}^{(K)} \) for \( \mbx^{(K)} = \mbA \mbh^{(K)} \).

\begin{multline}
|f_{\mbx}^{(K)}(\mbx) - f_{\mbx}(\mbx)| = \frac{1}{| \det(\mbA_i) |} \times \\ \left| \int_{\R^{N-m}} \left(  f_{\mbh}^{(K)}(t^{-1}(\mbx, \mbx')) - f_{\mbh}(t^{-1}(\mbx, \mbx')) \right) d \mbx' \right|,
\end{multline}

we can conjecture that due to the integral over \( \R^{N-m} \) the uniform convergence of \( f_{\mbh}^{(K)} \) to \( f_{\mbh} \) alone is generally not sufficient to infer the uniform convergence of \( f_{\mbx}^{(K)} \) to \( f_{\mbx} \).