SiML: Sieved Maximum Likelihood for Array Signal Processing

Consider an array of $L$ sensors with positions $\{\bm{p}_{1},\ldots,\bm{p}_{L}\}\subset\mathbb{R}^{3}$. Assuming the emitting sources are in the far-field [3] of the sensor array, they can be thought of as lying on the unit sphere $\mathbb{S}^{2}$. To allow for arbitrary numbers of complex sources, we consider a notional continuous source field covering the entire sphere, with an associated amplitude function, describing for each direction $\bm{r}\in\mathbb{S}^{2}$ the emission strength of the source field. In practice, source amplitudes fluctuate randomly [1, 14], and the amplitude function can be modelled as a complex random function $\mathcal{S}=\{S(\bm{r}):\Omega\rightarrow\mathbb{C},\bm{r}\in\mathbb{S}^{2}\}$, where $\Omega$ is some probability space. More precisely, we assume $\mathcal{S}$ to be a Gaussian random function [21], i.e., that all its finite marginals have distribution:

$$\left[S(\bm{r}_{1}),\cdots,S(\bm{r}_{n})\right]\stackrel{{\scriptstyle d}}{{\sim}}\mathbb{C}\mathcal{N}_{n}(0,\mathcal{K}_{\mathcal{I}_{n}}),\,\forall\,\mathcal{I}_{n}\subset\mathbb{S}^{2},\,\forall\,n\in\mathbb{N},$$

where $\mathbb{C}\mathcal{N}_{n}$ denotes the $n$-variate centrally symmetric, complex Gaussian distribution [27, 28], $\mathcal{I}_{n}:=\{\bm{r}_{1},\ldots,\bm{r}_{n}\}$ and $\mathcal{K}_{\mathcal{I}_{n}}\in\mathbb{C}^{n\times n}$ is some valid covariance matrix depending on the set ${\mathcal{I}_{n}}$.

From the Huygens-Fresnel principle [29], exciting the source field with a narrowband waveform of wavelength $\lambda\in\mathbb{R}$ results in a diffracted wavefront, which, after travelling through an assumed homogeneous medium is recorded by the sensor array. In a far-field context, the Fraunhofer equation [29, 3] permits this wavefront at each sensor position $\bm{p}_{i}\in\mathbb{R}^{3}$ to be approximated by:

$\displaystyle Y(\bm{p}_{i})$

$\displaystyle=\int_{\mathbb{S}^{2}}\,S(\bm{r})\,\mbox{exp}\left(-j\frac{2\pi}{\lambda}\langle\bm{r},\bm{p}_{i}\rangle\right)\,d\bm{r}\;+\;n_{i},$

(1)

where $i=1,\ldots,L,$ and $\bm{n}=\left[n_{1},\ldots,n_{L}\right]$ is an additive white noise term capturing the inaccuracies in measurement of each sensor, distributed as [1]

$$\bm{n}\stackrel{{\scriptstyle d}}{{\sim}}\mathbb{C}\mathcal{N}_{L}(\bm{0},\sigma I_{L}),\quad\sigma>0.$$

Noise across sensors is assumed to be identically and independently distributed and independent of the random amplitude function $\mathcal{S}$.

We assume that every realisation, or sample function [21], $s_{\omega}:~\mathbb{S}^{2}\rightarrow\mathbb{C}$ of $\mathcal{S}$ is an element of some Hilbert space $\mathcal{H}=\mathcal{L}^{2}(\mathbb{S}^{2},\mathbb{C})$ of finite-energy functions, and thus Eq. 1 can be written as:

$\displaystyle Y(\bm{p}_{i})$

$\displaystyle=\langle S,\phi_{i}\rangle\;+\;n_{i},\quad i=1,\ldots,L,$

where $\phi_{i}(\bm{r}):=\mbox{exp}\left(j2\pi\langle\bm{r},\bm{p}_{i}\rangle/\lambda\right),$ which in turn can be re-written more compactly using an analysis operator [19] $\Phi^{\ast}:\mathcal{H}\rightarrow\mathbb{C}^{L}$, mapping an element of $\mathcal{H}$ to a finite number $L$ of measurements:

