Heuristics for Efficient Sparse Blind Source Separation

TL;DR

This paper introduces a heuristic-enhanced optimization strategy for sparse blind source separation, improving practical results especially in astrophysical data analysis, by addressing limitations of existing algorithms like PALM.

Contribution

It proposes a new heuristic approach combined with PALM to enhance the practical performance of sparse BSS algorithms, particularly for complex astrophysical data.

Findings

Improved separation results on astrophysical data

Heuristic approach enhances PALM performance

Demonstrated relevance in realistic data scenarios

Abstract

Sparse Blind Source Separation (sparse BSS) is a key method to analyze multichannel data in fields ranging from medical imaging to astrophysics. However, since it relies on seeking the solution of a non-convex penalized matrix factorization problem, its performances largely depend on the optimization strategy. In this context, Proximal Alternating Linearized Minimization (PALM) has become a standard algorithm which, despite its theoretical grounding, generally provides poor practical separation results. In this work, we propose a novel strategy that combines a heuristic approach with PALM. We show its relevance on realistic astrophysical data.

Tables

Table 1. Table 1: Average C A subscript 𝐶 𝐴 C_{A} (10 different initializations) for 3 SNR values and 5 algorithms.

	10 dB	15dB	20 dB
2 step	11.57	16.92	22.09
PALM	9.38	10.94	11.01
GMCA	8.87	12.15	15.87
EFICA	-6.92	5.11	9.41
RNA	-6.68	-5.49	6.27

Equations

$$\underset{\mathbf{A},\mathbf{S}}{\min }\frac{1}{2}{\Vert \mathbf{X}-\mathbf{AS}\Vert }_F^2+{\iota }_{\{\forall i;\Vert {\mathbf{A}}^i{\Vert }_2^2\le 1\}}(\mathbf{A})+{\Vert {\mathbf{R}}_{\mathbf{S}}\odot (\mathbf{S}{\mathbf{\Phi }}_{\mathbf{S}}^{\mathbf{T}})\Vert }_1$$

$${\left({\lambda }_1,{\lambda }_2,...,{\lambda }_n\right)}^T=k\times \mathtt{MAD}(\hat{\mathbf{S}})$$

$${\left({\lambda }_1,{\lambda }_2,...,{\lambda }_n\right)}^T\simeq k\times \frac{\gamma }{{\Vert {\mathbf{A}}^{\ast T}{\mathbf{A}}^{\ast }\Vert }_2}\mathtt{MAD}({\mathbf{A}}^{\ast T}\mathbf{N})$$

Full text

1 Introduction

1.1 Blind Source Separation problem

In the BSS [1] framework, the data are composed of $m$ observations, each of which has $t$ samples. These observations are supposed to be some linear combinations of $n$ sources. The objective of BSS is to retrieve the sources as well as the mixing coefficients. In matrix form, the goal is therefore to find two matrices $\mathbf{S}$ (of size $n\times t$) and $\mathbf{A}$ (of size $m\times n$), called respectively the source and the mixing matrices, such that: $\mathbf{X=AS+N}$, where $\mathbf{X}$ (of size $m\times t$) is the observation matrix that is corrupted with some unkwown noise $\mathbf{N}$. Since it requires tackling an ill-posed unsupervised matrix factorization problem, further assumptions are needed, including the statistical independance of the sources (ICA - [1]), the non-negativity of $\mathbf{A}$ and $\mathbf{S}$ [2]. In this work, we will focus on the sparsity of the sources [3, 4, 5, 6]. In this framework, sparse BSS will aim at finding a (local) minimum of:

[TABLE]

