Interference Cancellation Information Geometry Approach for Massive MIMO Channel Estimation

From information geometry theory, we have that the manifold of complex Gaussian distribution is a dually flat manifold, which has an affine coordinate system and a dual affine coordinate system. The two coordinate systems are also called natural parameters and expectation parameters in the literature. In the following, we describe the two affine coordinate systems in detail.

The Gaussian distributions belong to the exponential family of distributions [30]. Let $\bm{\theta},\bm{\Theta}$ be the natural parameter of a complex Gaussian distribution of a random vector $\mathbf{x}\in\mathbb{C}^{N\times 1}$, then the PDF $p(\mathbf{x};\bm{\theta},\bm{\Theta})$ is defined as

where $\psi(\bm{\theta},\bm{\Theta})$ is the normalization factor, which is called free energy function and given by

Let $\mathcal{M}=\{p(\mathbf{x};\bm{\theta},\bm{\Theta})\}$ be the manifold of multivariate complex Gaussian distributions and $\bm{\theta},\bm{\Theta}$ is an coordinate system. The free energy function $\psi(\bm{\theta},\bm{\Theta})$ is a convex function of $\bm{\theta},\bm{\Theta}$ and introduces an affine flat structure, which means the $\bm{\theta},\bm{\Theta}$ is an affine coordinate system and each coordinate axis of $\bm{\theta},\bm{\Theta}$ is a straight line. Furthermore, the Bregman divergence [31] from $\bm{\theta},\bm{\Theta}$ to $\bm{\theta}^{\prime},\bm{\Theta}^{\prime}$ derived from $\psi(\bm{\theta},\bm{\Theta})$ is the same as the Kullback–Leibler (KL) divergence from $p(\mathbf{x};\bm{\theta}^{\prime},\bm{\Theta}^{\prime})$ to $p(\mathbf{x};\bm{\theta},\bm{\Theta})$.

The dual coordinate system of $\mathcal{M}$ is obtained by the Legendre transformation [32], i.e., the gradients of the free energy function $\psi(\bm{\theta},\bm{\Theta})$. From (2), we can obtain the dual coordinate $\bm{\mu},\mathbf{M}$ as

where $\bm{\Sigma}$ is the covariance matrix of $\mathbf{x}$. The dual coordinate is the combination of the first and second order moments of $\mathbf{x}$ and is also called the expectation parameter. The dual function of $\psi$ is given by

which is the negative entropy of the PDF $p(\mathbf{x};\bm{\theta},\bm{\Theta})$. By using the dual coordinate system, we have that

where $c$ is a constant. The dual function $\phi$ is a convex function of $\bm{\mu},\mathbf{M}$ and induces the dual affine flat structure. The Bregman divergence from $\bm{\mu},\mathbf{M}$ to $\bm{\mu}^{\prime},\mathbf{M}^{\prime}$ derived from $\phi$ is the KL divergence from $p(\mathbf{x};\bm{\mu},\mathbf{M})$ to $p(\mathbf{x};\bm{\mu}^{\prime},\mathbf{M}^{\prime})$.

The transformations from expectation parameters to natural parameters can also be obtained from the Legendre transformation as

where $\bm{\Sigma}=\mathbf{M}-\bm{\mu}\bm{\mu}^{H}$ is the covariance matrix and is used here for brevity. By using the expectation parameter, the PDF can also be written in the familiar form as

where $\psi(\bm{\mu},\bm{\Sigma})=\log(\pi^{N})+\log\det(\bm{\Sigma})+\bm{\mu}^{H}\bm{\Sigma}^{-1}\bm{\mu}$.

After introducing the affine and dual affine coordinate systems, we present the definitions of $e$-flat and $m$-projection, which are very important when describing the information geometry approach for channel estimation.

A submanifold $\mathcal{M}_{1}\subset\mathcal{M}$ is called $e$-flat if it has a linear constraint in the affine coordinate $\bm{\theta},\bm{\Theta}$. The term $m$-flat can be similarly defined by using the dual affine coordinate $\bm{\mu},\mathbf{M}$. An example of $e$-flat submanifold is the manifold of independent complex Gaussian distributions, defined as

where $\bm{\Theta}_{0}\in\mathbb{R}^{N\times N}$ is a diagonal matrix. Since $\bm{\Theta}_{0}$ are diagonal matrices, the expectation parameter is very easy to obtain as