$$\bm{Y}=\left[\begin{array}[]{c}Y(\bm{p}_{1})\\ \vdots\\ Y(\bm{p}_{L})\end{array}\right]=\left[\begin{array}[]{c}\langle S,\phi_{1}\rangle\\ \vdots\\ \langle S,\phi_{L}\rangle\end{array}\right]+\left[\begin{array}[]{c}n_{1}\\ \vdots\\ n_{L}\end{array}\right]=\Phi^{\ast}S+\bm{n}.$$

We call $\Phi^{\ast}$ the sampling operator [19] associated with the sensor array. As the sum of two independent centred complex Gaussian random vectors, the vector of measurements $\bm{Y}$ is also a centred complex Gaussian random vector with covariance matrix $\Sigma\in\mathbb{C}^{L\times L}:$

$\displaystyle(\Sigma)_{ij}=\iint_{\mathbb{S}^{2}\times\mathbb{S}^{2}}\kappa(\bm{r},\bm{\rho})\,\phi^{\ast}_{i}(\bm{r})\phi_{j}(\bm{\rho})\,d\bm{r}d\bm{\rho}\;+\;\sigma\,\delta_{ij},$

(2)

where $i,j=1,\ldots,L,$ $\delta_{ij}$ denotes the Kronecker delta and $\kappa:\mathbb{S}^{2}\times\mathbb{S}^{2}\rightarrow\mathbb{C}$ is the covariance kernel [21] of $\mathcal{S}$:

$$\kappa(\bm{r},\bm{\rho}):=\mathbb{E}\left[S(\bm{r})S^{\ast}(\bm{\rho})\right],\;\bm{r},\bm{\rho}\in\mathbb{S}^{2}.$$

Introducing the associated covariance operator $\mathcal{T}_{\kappa}:\mathcal{H}\rightarrow\mathcal{H}$:

$\displaystyle(\mathcal{T}_{\kappa}f)(\bm{r}):=\int_{\mathbb{S}^{2}}\kappa(\bm{r},\bm{\rho})f(\bm{\rho})d\bm{\rho},\quad f\in\mathcal{H},\,\bm{r}\in\mathbb{S}^{2},$

we can again reformulate Eq. 2 in terms of the sampling operator $\Phi^{\ast}$ and its adjoint, called the synthesis operator [19], $\Phi:\mathbb{C}^{L}\rightarrow\mathcal{H}$:

$$\Sigma=\Phi^{\ast}\mathcal{T}_{\kappa}\Phi\;+\;\sigma I_{L}.$$

(3)

By analogy with the finite dimensional case [30], it is customary to write $\kappa=\mbox{vec}(\mathcal{T}_{\kappa}),$ where the $\mbox{vec}(\cdot)$ operator maps an infinite-dimensional linear operator onto its associated kernel representation. Because of the Gaussianity assumption, the covariance kernel $\kappa$ (or equivalently the covariance operator $\mathcal{T}_{\kappa}$) completely determines the distribution of $\mathcal{S}$. Our goal is hence to leverage Eq. 3 in order to form an estimate of $\kappa$ from the covariance matrix $\Sigma$ of the instrument recordings. Often the source field is assumed to be spatially uncorrelated, in which case the random function $\mathcal{S}$ is Gaussian white noise [21] and $\kappa$ becomes diagonal. The diagonal part of $\kappa$

$$I(\bm{r}):=\kappa(\bm{r},\bm{r}),\quad\bm{r}\in\mathbb{S}^{2},$$

is called the intensity function of the source field, of crucial interest in many array signal processing applications.

In practice of course, the covariance matrix $\Sigma$ needs to be estimated from a finite number of i.i.d. observations of $\bm{Y}$, say $N$. Typically, the maximum likelihood estimate of $\Sigma$ is formed by $\hat{\Sigma}=\frac{1}{N}\sum_{i=1}^{N}\bm{y}_{i}\bm{y}_{i}^{H}.$ It follows a $L$-variate complex Wishart distribution [27, 31] with $N$ degrees of freedom and mean $\Sigma$:

$$N\hat{\Sigma}\stackrel{{\scriptstyle d}}{{\sim}}\mathbb{C}\mathcal{W}_{L}(N,\Sigma).$$

(4)

The density of a complex Wishart distribution can be found in [31]. In the next section, we use it to form the likelihood function of the data $\hat{\Sigma}$ and derive maximum likelihood estimates of the covariance kernel $\kappa$ and the noise level $\sigma$.

