Sparse component separation from Poisson measurements

TL;DR

This paper introduces pGMCA, a novel blind source separation algorithm designed specifically for sparse signals contaminated with Poisson noise, addressing challenges in low photon count and high-energy astronomical imaging.

Contribution

The paper presents a new BSS method tailored for Poisson measurements, filling a gap in existing Gaussian noise-based approaches.

Findings

Effective separation of sparse sources from Poisson data.

Applicable to low photon count and astronomical imaging.

Demonstrates improved performance over traditional methods.

Abstract

Blind source separation (BSS) aims at recovering signals from mixtures. This problem has been extensively studied in cases where the mixtures are contaminated with additive Gaussian noise. However, it is not well suited to describe data that are corrupted with Poisson measurements such as in low photon count optics or in high-energy astronomical imaging (e.g. observations from the Chandra or Fermi telescopes). To that purpose, we propose a novel BSS algorithm coined pGMCA that specifically tackles the blind separation of sparse sources from Poisson measurements.

Equations

$$\underset{\mathbf{A}\in \mathcal{C},\mathbf{S}\ge 0}{\min }{\Vert \mathbf{\Lambda }\odot \mathbf{S}{\mathbf{\Phi }}^T\Vert }_{{\ell }_1}+\mathcal{L}\left(\mathbf{X}|\mathbf{A},\mathbf{S}\right),$$

$$\underset{\mathbf{A}\in \mathcal{C},\mathbf{S}\ge 0}{\min }{\Vert \mathbf{\Lambda }\odot \mathbf{S}{\mathbf{\Phi }}^T\Vert }_{{\ell }_1}+{\mathcal{L}}_{\mu }\left(\mathbf{X}|\mathbf{A},\mathbf{S}\right).$$

Full text

Introduction

Multichannel data are composed of $m$ observations ${\bf X}_{i}$, each of which is made of $t$ samples. According the standard instantaneous linear mixture, each observation is described as a linear combination of $n$ elementary sources or components ${\bf S}_{j}$. It is classically assumed that the data are corrupted with additive, generally Gaussian, noise, which leads to the following matrix formulation ${\bf X}={\bf AS}+{\bf N}$, where ${\bf X}\in\mathbb{R}^{m\times t}$ is the observations matrix, ${\bf S}\in\mathbb{R}^{n\times t}$ the sources matrix, ${\bf A}\in\mathbb{R}^{m\times n}$ the mixing matrix and ${\bf N}\in\mathbb{R}^{m\times t}$ for the noise contribution. In this context, BSS aims at recovering both the mixing matrix $\bf A$ and the sources $\bf S$ from the data $\bf X$ only. This is essentially an unsupervised matrix factorization; being ill-posed it requires additional assumptions about the sources and/or mixing matrix such as statistical independence [13], non-negativity [6] or sparsity [9, 2] to only name three. However, the above linear mixture model does not describe precisely data that are commonly found in low photon count imaging, such as in X-ray [1]. For that purpose, one needs to account for the exact statistics of the measurements, which precisely follow a Poisson distribution. Hence, the data $\bf X$ are only defined statistically from the ”pure” mixtures ${\bf AS}$; the probability for a given sample to take the value ${\bf X}_{i}[t]$ is given by the Poisson law: $\mathcal{P}({\bf X}_{i}[t]\lvert[{\bf AS}]_{i}[t])=\frac{e^{-[{\bf AS}]_{i}[t]}\,[{\bf AS}]_{i}[t]^{{\bf X}_{i}[t]}}{{\bf X}_{i}[t]!}$, where $[{\bf AS}]_{i}[t]$ is the sample of the matrix $\bf AS$ located at the $i$-th row and $t$-th column. To tackle BSS from Poisson measurements, a straightforward approach consists in maximizing the likelihood of mixture variable. In the case of Poisson statistics, this amounts to minimizing the Kullback-Leibler divergence between the data $\bf X$ and the mixture model ${\bf AS}$ with respect to $\bf A$ and $\bf S$, which has been investigated both in the scope of Independent Component Analysis (ICA - [7]) and Non-negative Matrix Factorization (NMF - [6]) where it generally refers to robust BSS. In the scope of sparse BSS, the sources are assumed to admit a sparse distribution in some signal representation $\bf\Phi$. While sparse BSS has been successful in various applications [4], to the best of our knowledge, this blind separation of sparse sources has not been investigated when the data follow a Poisson distribution. We therefore propose a new algorithm coined poisson Generalized Morphological Component Analysis (pGMCA) to tackle sparse BSS from Poisson measurements.