The projection to an $e$-flat manifold is called $m$-projection because the projection can be realized linearly in the dual affine coordinate system. Let $p(\mathbf{x};\bm{\theta}_{0},\bm{\Theta}_{0})$ and $p(\mathbf{x};\bm{\theta}_{1},\bm{\Theta}_{1})$ be two points in the manifold $\mathcal{M}$, refered as $P_{0}$ and $P_{1}$, and $P_{0}\in\mathcal{M}_{0}$ and $P_{1}\notin\mathcal{M}_{0}$.

The $m$-projection is unique and minimizes the KL divergence. Specifically, the $m$-projection from $P_{1}$ to $\mathcal{M}_{0}$ is defined as

By using the affine and dual affine coordinate systems, the $m$-projection is easy to obtain. First, rewrite the KL divergence as [17]

Then, the minimum can be obtained from the first-order optimal condition

where $\mathbf{I}\odot{\mathbf{M}}_{1}$ is obtained because $\bm{\Theta}_{0}$ is a diagnal matrix. Thus, we have $\bm{\mu}_{0}=\bm{\mu}_{1}$ and $\mathbf{I}\odot{\mathbf{M}}_{0}=\mathbf{I}\odot{\mathbf{M}}_{1}$ for the projection point, and $\bm{\theta}_{1}^{0},\bm{\Theta}_{1}^{0}$ are their dual coordiantes. The projection is simple in the dual affine coordinate system. Let $P_{1}^{0}$ be the projection point, the dual straight line connecting $P_{1}$ and $P_{1}^{0}$ is the shortest one among those dual straight lines from $P_{1}$ to $M_{0}$.

In massive MIMO channel estimation, a general received signal model is given by

where $\mathbf{A}\in\mathbb{C}^{M\times N}$ is the deterministic measurement matrix, $\mathbf{h}$ is a random Gaussian vector distributed as $\mathbf{h}\sim\mathcal{CN}(\mathbf{0},\mathbf{D})$, and $\mathbf{z}$ is the complex Gaussian noise vector distributed as $\mathbf{h}\sim\mathcal{CN}(\mathbf{0},\sigma_{z}^{2}\mathbf{I})$.

For the received signal model in (13), the posterior distribution can be written as

According to (1), the natural parameters of $p(\mathbf{h}|\mathbf{y})$ are

whereas the dual affine coordinate can be obtained from (3a) as

The dual affine coordinate $\bm{\mu}$ is the posterior mean of $\mathbf{h}$, and thus is also the MMSE estimation. For massive MIMO systems, the complexity of $\bm{\mu}$ is often too high due to the inversion of the large dimensional matrix. Thus, one of the most important problems for massive MIMO is to derive low-complexity channel estimation methods.

From Section II-B, we know that if we can $m$-project $p(\mathbf{h};\bm{\theta},\bm{\Theta})$ onto the $e$-flat submanifold ${\cal{M}}_{0}$, then $\bm{\mu}$ is equal to the mean at the projection point, which is easy to obtain. However, the $m$-projection still involves the matrix inversion of the large dimensional matrix, and thus can not provide a low-complexity solution.

Information geometry provides other low-complexity ways to find a point in ${\cal{M}}_{0}$ whose dual coordinate is approximation or the same as that of the $m$-projection point instead of using the $m$-projection. The IGA algorithm proposed in [11] is a specific algorithm derived based on information geometry, but its complexity can be further reduced. To extend the information geometry approach to derive new low-complexity algorithms, we provide a framework of information geometry for massive MIMO channel estimation.

We call the submanifold ${\cal{M}}_{0}$ the target manifold since we want to find a target point in it. Since the target point can not be obtained directly by using the $m$-projection, the auxiliary manifolds and PDFs are needed to find a way to approximate the $m$-projection. The natural parameters or affine coordinates of original posterior PDF are $\bm{\theta}_{or}=\sigma_{z}^{-2}\mathbf{A}^{H}\mathbf{y}$, $\bm{\Theta}_{or}=-(\sigma_{z}^{-2}\mathbf{A}^{H}\mathbf{A}+\mathbf{D}^{-1})$. With the auxiliary manifolds, the process of the information geometry framework for channel estimation is summarized below:

1.  (1)

   The natural parameters $\bm{\theta}_{or}$ and $\bm{\Theta}_{or}$ are split to construct $Q$ auxiliary manifolds of PDFs and one target manifold of PDFs;

2.  (2)

   Initialize the auxiliary points and the target point in the auxiliary manifolds and target manifold, respectively;