The log-likelihood function for $\kappa$ and $\sigma$ given the sufficient statistic $\hat{\Sigma}$ [32] can be written in terms of the density function of the complex Wishart distribution [31],

$$\ell\left(\kappa,\sigma|\hat{\Sigma}\right)=-\mbox{Tr}\left[\left(\Phi^{\ast}\mathcal{T}_{\kappa}\Phi+\sigma I_{L}\right)^{-1}\hat{\Sigma}\right]-\log\left|\Phi^{\ast}\mathcal{T}_{\kappa}\Phi+\sigma I_{L}\right|,$$

(5)

where the terms independent of $\kappa$ and $\sigma$ have been dropped. As $\sigma>0$, it is guaranteed that the matrix $\Phi^{\ast}\mathcal{T}_{\kappa}\Phi+\sigma I_{L}$ is invertible and that the log-likelihood function is hence well-defined. Maximum likelihood estimates for $\kappa$ and $\sigma$ are then obtained by maximising Eq. 5 with respect to $\kappa\in\mathcal{L}^{2}(\mathbb{S}^{2}\times\mathbb{S}^{2})$ and $\sigma>0$. Since the sampling operator $\Phi^{\ast}$ has finite rank and consequently a non-trivial kernel, the log-likelihood function admits infinitely many local maxima. Indeed, for $f\in\mathcal{N}(\Phi^{\ast})$, adding a kernel of the form¹¹1The tensor product $\otimes$ is defined as $\left(\bar{f}\otimes f\right)g:=\langle g,f\rangle f,\quad\forall f,g\in\mathcal{L}^{2}(\mathbb{S}^{2},\mathbb{C})$. $\bar{f}\otimes f$ to $\mathcal{T}_{\kappa}$ in (5) does not change the value of the log-likelihood function. We thus choose to impose a unique maximum by restricting the search space for $\kappa$ to a lower dimensional subspace, and look for solutions in the range of some synthesis operator $\bar{\Psi}\otimes\Psi$, which will be specified in Section 3.4:

$$\kappa=\left(\bar{\Psi}\otimes\Psi\right)\mbox{vec}(R)=\sum_{i,j=1}^{M}R_{ij}\;\bar{\psi}_{j}\otimes\psi_{i},\;\Leftrightarrow\;\mathcal{T}_{\kappa}=\Psi R\Psi^{\ast},$$

where $R\in\mathbb{C}^{M\times M}$ is a Hermitian symmetric matrix and $\Psi^{\ast}:\mathcal{H}\rightarrow\mathbb{C}^{M}$, $\Psi:\mathbb{C}^{M}\rightarrow\mathcal{H}$ are the analysis and synthesis operators associated with the family of functions $\{\psi_{1},\ldots,\psi_{M}\}\subset\mathcal{H}$. This regularisation of the likelihood problem by restricting the parameter space to a lower dimensional subspace is generally known as the method of sieves [22, 23]. The maximum likelihood estimates of $R$ and $\sigma$ are then given by minimising the negative log-likelihood:

$$\hat{R},\hat{\sigma}=\mbox{arg}\min_{\begin{subarray}{c}R\in\mathbb{C}^{M^{2}}\\ \sigma>0\end{subarray}}\mbox{Tr}\left[\left(GRG^{H}+\sigma I_{L}\right)^{-1}\hat{\Sigma}\right]+\log\left|GRG^{H}+\sigma I_{L}\right|,$$

(6)

where $G=\Phi^{\ast}\Psi\in\mathbb{C}^{L\times M}$ is the so-called Gram matrix [19], given by $(G)_{ij}=\langle\psi_{j},\phi_{i}\rangle$. For Eq. 6 to admit a unique solution, it is necessary to have at least as many measurements as unknowns. When the noise power is unknown a priori, this requires that $M<L$. When the noise power is known, there is one less unknown, leading to $M\leq L$. This is however not a sufficient condition for identifiability, and we must further assume $G$ to be of full column-rank. If the latter condition is verified, we say that the two families of functions $\{\phi_{1},\ldots,\phi_{L}\}$ and $\{\psi_{1},\ldots,\psi_{M}\}$ are coherent with one another.

Suppose the noise power $\sigma$ is known. Then $R$ becomes the only variable in Eq. 6, and a solution can easily be obtained by cancelling the derivative. This yields