The pGMCA algorithm

In the next, we assume each source admits a sparse representation in some signal representation $\bf\Phi$, which will be assumed to be an orthogonal basis for the sake of simplicity. Next, the sparsity of the sources will be enforced by minimizing their re-weighted $\ell_{1}$ norm as follows:

[TABLE]

where the matrix $\bf\Lambda$ contains the regularization parameters as well as weigths for re-weighted $\ell_{1}$ regularization [15]. The second term is the anti-log likelihood of variables $\bf A$ and $\bf S$: $\mathcal{L}\left({\bf X}|{\bf A},{\bf S}\right)=\sum_{i,t}[{\bf AS}]_{i}[t]-{\bf X}_{i}[t]\log\left([{\bf AS}]_{i}[t]\right)$. The convex set $\mathcal{C}$ is the intersection of the positivity constraint ${\bf A}\geq 0$ and the multi-dimensional $\ell_{2}$-ball constraint that enforces each column of the mixing matrix to have a $\ell_{2}$ norm lower than $1$ so as to alleviate the scale indeterminacy of the mixture model. The norm $\|\,.\,\|_{F}$ is the Frobenius norm. Optimizing the problem in Eq. 1 raises several challenges:

**- A multi-convex problem: ** the problem is non-convex but convex with respect to each variable $\bf A$ and $\bf S$ assuming the other one is fixed. Hopefully, several optimization strategies have been proposed so far to tackle multi-convex problem such as the Block-Coordinate-Descent algorithm [16, 12], which sequentially optimizes over each variable independently.

**- Non-differentiability and curvature of the data fidelity term: ** the data fidelity term $\mathcal{L}\left({\bf X}|{\bf A},{\bf S}\right)$ is not smooth about [math], which rigorously excludes the use of BCD. Furthermore, its curvature at the vicinity of [math] increases at a hyperbolic rate. Consequently, neglecting the non-differentiability of $\mathcal{L}$ at [math] to use the BCD would yield to a dramatically slow convergence rate. This would be especially true for measurements that correspond to low flux values.

To make use of the BCD, we rather propose to resort to a smooth approximation of $\mathcal{L}$ based on the smoothing technique introduced by Nesterov in [8]. Following Nesterov’s approach, a smooth approximate $\mathcal{L}_{\mu}$ of $\mathcal{L}$ can be built from its Fenchel dual $\mathcal{L}^{\star}$ [14]:$\mathcal{L}_{\mu}({\bf Y})=\mbox{inf}_{\bf U}\,\langle{\bf Y},{\bf U}\rangle-\mathcal{L}^{\star}({\bf U})-\mu\|{\bf U}\|_{F}^{2}$. The resulting approximate $\mathcal{L}_{\mu}({\bf Y})$ is differentiable and its gradient is $\mu$-Lipschitz, which makes the application of BCD possible. The approximate sparse BSS problem to tackle is the following:

[TABLE]

Building upon the BCD algorithm, the pGMCA algorithm (e.g. poisson-GMCA) sequentially updates each variable $\bf S$ and $\bf A$; it is described in Alg.1.

Step $1$ of the pGMCA algorithm requires minimizing the sum of a smooth data fidelity term with two extra non-smooth penalization, namely the $\ell_{1}$ norm of the sources in $\bf\Phi$ and the positivity constraint. In the general case, this step does not admit a closed-form expression. Therefore, at each iteration $k$, the estimate ${\bf S}^{(k+1)}$ is obtained numerically using an implementation of the Generalised Forward-Backward Splitting (G-FBS) algorithm. Similarly, step $2$ involves the minimization of $\mathcal{L}_{\mu}$ subject to ${\bf A}\in\mathcal{C}$, which does not admit a closed-form expression. Since $\mathcal{C}$ is the intersection of a $\ell_{2}$-ball with unit norm and the positivity orthant, it can be shown that its proximal operator [5] is explicit. Therefore, step $2$ is solved using a FISTA-like algorithm [18]. Practical details of the implementation are described below:

**- Initialization: ** The problem in Eq. 2 being multi-convex, the BCD algorithm can be quite sensitive to the initial point. In the next section, the pGMCA algorithm is initialized with the output of GMCA, which generally provides a robust first guess estimate.

**- Regularization parameters: ** In practice, the BCD turns out to be quite sensitive to the choice of the regularization parameter. In these preliminary investigations, we did not implement any $\ell_{1}$ re-weighing scheme; the matrix $\bf\Lambda$ is therefore constant for all the entries of a given source and varies from one source to the other. In practice the regularisation parameter $\lambda_{i}$ for the $i$-th source is fixed based on the level of the noise that propagates in a single iteration of the G-FBS. More precisely, it is chosen as $\lambda_{i}=\tau.\mbox{MAD}\left(\nabla_{\bf S}\mathcal{L}_{\mu}({\bf A}^{(0)},{\bf S}^{(0)}){\bf\Phi}^{T}\right)$ where the Median-Absolute-Deviation (MAD) is an empirical estimate of the standard deviation of the noise measured in the gradient in the sparsifying representation. Interestingly, while this is suited for Gaussian noise, such a strategy works quite fine in Poisson case for $\tau=1$ thanks to the use of a smooth approximation of $\mathcal{L}$. The smoothing parameter mainly depends on the Poisson density $\bf AS$; a small value will yield a more precise approximation but at the cost of a slower convergence (i.e. the gradient step size in step $1$ and $2$ scales like $\mu$). In the next, it has been set to the mean number of measures counts.

**- Stopping criterion and number of iterations: ** The number of iterations is fixed $10000$. The algorithm steps when the Frobenius norm between two consecutive estimates of the mixing matrix $\left\|{\bf A}^{(k+1)}-{\bf A}^{(k)}\right\|_{F}$ is lower than $\epsilon=10^{-6}$.

Numerical experiments

In this section, numerical experiments are carried out on simulations of astrophysical data that have been generated from real Chandra observations of the Cassiopea A supernova remnants. These data are composed of a linear combination of $3$ components ($128\times 128$ pixels): synchrotron emission, and $2$ redshifted iron emission lines. They are composed of $50$ observations. The electromagnetic spectrum of the synchrotron emission follows a power law while the two iron components have Gaussian-shaped spectra, which are the instrumental response of lines at different redshifted, which are representative of these data in the energy band $5000-6000$ eV (electron-volt). The sources and their spectra are displayed in Fig. 1.

In this abstract, we investigate the performances of the various component separation methods when the mean flux (mean number of counts) evolves. These methods are GMCA, pGMCA, $\beta$-NMF [6] and $\beta$-ICA [7]. In the last two ones, the parameter $\beta$ of the $\beta$-divergence is set to $1$ to minimize a Kullback-Leibler divergence, which is well suited for Poisson statistics. $\bf\Phi$ is chosen as isotropic undecimated wavelet [19]. Since these methods do no impose similar regularization on the sources, they are more fairly compared based on the quality of estimation of the mixing matrix. For that purpose , we make use of the spectral angular distance (SAD) between the estimated $[\hat{\bf A}]^{j}$ and input $[{\bf A}^{o}]^{j}$ columns of the mixing matrix: $\mbox{SAD}=\sum_{j=1}^{n}\cos^{-1}(|\langle[\hat{\bf A}]^{j},[{\bf A}^{o}]^{j}\rangle|)/n$. The results are displayed in Fig.2; each point is the mean valued of $10$ Monte-Carlo simulations with different noise realisations. This figure first shows that the robust ICA and robust NMF methods perform badly; the range of mean flux values (i.e. $0.5$ to $35$) correspond to a rather severe noise level for component separation, which could explain these poor results. For mean flux values larger than $2$ (i.e. on average a number counts equal to $2$ per pixel), GMCA performs satisfactorily and pGMCA yields a gain of about $1$ order of magnitude with respect to GMCA. For lower number of counts, it is likely that the noise level is too large to perform a reasonable separation process.

Conclusion

We investigated the development of a sparsity enforcing method to tackle BSS from measurements with Poisson noise. The proposed pGMCA algorithm builds upon a BCD-like algorithm that makes use a smooth approximation of the Poisson log-likelihood. Preliminary results show a clear improvement of the separation quality with respect to existing methods when the number of counts is low. During the workshop, we will show more extensive results that better highlight the benefits of the proposed pGMCA algorithm.