The first data fidelity term promotes a faithful reconstruction of the data. The use of the Froebenius norm $\left\|{.}\right\|_{F}$ stems from the assumption of a Gaussian noise $\mathbf{N}$. The second term ensures that the columns of $\mathbf{A}$ are all in the $\ell_{2}$ ball. This avoids degenerated $\mathbf{A}$ and $\mathbf{S}$ ($\left\|{\mathbf{A}}\right\|\rightarrow\infty$ and $\left\|{\mathbf{S}}\right\|\rightarrow 0$). The last term involving the Hadamard product $\odot$ enforces a $\ell_{1}$ sparsity constraint in a transformed domain $\mathbf{\Phi_{S}}$ (here, $\mathbf{\Phi_{S}}$ will be the starlet transform [7]). $\mathbf{R_{\mathbf{S}}}$ controls the trade-off between the data fidelity and the sparsity terms. It can be decomposed into $\mathbf{R_{\mathbf{S}}}=\mathbf{\mathbf{\Lambda_{S}}}\mathbf{W}$ where $\mathbf{\mathbf{\Lambda_{S}}}$ is a diagonal matrix of the regularization parameters $\lambda_{1},\lambda_{2},...,\lambda_{n}$ and $\mathbf{W}$ is a matrix used to introduce individual penalization coefficients in the context of reweighted $\ell_{1}$ [8].

1.2 Sparse BSS in practice

Since sparse BSS requires solving a penalized matrix factorization problem, the separation quality strongly depends on the optimization strategy. Different algorithmic frameworks have been used so far: projected Alternate Least-Square (ALS - [2]), PALM [10] and Block-Coordinate Descent (BCD - [9], which is not studied here due to a high computational burden):

• •

PALM algorithm: PALM is an iterative algorithm, which alternates at each iteration between a proximal gradient step on $\mathbf{A}$ and $\mathbf{S}$. PALM is proved to converge to a local minimum of (1) under mild conditions [10].

• •

GMCA algorithm: the Generalized Morphological Component Analysis (GMCA - [3]) algorithm is based on projected ALS. At each iteration (k), the gradient step appearing in PALM is replaced by a multiplication by a pseudo-inverse. Compared to PALM, GMCA cannot be proved to converge. Furthermore, even when stabilizing the output of GMCA is not really minimizing (1).

In practice, the solution of sparse BSS methods is highly sensitive to the initial point and the values of the regularization parameters, which are generally tricky to tune without any first guess of the solution. Thanks to heuristics, GMCA benefits from an automatic thresholding strategy (see Sec.2.1) which has been empirically shown to improve its robustness with respect to the initialization and give good estimations of $\mathbf{A^{*}}$ and $\mathbf{S^{*}}$. Such heuristics do not exist in PALM.

1.3 Contributions

While PALM is theoretically well rooted and yields rather fast minimization schemes (in contrast to BCD), it generally provides poor separation results. We show how PALM-based implementations can benefit from the information provided by heuristic approaches which are in practice more robust. The robustness of the proposed combined strategy is demonstrated on realistic astrophysical data.

2 Complexity of introducing heuristics in PALM

Building on the automatic thresholding strategy of GMCA, the goal of this part is to try to derive one for PALM.

2.1 Automatic parameter choice in GMCA

In GMCA, the threshold choice is performed computing the Median Absolute Deviation (MAD) for each currently estimated source $\mathbf{\hat{S}}$. The corresponding threshold (which changes during the iterations) is set to a multiple $k$ of this value:

[TABLE]

In this equation and in the remaining of this section, ${\bf\Phi_{S}^{T}}$ was supposed to be the identity matrix without loss of generality. $k$ is a positive number and the $\mathtt{MAD}$ operator is computed row-by-row with $\mathtt{MAD}(z)=\mathtt{median}(|z-\mathtt{median}(z)|)$ for $z\in\mathbb{R}^{t}$. Using this strategy, the threshold choice can be interpreted as a dense noise removal.

2.2 Introducing heuristics in PALM

Let us assume that the thresholds in PALM are computed the same way as in GMCA through Eq. (2) and that during the iterations the algorithm finds estimates that are close to both the true matrices $\mathbf{A}^{*}$ and $\mathbf{S}^{*}$. Then, due to the assumption of $\mathbf{S}^{*}$ sparsity:

[TABLE]

Therefore, using the MAD enables a thresholding of a projection of the noise $\mathbf{N}$, which yields a similar interpretation as in GMCA. However, this interpretation requires that the algorithm finds good estimates of $\mathbf{A^{*}}$ and $\mathbf{S^{*}}$ during the iterations (contrary to GMCA, which does not require to find a good estimate of $\mathbf{S^{*}}$ for the interpretation to be verified). In practice, using directly this strategy within a PALM algorithm therefore does not yield good estimates because when initialized with a random initialization the algorithm never becomes close to such good estimates.