$$\hat{R}=G^{\dagger}\tilde{\Sigma}\left(G^{\dagger}\right)^{H}=G^{\dagger}\left[\hat{\Sigma}-\sigma I_{L}\right]\left(G^{\dagger}\right)^{H},$$

where $G^{\dagger}\in\mathbb{C}^{L\times M}$ is the left pseudo-inverse²²2The left pseudo-inverse exists since $G$ is assumed full-column rank. [33] of $G$. Hence, when restricting the search space to $\mathcal{R}(\bar{\Psi}\otimes\Psi)$, the maximum likelihood of $\kappa$ is given by

	$\displaystyle\hat{\kappa}$	$\displaystyle=\sum_{i,j=1}^{M}\hat{R}_{ij}\;\bar{\psi}_{j}\otimes\psi_{i}$
		$\displaystyle=\left(\bar{\Psi}\otimes\Psi\right)\mbox{vec}\left(G^{\dagger}\left[\hat{\Sigma}-\sigma I_{L}\right]\left(G^{\dagger}\right)^{H}\right)$
		$\displaystyle=\left(\bar{\Psi}\otimes\Psi\right)\left[\bar{G}^{\dagger}\otimes G^{\dagger}\right]\left(\hat{\bm{\varsigma}}-\sigma\bm{\epsilon}\right),$		(7)

with $\hat{\bm{\varsigma}}=\mbox{vec}(\hat{\Sigma})\in\mathbb{C}^{L^{2}}$ and $\bm{\epsilon}=\mbox{vec}(I_{L})\in\mathbb{C}^{L^{2}}.$ The intensity function is then obtained by taking the diagonal part of $\hat{\kappa}:$

$$\hat{I}(\bm{r})=\sum_{i,j=1}^{M}\hat{R}_{ij}\psi_{i}(\bm{r})\bar{\psi}_{j}(\bm{r}),\quad\forall\bm{r}\in\mathbb{S}^{2}.$$

Using properties of the tensor product and the vec operator, we can re-write Eq. 3 as $\bm{\varsigma}=\left(\bar{\Phi}\otimes\Phi\right)^{\ast}\kappa\;+\;\sigma\bm{\epsilon}.$ Hence, since $\mathbb{E}[\hat{\bm{\varsigma}}]=~\bm{\varsigma}$, Eq. 7 becomes on expectation

$$\mathbb{E}[\hat{\kappa}]=\left(\bar{\Psi}\otimes\Psi\right)\left[\bar{G}^{\dagger}\otimes G^{\dagger}\right]\left(\bar{\Phi}\otimes\Phi\right)^{\ast}\kappa.$$

For $M=L$, $G$ is invertible and $G^{\dagger}=G^{-1}$, making $(\bar{\Psi}\otimes\Psi)[\bar{G}^{-1}\otimes G^{-1}](\bar{\Phi}\otimes\Phi)^{\ast}$ an oblique projection operator [19]. The operator $(\bar{\Psi}\otimes\Psi)[\bar{G}^{-1}\otimes G^{-1}]$ is indeed a right-inverse of $(\bar{\Phi}\otimes\Phi)^{\ast}$:

$$\left(\bar{\Phi}\otimes\Phi\right)^{\ast}\left(\bar{\Psi}\otimes\Psi\right)\left[\bar{G}^{-1}\otimes G^{-1}\right]=\overline{\Phi^{\ast}\Psi G^{-1}}\otimes\Phi^{\ast}\Psi G^{-1}=I_{L^{2}}.$$

(8)

In the specific case where $M=L$, the maximum likelihood estimate $\hat{\kappa}$ is hence an unbiased, consistent and asymptotically efficient estimator of the oblique projection of $\kappa$ onto $\mathcal{R}(\bar{\Psi}\otimes\Psi)$. When additionally setting $\Psi=\Phi$ the projection becomes orthogonal.

Suppose now the noise power is unknown. We must hence minimise Eq. 6 with respect to both $\sigma$ and $R$. Using the result from theorem 1.1 of [16], we can write explicit solutions for the unique minimisers of Eq. 6:

$\displaystyle\hat{\sigma}=\frac{\mbox{Tr}\left(\hat{\Sigma}-GG^{\dagger}\hat{\Sigma}\right)}{L-M},\qquad\hat{R}=G^{\dagger}\left[\hat{\Sigma}-\hat{\sigma}I\right]{\left(G^{\dagger}\right)}^{H}.$

