SDecGMCA

SDecGMCA (Spherical Deconvolution Generalized Morphological Component Analysis) is an algorithm aiming at solving joint Deconvolution and Blind Source Separation (DBSS) problems with spherical data.

The 2D adaptation of the algorithm can be found here.

Contents

1. Introduction
2. Procedure
3. Getting Started
4. Parameters
5. Example
6. Authors
7. Reference
8. License

Introduction

Let us consider the imaging model:

,

where:

• are the multiwavelength spherical data, stacked in a matrix,
• is measurement operator,
• is the mixing matrix,
• are the sources, stacked in a matrix,
• is the noise.

In the following, the measurement operator is assumed channel-dependant, linear and isotropic. Thus, for channel :

where denotes the convolution product on the sphere. The equation can be simplified in the spherical harmonics domain for each harmonic coefficient and for all channels:

.

The sources are assumed to be sparse in the starlet representation . SDecGMCA aims at minimizing the following objective function with respect to A and S:

,

where denotes the element-wise product, are the sparsity regularization parameters and is the oblique ensemble. Moreover, non-negativity constraints on and/or may be added.

Procedure

The algorithm is built upon a sparsity-enforcing projected alternate least-square procedure, which updates iteratively the sources and the mixing matrix. In a nutshell, when either or is updated, a first least-squares estimate is computed by minimizing the data-fidelity term. This step is then followed by the application of the proximal operator of the corresponding regularization term.

In contrast to standard BSS problems, the least-square update of the sources is not necessarily stable with respect to noise. Thus, an extra Tikhonov regularization is added.

The separation is comprised of two stages. The first stage estimates a first guess of the mixing matrix and the sources (warm-up); it is required to provide robustness with respect to the initial point. The second stage refines the separation by employing a more precise Tikhonov regularization strategy (refinement). Lastly, the sources are improved during a finale step with the output mixing matrix.

During the warm-up, the Tikhonov regularization coefficients are calculated from the current estimation of the mixing matrix. During the refinement, they are rather calculated from the current estimation of the sources spectra.

Getting Started

Requirements

SDecGMCA has been developed with Python 3.7. It depends on the Python Healpy library, which is only supported on Linux and MAC OS X, not Windows.

Prerequisites

The following Python libraries need to be installed to run the code:

• NumPy,
• Healpy (see https://healpy.readthedocs.io/en/latest/install.html for the installation procedure),
• Matplotlib [optional].

SDecGMCA class

SDecGMCA is implemented in a class. The data and the parameters of the separation are provided at the initialization of the object. The separation is performed by running the method run. The results are stored in the attributes A and S (the harmonic projection of the sources is in addition stored in Slm).

Parameters

Below is the list of the main attributes of the SDecGMCA class.

Parameter	Type	Information	Default value
X	(m,p) float or (m,t) complex numpy.ndarray	input data in Healpix representation or in spherical harmonic domain, each row corresponds to an observation	N/A
Hl	(m,lmax+1) float numpy.ndarray	convolution kernels in spherical harmonics domain (lmax : maximum frequency)	N/A
n	int	number of sources to be estimated	N/A
alm_in	bool	the data X are in the spherical harmonic domain (Healpix nside=lmax/3)	False
M	(p,) float numpy.ndarray	mask in Healpix representation	None
nnegA	bool	non-negativity constraint on A	False
nnegS	bool	non-negativity constraint on S	False
nneg	bool	non-negativity constraint on A and S. If not None, overrides nnegA and nnegS	None
c_wu	float or (2,) numpy.ndarray	Tikhonov regularization hyperparameter at warm-up	0.5
c_ref	float	Tikhonov regularization hyperparameter at refinement	0.5
nStd	float	noise standard deviation (compulsory for spectrum-based regularization & analytical estimation of the noise std in source domain)	N/A
useMad	bool	True: noise std estimated with MAD, False: noise std estimated analytically	False
nscales	int	number of starlet detail scales	3
k	float	parameter of the k-std thresholding	3
K_max	float	maximal L0 norm of the sources. Being a percentage, it should be between 0 and 1	0.5
thrEnd	bool	perform thresholding during the finale estimation of the sources	True
eps	(3,) float numpy.ndarray	stopping criteria of (1) the warm-up, (2) the refinement and (3) the finale refinement of S	[1e-2, 1e-4, 1e-4]
verb	int	verbosity level, from 0 (mute) to 5 (most talkative)	0

Below is the list of the other attributes of the SDecGMCA class, which can reasonably be equal to their default values.

Parameter	Type	Information	Default value
AInit	(m,n) float numpy.ndarray	initial value for the mixing matrix. If None, PCA-based initialization	None
keepWuRegStr	bool	keep warm-up regularization strategy during refinement (else: spectra-based coefficients)	False
cstWuRegStr	bool	use constant regularization coefficients during warm-up (else: mixing-matrix-based coefficients)	False
minWuIt	int	minimum number of iterations at warm-up	100
cwuDec	int	number of iterations for the decrease of c_wu (if c_wu is an array)	minWuIt/2
L1	bool	if False, L0 rather than L1 penalization	True
doRw	bool	do L1 reweighing during refinement (only if L1 penalization)	True
iterSH	int	number of iterations for the spherical harmonic transforms	3

Example

Perform a DBSS on the data X with 4 sources.

sdecgmca = SDecGMCA(X, Hl, n=4, k=3, nscales=3, c_wu=numpy.array([5, 0.5]), c_ref=0.5, nStd=1e-7)
sdecgmca.run()
S = sdecgmca.S.copy()  # estimated sources
A = sdecgmca.A.copy()  # estimated mixing matrix

Authors

• Rémi Carloni Gertosio
• Jérôme Bobin

Reference

R. Carloni Gertosio, J. Bobin, Joint deconvolution and unsupervised source separation for data on the sphere

License

This project is licensed under the LGPL-3.0 License.