wGMCA

The wGMCA algorithm aims at solving joint deconvolution and blind source separation (DBSS) problems from non-coplanar interferometeric data.

Contents

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

Introduction

The wGMCA builds upon the w-stacking framework [Offringa et all, 2014]. The w-axis is discretized uniformly into values, and for each channel, the interferometric samples are assigned to their nearest w-plane. Next, for each channel and for each w-plane, the interferometric samples are gridded along the (u,v) axes on a uniform grid of size and then flattened in a vector of size . This leads to the obtaining of a three-dimensional tensor , with the number of channels.

The mixing model is described by (see Reference for a more precise formulation):

,

where:

• is the interferometer response (in the visibility space),
• denotes the element-wise product,
• is a tensor operator based on the two-dimensional (fast) Fourier transform,
• is a tensor operator which accounts for the non-coplanar effect,
• is the mixing matrix,
• are the sources, flattened and stacked in a matrix,
• is a complex Gaussian noise with known variances.

The sources are assumed to be sparse in the starlet domain . Moreover, both the sources and the mixing matrix are supposed to be constituted of nonnegative elements.

The wGMCA algorithm aims at minimizing the following objective function with respect to and :

where is a quadratic form which depends on the noise variance, contains the sparsity regularization parameters, and are the nonnegative orthants for sources and mixing matrices, and ensures that the columns of the mixing matrix have a norm less than or equal to unity.

Procedure

The wGMCA algorithm is based on GMCA, which is a BSS procedure built upon a projected alternating least-square (pALS) minimization scheme. In brief, the updates of and comprise a least-square estimate, to minimize the data-fidelity term, 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 (i) generally intractable and (ii) not necessarily stable with respect to noise. Thus, (i) the least-square update is decomposed into two simpler updates and (ii) 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 provides robustness with respect to the initial point. The second stage refines the separation by employing more precise strategies (refinement). Lastly, the sources S are improved during a finale step with the output mixing matrix.

Getting started

Requirements

• Python (last tested with v3.7.10)
• NumPy (last tested with v1.20.2)
• SciPy (last tested with v1.7.0)

wGMCA class

The wGMCA algorithm is implemented in a class wGMCA. The data and the algorithm parameters are provided at the object initialization. The DBSS is performed by running the method run. The results are stored in the attributes A and S.

Parameters

Below are the five parameters of the wGMCA class that must always be provided at initialization.

Parameter	Type	Information	Default value
X	(m,w,p) complex numpy.ndarray	input data (gridded visibilities); 1st axis: channel, 2nd axis: w axis, 3rd axis: flattened (u,v) axes	N/A
H	(m,w,p) float numpy.ndarray	interferometer response in visibility domain for several w-values, with zero-frequency shifted to the center and flattened	N/A
n	int	number of sources to be estimated	N/A
Var	(m,w,p) float array or float	noise variance (wGMCA does not account for potential noise covariances)	N/A
G	(w,p) float array	w-term matrices	N/A

Below are the essential parameters of the wGMCA class. They may be assigned their default value.

Parameter	Type	Information	Default value
nnegA	bool	non-negativity constraint on	True
nnegS	bool	non-negativity constraint on	True
nneg	bool	non-negativity constraint on and , overrides nnegA and nnegS if not None	None
c_wu	float	Tikhonov regularization hyperparameter at warm-up	0.5
c_ref	float	Tikhonov regularization hyperparameter at refinement	0.5
c_end	float	Tikhonov regularization hyperparameter at finale refinement of S	0.5
itCG	int	maximum number of iterations of the conjugate gradient algorithm during refinement (useful for strong non-coplanar effects)	100
nscales	int	number of starlet detail scales	2
k	float	parameter of the k-std thresholding (~1 for approximately sparse sources, 3 for very sparse sources)	3
K_max	float	maximal L0 norm of the sources. Being a percentage, it should be between 0 and 1 (small for very sparse sources, 1 for approximately sparse sources)	0.5
K_end	float	maximal L0 norm of the sources during finale refinement of S. Being a percentage, it should be between 0 and 1	1
thr_end	bool	perform thresholding during the finale refinement of S (consider False if significant information lies in the coarse scales)	True
doRw	bool	do L1 reweighing during refinement (consider False if algorithm unstable)	True
doRw_end	bool	do L1 reweighing during finale refinement of S	doRw
H_reconv	(,p) float numpy.ndarray	kernel in Fourier domain, with zero-frequency shifted to the center and flattened, by which the sources are reconvolved for the estimation of (typically to get the resolution of the first channel)	None
removeCoarseScaleData	bool	remove coarse scale from data for the estimation of (implemented very approximately, consider True only if resolution ratio is ≳ 0.5)	False
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 are other parameters of the wGMCA class, which can reasonably be assigned their default value.

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	False
cstWuRegStr	bool	use constant regularization coefficients during warm-up	False
minWuIt	int	minimum number of iterations at warm-up	50
useMad	bool	estimate noise std in source domain with MAD (else: analytical estimation, use with caution)	False
useMad_end	bool	estimate noise std in source domain with MAD during finale refinement of S (else: analytical estimation, use with caution)	False
L1	int	L1 penalization (else: L0 penalization)	True
S0	(n,p) float array	ground truth sources (for testing purposes)	None
A0	(m,n) float array	ground truth mixing matrix (for testing purposes)	None

Example

Perform a DBSS on the gridded data X with four sources. The interferometer response is stored in H and the data variance inVar.

    wgmca = wGMCA(X, H, 4, Var, K_max=0.5, nscales=3)
    wgmca.run()
    S = wgmca.S.copy()
    A = wgmca.A.copy()

Authors

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

Reference

TODO: add ref

License

This project is licensed under the LGPL-3.0 License.