(9)

Again, the constrained maximum likelihood estimate of $\kappa$ is given by

$$\hat{\kappa}=\left(\bar{\Psi}\otimes\Psi\right)\left[\bar{G}^{\dagger}\otimes G^{\dagger}\right]\left(\hat{\bm{\varsigma}}-\hat{\sigma}\bm{\epsilon}\right),$$

(10)

with intensity function $\hat{I}(\bm{r})=\sum_{i,j=1}^{M}\hat{R}_{ij}\psi_{i}(\bm{r})\bar{\psi}_{j}(\bm{r}).$ This time, since $M<L$ the consistency condition Eq. 8 cannot be met, and $\mathbb{E}[\hat{\kappa}]$ can no longer be interpreted as an oblique projection of $\kappa$. For values of $M$ comparable to $L$ though, the consistency condition should still hold approximately³³3More precisely, the consistency condition will hold on a subspace of $\mathbb{C}^{L}$ of dimension $M$. , and this geometrical interpretation provides intuition.

We have thus far only required the synthesis operator $\Psi$ to be identifiable, with the coherency condition requiring $G=\Phi^{\ast}\Psi$ to be full column-rank. This still leaves plenty of potential candidates. For practical purposes, we recommend taking $\Psi=\Phi W$ where $W\in\mathbb{C}^{L\times M}$ is a tall matrix, with columns containing the first $M$ eigenvectors of $\hat{\Sigma}$ (assuming eigenvalues sorted in descending order). Such a choice presents numerous advantages. First, since $\mathcal{R}(\Phi)^{\perp}=\mathcal{N}(\Phi^{\ast})$, the instrument can only sense functions within the range of $\Phi$, and it is hence natural to choose $\mathcal{R}(\Psi)=\mathcal{R}(\Phi)$. This canonical choice moreover yields an analytically computable Gram matrix $G$. Indeed, we have $G=\Phi^{\ast}\Phi W=HW$, where $H\in\mathbb{C}^{L\times L}$ is given by (see of [7, Chapter 4 section 1.1]):

$$(H)_{ij}=~4\pi\;\mbox{sinc}(2\pi\|\bm{p}_{i}-\bm{p}_{j}\|_{2}/\lambda),\quad i,j=1,\ldots,L.$$

Finally, by choosing the columns of $W$ as the first $M$ eigenvectors of $\hat{\Sigma}$, $M$ acts as a regularisation parameter. Indeed, the eigenvectors associated to the smallest eigenvalues of $\hat{\Sigma}$ are usually the most polluted by noise. Hence, truncating to the $M$ largest eigenvalues reduces the contribution of the noise in the final estimate (see Figs. 2f, 2g and 2h). Moreover, small values of $M$ will increase the chances of $(G^{H}G)\in\mathbb{C}^{M\times M}$ in the left pseudo-inverse $G^{\dagger}=(G^{H}G)^{-1}G^{H}$ being well-conditioned, thus improving the overall numerical stability of the algorithm. Suitable values of $M$ can be obtained by minimising the Bayesian Information Criterion (BIC) [24], often used in model selection: $BIC(M)=-2\hat{\ell}_{M}+2M^{2}\log(L),$ where $\hat{\ell}_{M}$ is the maximised log-likelihood function for a specific choice of $M$. Example of a BIC profile and evolution of the BIC-selected $M$ with the signal-to-noise ratio are depicted in Fig. 1.

Fig. 2 compares the performance of the proposed Sieved Maximum Likelihood (SiML) method in a radio astronomy setup to three popular spectral-based methods, namely Matched Beamforming (MB), Maximum Variance Distortionless Response (MVDR) and the Adapted Angular Response (AAR) [34]. For this experiment, we generated randomly a layout $L=300$ antennas and simulated $N=2000$ random measurements from the ground truth intensity field Fig. 2a. Furthermore, we considered two metrics to assess the quality of the recovered images: the traditional relative Mean Squared Error (MSE) and the Root Mean Squared (RMS) metric, which measures the contrast of an image by computing its standard deviation over all pixels. The simulations reveal that the SiML outperforms all the traditional algorithms for the considered SNR range in both metrics, except for large SNRs where MVDR exhibits a slightly better contrast. As for the traditional SML method, SiML performs particularly well for challenging scenarios with very low SNR.