3.  (3)

   Calculate the $m$-projections of the auxiliary points to the target manifold and compute the beliefs in the affine coordinate system;

4.  (4)

   Update the natural parameters of the auxiliary and target points;

5.  (5)

   Repeat (3) and (4) until the algorithm converges or fixed iterations, output the mean and variance of the target point.

With the framework, the $m$-projection of the original point to the target manifold is approximated by the $m$-projections from the auxiliary points to the target manifold. To make the approximation well enough, two important conditions described in the following subsections are needed.

In the general information geometry framework, the most important thing is to define the auxiliary points and manifolds. These definitions depend on the way of splitting natural parameter and determine what specific algorithm can be derived.

The split of $\bm{\theta}_{or}$ and $\bm{\Theta}_{or}$ into $Q$ items is given as

where the setting of $\mathbf{b}_{q}$, $\mathbf{C}_{q}$, $Q$, and $\bm{\Lambda}_{c}$ depends on specific algorithms, and $\bm{\Lambda}_{c}$ is usually a diagonal matrix. Let $\mathbf{a}_{q}=[\mathbf{A}^{H}]_{:,q}$ be the $q$-th column of $\mathbf{A}^{H}$. In the IGA algorithm proposed in [11], the split $\mathbf{b}_{q}=\mathbf{a}_{q}y_{q}$, $\mathbf{C}_{q}=\mathbf{a}_{q}\mathbf{a}_{q}^{H}$, $Q=M$, and $\bm{\Lambda}_{c}=\mathbf{D}^{-1}$ is used.

Based on the split, we define $Q$ auxiliary points or PDFs. Let $\bm{\Lambda}_{q}$ be a diagonal matrix. The natural parameter of the $q$-th auxiliary point $p(\mathbf{h};{\bm{\theta}}_{q},{\bm{\Theta}}_{q})$ is defined by

where $\bm{\lambda}_{q}$ and $\bm{\Lambda}_{q}$ are variables in the natural parameter. The corresponding auxiliary PDF and manifold are then given as

and

These auxiliary manifolds $\mathcal{M}_{q}$s are parallel to each other when $\mathbf{C}_{q}$s are not diagonal since the points in different auxiliary manifolds never intersect.

The target manifold is still the manifold of independent complex Gaussian distributions $\mathcal{M}_{0}=\left\{p(\mathbf{h};{\bm{\theta}}_{0},{\bm{\Theta}}_{0})\right\}$. To be consistent with the auxiliary points, the natural parameter of the target point is defined as

where ${\bm{\Lambda}}_{0}$ is diagonal. The corresponding PDF $p(\mathbf{h};{\bm{\theta}}_{0},{\bm{\Theta}}_{0})$ is

The target manifold is also parallel to all the auxiliary manifolds.

Now, we introduce the first condition of the information geometry framework as

It means the original point, the target point and the auxiliary points are on a hyperplane in the affine coordinate system. This is called the $e$-condition since it is a linear condition in the coordinate system $\bm{\theta},\bm{\Theta}$ as shown in Fig. 1. Because of the split, the $e$-condition always holds if

Thus, the above formula is also the $e$-condition. The $e$-condition makes sure the target point and auxiliary points are related to the original point, which is important to finally get an approximation of the $m$-projection point.

The $e$-condition only states the relation between the natural parameters of the original point, the auxiliary points, and the target point, but what we need is part of the expectation parameters of the original point. To obtain an approximation of the $m$-projection point, we need another condition, i.e, the $m$-condition.

The $m$-condition is named because it is a linear condition in the dual affine coordinate system $\bm{\mu},\mathbf{M}$, and is given by

From Section II-B, it means all the $m$-projection points of the auxiliary points are the same and equal to the target point. By combining the $m$-condition with the $e$-condition, a good approximation of the $m$-projection point on the target manifold of the original point can be obtained.

From (3a), the expectation parameters ${\bm{\mu}}_{q},{\mathbf{M}}_{q}$ of auxiliary points can be obtained as

Then, the expectation parameters of the $m$-projection points on the target manifold $\mathcal{M}_{0}$ satisfies

as shown in Section II-B, and further we have the covariance of the $m$-projection ${\bm{\Sigma}}_{q}^{0}=\mathbf{I}\odot{\bm{\Sigma}}_{q}$. The natural parameters of the $m$-projection points can be obtained as