3 Combining GMCA and PALM: a hybrid strategy

In this part, we propose to combine the best of GMCA and PALM in a two step approach. The algorithm comprehends a warm-up stage, in which GMCA is performed, followed by a refinement stage during which PALM is performed retaining as much information as possible coming from the warm-up stage.

3.1 Motivation and full description of the algorithm

Our approach is motivated by several remarks: i) PALM theoretical background: in particular, the 2-step algorithm is thus proved to converge and, once the weights $\mathbf{W}$ estimated by the warm-up stage, it truly looks for a critical point of Eq. (1); ii) GMCA robustness with regards to initialization, which will lead to a robust 2-step algorithm; iii) Benefit from GMCA solution:

• •

Since both $\mathbf{A_{GMCA}}$ and $\mathbf{S_{GMCA}}$ are close to $\mathbf{A}^{*}$ and $\mathbf{S}^{*}$, they can be used to derive the thresholds using the MAD according to the previous section.

• •

$\mathbf{S_{GMCA}}$ should already give a good approximation of the most prominent peaks. This can be exploited in the refinement stage through the introduction of reweighted L1 [8] using the reweighting matrix $\mathbf{W}$ in problem (1): $\mathbf{W}_{i}^{j}=\frac{\epsilon}{\epsilon+\frac{|\mathbf{\hat{S}}^{j}_{i}|}{\left\|{\mathbf{\hat{S}}_{i}}\right\|_{\infty}}}$, with $\epsilon$ a small constant, $\mathbf{W}_{i}^{j}$ the coefficient of $\mathbf{W}$ corresponding to the ith line and jth column and $\mathbf{\hat{S}}_{i}$ the ith line of $\mathbf{\hat{S}}$.

These remarks lead to the following 2-step algorithm:

Input : $\mathbf{X}$ (data matrix)

•

Random initialization $\mathbf{A_{0}}$ and $\mathbf{S_{0}}$

•

Warm-up stage:

$\mathbf{A_{GMCA}}$, $\mathbf{S_{GMCA}}$ = GMCA($\mathbf{X}$,$\mathbf{A_{0}}$,$\mathbf{S_{0}}$)

•

Refinement stage:

$\mathbf{A_{PALM}}$, $\mathbf{S_{PALM}}$ = PALM($\mathbf{X}$,$\mathbf{A_{GMCA}}$,$\mathbf{S_{GMCA}}$)

The initialization, thresholding strategy and reweighting information come from the warm-up stage.

4 Experiment on realistic data

4.1 Data description and experimental protocol

The goal of this part is to apply our algorithm on realistic data to show its efficiency. The $n=2$ sources come from simulations obtained from real data of Cassiopeia A supernova remnant. They each consists in a 2D image of resolution $t=$ 128 $\times$ 128 pixels, supposed to be approximately sparse in the starlet domain. The mixing is performed through a $\mathbf{A}$ matrix drawn randomly following a standard normal distribution and modified to have unit columns. Its condition number is $C_{d}=10$. There is $m=n$ observations. To increase the realism of the data and further test the algorithm, we tried three relatively low SNR values: 10, 15 and 20 dB. $k$ is set to 3, which corresponds to a classical hypothesis in terms of Gaussian noise removal.

4.2 Empirical results

The mean of of the mixing matrix criterion $C_{A}$ [11] is displayed in Table 1. The 2-step approach always achieve better results than the two classical BSS algorithms with which we performed the comparison, namely Relative Newton Algorithm (RNA) and Eficient FastICA (EFICA). It also outerpeforms both GMCA and a PALM using directly the MAD heuristic, being always better by at least 2 dB than the best of them.

In addition to the results displayed in Fig. 1, the standard deviation of $C_{A}$ over different initializations is almost 0, which shows the robustness of the algorithm.

Conclusion

In this work, we introduce a 2-step strategy combining PALM with robust heuristic methods such as GMCA. Beyond improving the robustness of PALM-based implementations with respect to initialization, the regularization parameters can be automatically set in the proposed approach. Numerical experiments on realistic data demonstrate a high separation quality and good robustness on mixings with low SNR.

Acknowledgment

This work is supported by the European Community through the grant LENA (ERC StG - contract no. 678282).