To make both the $e$-condition and the $m$-condition hold, the auxiliary points need to exchange beliefs. Define ${\bm{\lambda}}_{q}^{0}$ and ${\bm{\Lambda}}_{q}^{0}$ to make ${\bm{\theta}}_{q}^{0}={\bm{\lambda}}_{q}^{0}$ and ${\bm{\Theta}}_{q}^{0}=-({\bm{\Lambda}}_{q}^{0}+\bm{\Lambda}_{c})$ hold. The beliefs are defined as

Then, the natural parameters of the $m$-projection points can be expressed as

By comparing them with the natural parameters of auxiliary point $p(\mathbf{h};{\bm{\theta}}_{q},{\bm{\Theta}}_{q})$, it can be observed that $\bm{\xi}_{q}$, $\bm{\Xi}_{q}$ are approximations of $\mathbf{b}_{q}$, $\mathbf{C}_{q}$ in ${\bm{\theta}}_{q}$ and ${\bm{\Theta}}_{q}$, respectively, where $\bm{\Xi}_{q}$ is also diagonal.

After defining the beliefs, the iterative update of natural parameters of the target point and auxiliary points are constructed as

and

From the above two equations, it is observed that the $e$-condition always holds. An information geometry based algorithm is obtained by iteratively calculating (28a), (29a), (31a) and (32a). When the algorithm converges, it is easy to obtain that

which means all the $m$-projection points are equal to the target point, and is equivalent to the $m$-condition. In conclusion, both the $e$-condition and the $m$-condition are satisfied when the algorithm converges.

Finally, the output of this algorithm is the mean and covariance of the target point

It is regarded as the approximated mean and covariance of the marginal PDF corresponding to the original PDF. When the algorithm does not converge, damping can be introduced in the updating of beliefs to ensure the convergence of the algorithm without changing its equilibrium. Although the process of the $m$-projection in (28a) still involves the matrix inversion, i.e., $({\bm{\Lambda}}_{q}+\mathbf{C}_{q}+\bm{\Lambda}_{c})^{-1}$, it can be implemented with low complexity by properly setting of $\mathbf{C}_{q}$ since ${\bm{\Lambda}}_{q}$ and $\bm{\Lambda}_{c}$ are diagonal matrices. For example, in the IGA algorithm proposed in [11], the matrix $\mathbf{C}_{q}$ is a rank-$1$ matrix.

This iterative process and the stabilization point can be explained geometrically. Define the $m$-flat manifold $\mathcal{M}^{\star}$ and the $e$-flat manifold $\mathcal{E}^{\star}$ as

where $c_{q}$ are the positive coefficients. The geometric interpretation of the $m$-condition is given by Fig. 2. The $e$-flat auxiliary and target manifolds are perpendicular to the $m$-flat manifold $\mathcal{M}^{\star}$ under the metric induced by the KL divergence.

From the $e$-condition and the $m$-condition shown in Figs. 1 and 2, it can be seen that when the algorithm converges, all the auxiliary points, the target points and the original points are also on the same $e$-flat manifold $\mathcal{E}^{\star}$. Meanwhile, all the auxiliary points and the target point are on the same $m$-flat manifold $\mathcal{M}^{\star}$. The $m$-condition makes sure the $m$-projection points of the auxiliary points and the target point are the same. The $e$-condition relates the auxiliary points and the target point to the original point. By combining these two conditions, the original point is close to the manifold $\mathcal{M}^{\star}$. Theorem 11.8 of [17] proves that $\mathcal{M}^{\star}$ contains the original point when the factor graph is acyclic, i.e., which means the exact mean of the original distribution is obtained. However, this property is not guaranteed when the factor graph is cyclic, which is the common case for the channel estimation problem. Fortunately, for the channel estimation problem, it is usually able to prove the mean $\hat{\bm{\mu}}_{0}^{\star}$ is equal to the mean $\bm{\mu}_{or}$ when the algorithm converges as shown in [11].

In this subsection, we provide a new way of splitting the natural parameters based on a modified MMSE form. Then, the corresponding auxiliary manifolds are constructed.

To derive a new low-complexity IGA algorithm, we want to make each auxiliary manifold focus on the computation of one element of $\mathbf{h}$. According to the framework of information geometry methods, the $n$-th auxiliary PDF can be defined as

where ${\bm{\theta}}_{n}={\bm{\lambda}}_{n}+\mathbf{b}_{n}$ , ${\bm{\Theta}}_{n}={\bm{\Lambda}}_{n}+\mathbf{C}_{n}$ and $\bm{\Lambda}_{c}$ is set to be a zero matrix.

Let $\mathbf{K}=\sigma_{z}^{-2}\mathbf{A}^{H}\mathbf{A}$. To make the $n$-th auxiliary PDF only compute the information of the $n$-th element of $\mathbf{h}$, a natural idea is to set $\mathbf{C}_{n}$ and $\mathbf{b}_{n}$ as

where $\mathbf{a}_{n}$ is the $n$-th column of $\mathbf{A}$, $\bar{\mathbf{k}}_{n}$ is the vector obtained by deleting the $n$-element $k_{nn}$ from the $n$-column of $\mathbf{K}$, $c_{n}=\sigma_{z}^{-2}\mathbf{a}_{n}^{H}\mathbf{a}_{n}+d_{n}^{-1}$, and $\mathbf{P}_{1n}\in\mathbb{R}^{N\times N}$ is the ordering matrix obtained by extracting $n$-th row of $\mathbf{I}_{N}$ and put it in the first row. The matrix $\mathbf{P}_{1n}$ has the property that $\mathbf{P}_{1n}\mathbf{P}_{1n}^{H}=\mathbf{P}_{1n}^{H}\mathbf{P}_{1n}=\mathbf{I}_{N}$. By left-multiplying $\mathbf{P}_{1n}$ or right-multiplying $\mathbf{P}_{1n}^{H}$, one can extract the $n$-th row or column of a given matrix and put it in the first row or column.

However, there is a flaw in this setup. The items associated with elements other than $h_{n}$ in $\mathbf{h}$ are missing in the auxiliary function, and thus might not form a Gaussian-distributed PDF since ${\bm{\Lambda}}_{n}+\mathbf{C}_{n}$ are not always positive semidefinite. To overcome this issue, we present the following theorem.

Based on Theorem 1, we can use the following PDF

to compute the original mean. With the modified MMSE form, we propose a new way of splitting natural parameter as

where $\mathbf{b}_{n}$ and $\mathbf{C}_{n}$ are defined as

The natural parameter of the $n$-th auxiliary PDF is defined as

where $\bm{\lambda}_{-n},\bm{\Lambda}_{-n}$ are given by

The subscript ${-n}$ instead of $n$ is used because we impose more constraints on $\bm{\lambda}_{n},\bm{\Lambda}_{n}$, and the $n$-th element has been set to zero. The auxiliary PDF $p(\mathbf{h};\bm{\theta}_{n},\bm{\Theta}_{n})$ can be constructed as

and the corresponding auxiliary manifolds are $\mathcal{M}_{n}=\left\{p(\mathbf{h};\bm{\theta}_{n},\bm{\Theta}_{n})\right\},\forall n\in\mathbb{Z}_{N}^{+}$.

The natural parameter of the target point is defined as $\bm{\theta}_{0}=\bm{\lambda}$, $\bm{\Theta}_{0}=-\bm{\Lambda}$, where

The target PDF becomes

and the target manifold is still $\mathcal{M}_{0}=\left\{p(\mathbf{h};\bm{\theta}_{0},\bm{\Theta}_{0})\right\}$. It is easy to check the $e$-condition always holds due to

After constructing the auxiliary manifolds, the estimation problem can be solved by applying the information geometry framework of section III. For convenience, we define $\bm{\lambda}_{n},\bm{\Lambda}_{n}$ as

The following theorem gives the beliefs $\bm{\xi}_{n}$, $\bm{\Xi}_{n}$ corresponding to the $n$-th auxiliary point at the $t$-th iteration.

Before proceeding to the derivation of IC-IGA, we present more insights about Theorem 2. Let $\bar{\mathbf{h}}_{n}$ and $s_{n}$ be defined as $\bar{\mathbf{h}}_{n}=(h_{1},\cdots,h_{n-1},h_{n+1},\cdots,h_{N})^{T}$ and

we have that

where $Z_{n}$ is the normalization constant. It means this part of $p(\mathbf{h};\bm{\theta}_{n},\bm{\Theta}_{n})$ can be viewed as a conditional PDF of $h_{n}$ as $p(h_{n}|h_{1},\cdots,h_{n-1},h_{n+1},\cdots,h_{N},\mathbf{y})$. Furthermore, the remain part of $p(\mathbf{h};\bm{\theta}_{n},\bm{\Theta}_{n})$ can be viewed as

The corresponding PDF $p_{n}(\mathbf{h};\bm{\theta}_{n},\bm{\Theta}_{n})$ can be rewritten as