Non-convex optimization for self-calibration of direction-dependent effects in radio interferometric imaging

TL;DR

This paper introduces a novel non-convex optimization method for joint calibration and imaging in radio interferometry, effectively handling direction-dependent effects with convergence guarantees, improving imaging quality for future telescopes.

Contribution

It presents the first joint calibration and imaging algorithm for radio interferometry with proven convergence, addressing non-convexity due to direction-dependent effects.

Findings

Efficient reconstruction of point source images.

Effective imaging of complex extended sources.

Method outperforms traditional approaches in simulations.

Abstract

Radio interferometric imaging aims to estimate an unknown sky intensity image from degraded observations, acquired through an antenna array. In the theoretical case of a perfectly calibrated array, it has been shown that solving the corresponding imaging problem by iterative algorithms based on convex optimization and compressive sensing theory can be competitive with classical algorithms such as CLEAN. However, in practice, antenna-based gains are unknown and have to be calibrated. Future radio telescopes, such as the SKA, aim at improving imaging resolution and sensitivity by orders of magnitude. At this precision level, the direction-dependency of the gains must be accounted for, and radio interferometric imaging can be understood as a blind deconvolution problem. In this context, the underlying minimization problem is non-convex, and adapted techniques have to be designed. In this work, leveraging recent developments in non-convex optimization, we propose the first joint calibration and imaging method in radio interferometry, with proven convergence guarantees. Our approach, based on a block-coordinate forward-backward algorithm, jointly accounts for visibilities and suitable priors on both the image and the direction-dependent effects (DDEs). As demonstrated in recent works, sparsity remains the prior of choice for the image, while DDEs are modelled as smooth functions of the sky, i.e. spatially band-limited. Finally, we show through simulations the efficiency of our method, for the reconstruction of both images of point sources and complex extended sources. MATLAB code is available on GitHub.

Tables

Table 1. Table 1: Parameters considered in Section  7.2 for Algorithm  1 and its initialization.

Parameters		Initialization	Algorithm 1
$J_{\text{cyc}}$		101	11
$J_{{𝗨}_1}^{\left(i\right)}=J_{{𝗨}_2}^{\left(i\right)}$	if $\left(i+1\right)=0\left[\text{mod }{J}_{\text{cyc}}\right]$	0	0
	otherwise	2	5
$J_{ϵ}^{\left(i\right)}$	if $\left(i+1\right)=0\left[\text{mod }{J}_{\text{cyc}}\right]$	1000	1000
	otherwise	0	0
${\xi }_{{𝗨}_{\text{tot}}}$		${10}^{-6}$
${\xi }_{ϵ}$		${10}^{-5}$	${10}^{-5}$
$\zeta$		$2\times {10}^{-2}$	$2\times {10}^{-2}$
$\underset{¯}{\tau }$		0.6
$\overline{\tau }$		0.9

Table 2. Table 2: Parameters considered in Section  7.3 for Algorithm  1 and its initialization.

Parameters		Initialization	Algorithm 1
$J_{\text{cyc}}$		101	11
$J_{{𝗨}_1}^{\left(i\right)}=J_{{𝗨}_2}^{\left(i\right)}$	if $\left(i+1\right)=0\left[\text{mod }{J}_{\text{cyc}}\right]$	0	0
	otherwise	2	5
$J_{ϵ}^{\left(i\right)}$	if $\left(i+1\right)=0\left[\text{mod }{J}_{\text{cyc}}\right]$	1000	1000
	otherwise	0	0
${\xi }_{{𝗨}_{\text{tot}}}$		${10}^{-6}$
${\xi }_{ϵ}$		${10}^{-5}$	${10}^{-5}$
$\zeta$		${10}^{-2}$	${10}^{-2}$
$\underset{¯}{\tau }$		0.2
$\overline{\tau }$		1.1

Table 3. Table 3: Notations used to define the radio interferometric inverse problem.

$N$	Image dimension
$n_a$	Number of antennas
$M=n_a\left(n_a-1\right)/2$	Number of measurements at each time instant
$T$	Number of time instants
$\overline{𝒙}\in {ℝ}^N$	Original unknown sky intensity image
${𝒙}_o\in {ℝ}^N$	Image containing the known brightest sources of $\overline{𝒙}$
$\overline{ϵ}\in {ℝ}^N$	Image containing the unknown sources of $\overline{𝒙}$
$𝗙\in {ℂ}^{N\times N}$	Fourier matrix
${\overline{𝒅}}_{t,\alpha }\in {ℂ}^N$	Original unknown DDE related to antenna $\alpha$ at time instant $t$
$\overline{𝗚}\in {ℂ}^{TM\times N}$	Convolution matrix containing the antenna-based gains
$𝒚\in {ℂ}^{TM}$	Complex visibility vector
$y_{t,\alpha ,\beta }\in ℂ$	Element of the vector $𝒚$ associated with the antenna pair $(\alpha ,\beta )$ at time instant $t$
${𝗬}_t\in {ℂ}^{n_a\times n_a}$	Complex visibility matrix at instant $t$
${𝖸}_{t,(\alpha ,\beta )}\in ℂ$	Element of the matrix ${𝗬}_t$ associated with the antenna pair $(\alpha ,\beta )$ at time instant $t$
${𝗬}_{t,\alpha }\in {ℂ}^{n_a-1}$	Column of ${𝗬}_t$ without its $\alpha$-th element, i.e. ${𝖸}_{t,\alpha }={\left({𝗬}_{t,(\alpha ,\beta )}\right)}_{\underset{\beta \ne \alpha }{\beta \in \{1,\dots ,n_a\}}}$
${\overline{𝗗}}_{t,1}\in {ℂ}^{n_a\times N}$	Matrix containing on each row the reordered Fourier transforms of the DDEs at time instant $t$
${\overline{𝗗}}_{t,2}\in {ℂ}^{n_a\times N}$	Matrix containing on each row the complex conjugates of the Fourier transforms of the DDEs at time instant $t$
$\overline{𝗫}\in {ℂ}^{N\times N}$	Convolution matrix containing the Fourier transform of the original image

Table 4. Table 4: Notations used to define the proposed minimization problem associated with the imaging and calibration problem.

$𝕊$	Support of the Fourier transform of the DDEs
$S$	Cardinal of $𝕊$, i.e. number of non-zero Fourier coefficients of each DDE
$𝔻$	Constraint set for the DDEs defined in eq. (22)
$({\vartheta }_1,{\vartheta }_2)\in ]0,+\infty {[}^2$	Bounds on the amplitude of the Fourier coefficients of the DDEs defined in eq (16)-(17)
$𝔼$	Constraint set for the image defined in eq. (32)
${\overline{𝗨}}_{t,1}\in {ℂ}^{n_a\times S}$	Matrix containing on each row the non-zero reordered Fourier coefficients of the DDEs at time instant $t$
${\overline{𝗨}}_{t,2}\in {ℂ}^{n_a\times S}$	Matrix containing on each row the non-zero reordered Fourier coefficients of the DDEs at time instant $t$
${\overline{𝗨}}_1\in {ℂ}^{T{n}_a\times S}$	Concatenation of matrices ${\left({\overline{𝗨}}_{t,1}\right)}_{1⩽t⩽T}$
${\overline{𝗨}}_2\in {ℂ}^{T{n}_a\times S}$	Concatenation of matrices ${\left({\overline{𝗨}}_{t,2}\right)}_{1⩽t⩽T}$
${\overline{𝒖}}_{t,\alpha ,1}\in {ℂ}^S$	Row $\alpha$ of matrix ${\overline{𝗨}}_{t,1}$
${\overline{𝒖}}_{t,\alpha ,2}\in {ℂ}^S$	Row $\alpha$ of matrix ${\overline{𝗨}}_{t,2}$
${𝗛}_{t,\alpha ,1}\in {ℂ}^{\left(n_a-1\right)\times S}$	Convolution matrix defined in eq. (19)
${𝗛}_{t,\alpha ,2}\in {ℂ}^{\left(n_a-1\right)\times S}$	Convolution matrix defined in eq. (20)
$𝚿\in {ℝ}^{N\times D}$	Sparsity basis used for the image estimation
$𝒢:{ℂ}^{T{n}_a\times S}\times {ℂ}^{T{n}_a\times S}\to {ℂ}^{TM\times N}$	Function defined such that $𝒢({\overline{𝗨}}_1,{\overline{𝗨}}_2)=\overline{𝗚}$
${𝒟}_{t,1}:{ℂ}^{n_a\times S}\to {ℂ}^{n_a\times N}$	Function defined such that ${𝒟}_{t,1}\left({\overline{𝗨}}_{t,1}\right)={\overline{𝗗}}_{t,1}$
${𝒟}_{t,2}:{ℂ}^{n_a\times S}\to {ℂ}^{n_a\times N}$	Function defined such that ${𝒟}_{t,2}\left({\overline{𝗨}}_{t,2}\right)={\overline{𝗗}}_{t,2}$
${𝒟}_{t,\alpha ,1}:{ℂ}^{n_a\times S}\to {ℂ}^{\left(n_a-1\right)\times N}$	Selection operator to build matrix ${𝗛}_{t,\alpha ,1}$ in eq. (19)
${𝒟}_{t,\alpha ,2}:{ℂ}^{n_a\times S}\to {ℂ}^{\left(n_a-1\right)\times N}$	Selection operator to build matrix ${𝗛}_{t,\alpha ,2}$ in eq. (20)
${ℒ}_{t,\alpha ,1}:{ℂ}^{N\times N}\to {ℂ}^{N\times S}$	Selection operator to build matrix ${𝗛}_{t,\alpha ,1}$ in eq. (19)
${ℒ}_{t,\alpha ,2}:{ℂ}^{N\times N}\to {ℂ}^{N\times S}$	Selection operator to build matrix ${𝗛}_{t,\alpha ,2}$ in eq. (20)
${𝒳}_{t,\alpha ,1}:{ℂ}^N\to {ℂ}^{N\times S}$	Function defined such that ${𝒳}_{t,\alpha ,1}\left(𝗙\overline{𝒙}\right)={ℒ}_{t,\alpha ,1}\left(\overline{𝗫}\right)$
${𝒳}_{t,\alpha ,2}:{ℂ}^N\to {ℂ}^{N\times S}$	Function defined such that ${𝒳}_{t,\alpha ,2}\left(𝗙\overline{𝒙}\right)={ℒ}_{t,\alpha ,2}\left(\overline{𝗫}\right)$
$\sigma :{ℂ}^S\to {ℂ}^S$	Bijection defined such that ${\overline{𝒖}}_{t,\alpha ,2}=\sigma \left({\overline{𝒖}}_{t,\alpha ,1}\right)$
$h:{ℝ}^N\times {ℂ}^{T{n}_a\times S}\times {ℂ}^{T{n}_a\times S}\to ℝ$	Data-fidelity term defined in eq. (24), (25) and (26)
$r:{ℝ}^N\to ]-\infty ,+\infty ]$	Regularisation term for the image defined in eq. (31) or (34)
$p:{ℂ}^{T{n}_a\times S}\times {ℂ}^{T{n}_a\times S}\to ]-\infty ,+\infty ]$	Regularisation term for the DDEs defined in eq. (21)

Table 5. Table 5: Notations used for Algorithm  1 .

$\eta \in ]0,+\infty [$	Regularization parameter for the image
$\nu \in ]0,+\infty [$	Regularization parameter for the DDEs
$i\in \mathbb{N}$	Iteration number
$J_{ϵ}^{\left(i\right)}\in \mathbb{N}$	Number of sub-iterations to update the image at iteration $i$
$J_{{𝗨}_1}^{\left(i\right)}\in \mathbb{N}$	Number of sub-iterations to update ${𝗨}_1$ at iteration $i$
$J_{{𝗨}_2}^{\left(i\right)}\in \mathbb{N}$	Number of sub-iterations to update ${𝗨}_2$ at iteration $i$
$J_{\text{cyc}}\in \mathbb{N}$ ($J_{\text{cyc}}⩾2$)	Number of iterations to perform a cycle of Algorithm 1
${\gamma }_{ϵ}^{\left(i\right)}\in [0,+\infty [$	Step-size for the update of the image at iteration $i$
${\gamma }_{{𝒖}_{t,\alpha ,1}}^{\left(i\right)}\in [0,+\infty [$	Step-size for the update of ${𝒖}_{t,\alpha ,1}$ at iteration $i$
${\gamma }_{{𝒖}_{t,\alpha ,2}}^{\left(i\right)}\in [0,+\infty [$	Step-size for the update of ${𝒖}_{t,\alpha ,2}$ at iteration $i$
${\xi }_{{𝗨}_{\text{tot}}}\in ]0,+\infty [$	Stopping criterion for the DDE estimation defined in eq. (38)
${\xi }_{ϵ}\in ]0,+\infty [$	Stopping criterion for the image estimation defined in eq. (39)
$\zeta \in ]0,+\infty [$	Stopping criterion for the global algorithm defined in eq. (40)
$(\underset{¯}{\tau },\overline{\tau })\in ]0,+\infty {[}^2$	Parameters used to determine the initialisation of the image

Equations

$$\mathbf{y}=\mathsf{G}\mathsf{F}\mathbf{x}+\mathbf{b},$$

$$\mathsf{Y}=\mathsf{C}\mathsf{X}{\mathsf{C}}^{\ast }+\mathsf{B},$$

$${\mathsf{C}}_1^{+}=\underset{{\mathsf{C}}_1}{\mathrm{argmin}}\Vert {\mathsf{C}}_1\mathsf{X}{\mathsf{C}}_2^{\ast }-\mathsf{Y}{\Vert }_2^2,$$

$$\underset{\begin{array}{c}\mathbf{ϵ},{\mathsf{U}}_1,{\mathsf{U}}_2\end{array}}{\mathrm{minimize}}h(\mathbf{ϵ},{\mathsf{U}}_1,{\mathsf{U}}_2)+r(\mathbf{ϵ})+p({\mathsf{U}}_1,{\mathsf{U}}_2),$$

$$y_{t,\alpha ,\beta }=\int {\overline{d}}_{t,\alpha }(\mathbf{l}){\overline{d}}_{t,\beta }(\mathbf{l}{)}^{\ast }\overline{x}(\mathbf{l}){\mathsf{e}}^{-2\mathsf{i}\pi {\mathbf{v}}_{t,\alpha ,\beta }\cdot \mathbf{l}}\mathsf{d}\mathbf{l},$$

$$\begin{gathered} y_{t,\alpha ,\beta }=\sum _{n=-N/2}^{N/2-1}{\overline{d}}_{t,\alpha }(n){\overline{d}}_{t,\beta }(n{)}^{\ast }\overline{x}(n){\mathsf{e}}^{-2\mathsf{i}\pi k_{t,\alpha ,\beta }\frac{n}{N}} \\ +b_{t,\alpha ,\beta }, \end{gathered}$$

$$\mathbf{y}=\overline{\mathsf{G}}\mathsf{F}\overline{\mathbf{x}}+\mathbf{b},$$

$$(\forall \mathbf{x}\in {\mathbb{R}}^N)\tilde{h}(\mathbf{x})=\frac{1}{2}\Vert \overline{\mathsf{G}}\mathsf{F}\mathbf{x}-\mathbf{y}{\Vert }_2^2,$$

$$\underset{\begin{array}{c}\mathbf{x}\in {\mathbb{R}}^N\end{array}}{\mathrm{minimize}}\tilde{h}(\mathbf{x})+\eta \tilde{r}(\mathbf{x}),$$

$$\underset{\begin{array}{c}\mathbf{x}\in {\mathbb{R}}^N\end{array}}{\mathrm{minimize}}\tilde{r}(\mathbf{x})\text{ subject to }\Vert \overline{\mathsf{G}}\mathsf{F}\mathbf{x}-\mathbf{y}{\Vert }_2⩽\theta ,$$

$$(\forall \mathbf{x}\in {\mathbb{R}}^N)\tilde{r}(\mathbf{x})=\Vert \mathsf{W}{\mathbf{\Psi }}^{†}\mathbf{x}{\Vert }_1,$$

$${\iota }_{\mathbb{E}}(\mathbf{x})=\{\begin{array}{ll}0,& \text{if }\mathbf{x}\in \mathbb{E},\\ +\infty ,& \text{otherwise.}\end{array}$$

$$(\forall t\in\{1,\ldots,T\})\quad\boldsymbol{\mathsf{Y}}_{t}=\mathcal{S}\big{(}\overline{\boldsymbol{\mathsf{D}}}_{t,1}\,\overline{\boldsymbol{\mathsf{X}}}\,\overline{\boldsymbol{\mathsf{D}}}_{t,2}^{\top}\big{)}+\boldsymbol{\mathsf{B}}_{t},$$

$$(\forall t\in\{1,\ldots,T\})\quad\boldsymbol{\mathsf{Y}}_{t}=\mathcal{S}\big{(}\overline{\boldsymbol{\mathsf{C}}}_{t}\,\widetilde{\boldsymbol{\mathsf{X}}}\,\overline{\boldsymbol{\mathsf{C}}}_{t}^{*}\big{)}+\boldsymbol{\mathsf{B}}_{t},$$

$$\underset{\begin{subarray}{c}{(\boldsymbol{\mathsf{C}}_{t})_{1\leqslant t\leqslant T}\in(\mathbb{C}^{n_{a}\times n_{a}})^{T}}\end{subarray}}{\text{\rm minimize}}\;\;\sum_{t=1}^{T}\frac{1}{2}\|\mathcal{S}\big{(}\boldsymbol{\mathsf{C}}_{t}\,\widetilde{\boldsymbol{\mathsf{X}}}\,\boldsymbol{\mathsf{C}}_{t}^{*}\big{)}-\boldsymbol{\mathsf{Y}}_{t}\|_{2}^{2}.$$

$$\big{|}\real\big{(}\widehat{\overline{d}}_{t,\alpha}(0)\big{)}\big{|}\leqslant\vartheta_{1},\;\big{|}\text{\rm Im}\big{(}\widehat{\overline{d}}_{t,\alpha}(0)\big{)}\big{|}\leqslant\vartheta_{1},$$

$$(\forall k\in\mathbb{S}\setminus\{0\})\quad\big{|}\real\big{(}\widehat{\overline{d}}_{t,\alpha}(k)\big{)}\big{|}\leqslant\vartheta_{2},\;\big{|}\text{\rm Im}\big{(}\widehat{\overline{d}}_{t,\alpha}(k)\big{)}\big{|}\leqslant\vartheta_{2}.$$

$$\underset{\begin{array}{c}({\mathsf{U}}_1,{\mathsf{U}}_2)\in ({\mathbb{C}}^{T{n}_a\times S}{)}^2\end{array}}{\mathrm{minimize}}\overline{h}({\mathsf{U}}_1,{\mathsf{U}}_2)+\nu p({\mathsf{U}}_1,{\mathsf{U}}_2),$$

$$\overline{h}({\mathsf{U}}_1,{\mathsf{U}}_2)$$

$$=\sum _{t=1}^T\sum _{\alpha =1}^{n_a}\frac{1}{2}\Vert {\mathsf{H}}_{t,\alpha ,2}{\mathbf{u}}_{t,\alpha ,2}-{\mathsf{Y}}_{t,\alpha }{\Vert }_2^2,$$

$$\begin{gathered} p({\mathsf{U}}_1,{\mathsf{U}}_2)=\sum _{t=1}^T\sum _{\alpha =1}^{n_a}\frac{1}{2}\Vert {\mathbf{u}}_{t,\alpha ,2}-\sigma ({\mathbf{u}}_{t,\alpha ,1}){\Vert }_2^2 \\ +{\iota }_{\mathbb{D}}({\mathbf{u}}_{t,\alpha ,1})+{\iota }_{\mathbb{D}}({\mathbf{u}}_{t,\alpha ,2}) \end{gathered}$$

$$\displaystyle\mathbb{D}=\Big{\{}\boldsymbol{u}\in\mathbb{C}^{S}$$

$$\text{and }(\forall s\in \{-S/2,\dots ,S/2-1\}\setminus \{0\})$$

$$\displaystyle|\real\big{(}u(s)\big{)}|\leqslant\vartheta_{2},\,|\text{\rm Im}\big{(}u(s)\big{)}|\leqslant\vartheta_{2}\Big{\}}.$$

$$\underset{\begin{array}{c}\underset{({\mathsf{U}}_1,{\mathsf{U}}_2)\in ({\mathbb{C}}^{T{n}_a\times S}{)}^2}{\mathbf{ϵ}\in {\mathbb{R}}^N}\end{array}}{\mathrm{minimize}}h(\mathbf{ϵ},{\mathsf{U}}_1,{\mathsf{U}}_2)+\eta r(\mathbf{ϵ})+\nu p({\mathsf{U}}_1,{\mathsf{U}}_2),$$

$$h(\mathbf{ϵ},{\mathsf{U}}_1,{\mathsf{U}}_2)$$

$$=\Vert \mathcal{G}({\mathsf{U}}_1,{\mathsf{U}}_2)\mathsf{F}({\mathbf{x}}_o+\mathbf{ϵ})-\mathbf{y}{\Vert }_2^2$$

$$\displaystyle\!=\!\!\!\!\sum_{\underset{1\leqslant\alpha\leqslant n_{a}}{1\leqslant t\leqslant T}}\!\frac{1}{2}\|\mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)}\mathcal{X}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{F}}(\boldsymbol{x}_{o}+\boldsymbol{\epsilon})\big{)}\boldsymbol{u}_{t,\alpha,1}-\boldsymbol{\mathsf{Y}}_{t,\alpha}\|_{2}^{2}$$

$$\displaystyle\!=\!\!\!\!\sum_{\underset{1\leqslant\alpha\leqslant n_{a}}{1\leqslant t\leqslant T}}\!\frac{1}{2}\|\mathcal{D}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{U}}_{t,1}\big{)}\mathcal{X}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{F}}(\boldsymbol{x}_{o}+\boldsymbol{\epsilon})\big{)}\boldsymbol{u}_{t,\alpha,2}-\boldsymbol{\mathsf{Y}}_{t,\alpha}\|_{2}^{2}.$$

$${\text{prox}}_{\psi }(\mathbf{z})=\underset{\mathbf{u}\in {\mathbb{R}}^Q}{\mathrm{argmin}}\psi (\mathbf{u})+\frac{1}{2}\Vert \mathbf{u}-\mathbf{z}{\Vert }_2^2.$$

$$\left[{\text{prox}}_{\eta {\ell }_1}(\mathbf{z})\right](n)=\{\begin{array}{ll}-z(n)+\eta ,& \text{if }z(n)<-\eta ,\\ 0,& \text{if }-\eta ⩽z(n)⩽\eta ,\\ z(n)-\eta ,& \text{otherwise},\end{array}$$

$$\{\begin{array}{l}\Re (v(0))=\min \{\max \{\Re (u(0)),-{\vartheta }_1\},{\vartheta }_1\},\\ \text{Im}(v(0))=\min \{\max \{\text{Im}(u(0)),-{\vartheta }_1\},{\vartheta }_1\},\end{array}$$

$$\{\begin{array}{l}\Re (v(s))=\min \{\max \{\Re (u(s)),-{\vartheta }_2\},{\vartheta }_2\},\\ \text{Im}(v(s))=\min \{\max \{\text{Im}(u(s)),-{\vartheta }_2\},{\vartheta }_2\}.\end{array}$$

Full text

††pubyear: 2017††pagerange: Non-convex optimization for self-calibration of direction-dependent effects in radio interferometric imaging– Non-convex optimization for self-calibration of direction-dependent effects in radio interferometric imaging

1 Introduction

Radio interferometry (RI) aims to observe an area of the sky at high angular resolution through an array of antennas. More precisely, each measurement is acquired by a pair of antennas, in the Fourier domain of the intensity image of interest. These measurements consist of complex noisy visibilities, and are related to an undersampled selection of the Fourier coefficients of the intensity image degraded by antenna gains (Thompson et al., 2001).

When the antenna gains are perfectly calibrated, the imaging problem relies on estimating the observed sky intensity from the acquired visibilities. To this aim, several efficient methods have been developed since the early 1970’s, starting with the celebrated CLEAN algorithm (Högbom, 1974; Schwarz, 1978; Thompson et al., 2001). Essentially, this algorithm is a non-linear deconvolution method based on local iterative beam removal. At each iteration, it involves computing the residual dirty image and removing the dirty beam at the point with maximum absolute value in the residual image. The beam removal step is dependent on a loop gain factor, controlling the fraction of the value of the brightest point to be removed. This whole process implicitly assumes sparsity of the sought image. In recent years, several variants of CLEAN have been proposed (e.g. Multi-Scale CLEAN (Cornwell, 2008), Adaptive Scale Pixel- CLEAN (Bhatnagar & Cornwell, 2004)). Furthermore, it is interesting to note that it has been shown recently that the CLEAN methods can be described within the framework of optimization algorithms (Cornwell, 2008; Wiaux et al., 2009a). In particular, it has similarities with the matching pursuit algorithm (Mallat & Zhang, 1993), and it can be seen as a gradient descent algorithm regularized with sparsity constraint on the image (Rau et al., 2009; Carrillo et al., 2014; Onose et al., 2016). The latter comes from the fact that the computation of the residual image is analogous to a gradient step, while the local beam removal by the loop gain factor resembles the thresholding operation to impose sparsity.

Recently, new imaging approaches have been developed to tackle the RI imaging problem leveraging optimization methods (Rockafellar, 1970; Hiriart-Urruty & Lemaréchal, 1993; Boyd & Vandenberghe, 2004), and compressive sensing theory (Candès et al., 2006a, b; Candès, 2006; Donoho, 2006; Baraniuk, 2007). The additive noise on the visibilities being assumed to be independent and identically distributed (i.i.d.) Gaussian, these approaches rely on solving a regularized least squares (LS) minimization problem. The regularization term is then chosen to incorporate prior information on the target intensity image. In particular, in (Wiaux et al., 2009a), the authors exploit the fact that many images are sparse in some dictionary (e.g. image domain, possibly redundant wavelet basis (Mallat, 2009), gradient domain (Rudin et al., 1992; Gilboa & Osher, 2008), etc.) by using an $\ell_{1}$ regularization term. Note that this approach has subsequently been investigated in several works (Wiaux et al., 2009b; Wenger et al., 2010; Li et al., 2011; McEwen & Wiaux, 2011; Carrillo et al., 2012; Dabbech et al., 2015; Onose et al., 2017). Furthermore, in the field of RI imaging, the resulting minimization problem involves large dimensional variables, especially for the big data problems encountered in the era of modern radio telescopes such as the Square Kilometer Array (SKA)111http://www.skatelescope.org/ and LOw Frequency ARray (LOFAR)222http://www.lofar.org/. In this context, it has been shown that proximal methods (Combettes & Pesquet, 2010; Parikh & Boyd, 2013; Komodakis & Pesquet, 2015) are particularly well-suited, since they are designed to solve large scale problems. Notably, a first proximal method for RI imaging has been designed by Carrillo et al. (2014) and developed in the PURIFY package (Pratley et al., 2016), using an simultaneous direction method of multipliers (SDMM) based algorithm (Setzer et al., 2010). This method has been further improved by Onose et al. (2016), using the alternating direction method of multipliers (ADMM) algorithm (Boyd et al., 2011) and a primal-dual forward-backward-based algorithm (Condat, 2013; Vũ, 2013; Pesquet & Repetti, 2014). It is worth noting that the image reconstruction quality obtained by these methods was shown to have the potential to supersede CLEAN for recovering extended emission (Onose et al., 2016).

In practice, however, antenna-based gains are unknown. These gains, related to each antenna of the interferometer, are classified as direction-independent and dependent effects (DIEs and DDEs, respectively). On the one hand, DDEs are time-variable complex gains, different for each antenna and varying across the field of view. They correspond to a multiplication in the image domain, and hence a convolution in the Fourier domain. On the other hand, DIEs stand for a particular case of the DDEs, when the spatial dependency can be ignored at each instant of the integration (i.e. only a scalar complex gain factor is considered for each antenna). Traditionally, the calibration process, i.e. estimation of DIEs/DDEs, focusses on calibration of DIEs only. However, the new generation radio telescopes are being designed to produce orders of magnitude of improvement in both resolution and dynamic range. In order to match the imaging capabilities of these instruments, estimating only the DIEs is no longer sufficient and the DDEs must be incorporated in the calibration process. In general, the calibration problem, either to calibrate for DIEs or independently for DDEs, can be formulated as a non-linear LS minimization problem (Mitchell et al., 2008; Salvini & Wijnholds, 2014; Smirnov & Tasse, 2015). In particular, traditional DIEs self-calibration methods solve this problem using gradient-based techniques such as Levenberg-Marquardt (LM) algorithm (Levenberg, 1944; D.W, 1963). However, Mitchell et al. (2008); Salvini & Wijnholds (2014) pointed out its high computational cost and proposed a fast solver, namely StEFCal, based on an alternating direction implicit (ADI) method. Nevertheless, these methods are only designed to solve the DIE calibration problem, and with an upsurge in the number of unknown parameters, the calibration of DDEs is a highly challenging task. Keeping this in mind, few methods have been developed recently for DDE calibration. On the one hand, Intema et al. (2009) have proposed a DDE calibration method that attempts to iteratively solve and correct for ionospheric phase errors, namely the Source Peeling and Atmospheric Modeling (SPAM) method. Basically, SPAM solves for complex gains towards a number of bright sources, then approximates the DDEs by a sum of low-order Zernike polynomials, fitting the phase component of the solutions. On the other hand, traditional methods for calibration have been improved by the SAGE algorithm in (Yatawatta et al., 2009; Kazemi et al., 2011), while another general framework has been developed by Smirnov & Tasse (2015) solving for non-linear LS, considering the complex Jacobian formalism. Furthermore, a different approach based on non-linear Kalman filter is adopted by Tasse (2014b). Finally, a new calibration scheme based on facet calibration, has been developed in the last years (see e.g. Tasse (2014a); Smirnov & Tasse (2015); van Weeren et al. (2016)). To apply this method, it is assumed that the DDEs are piece-wise constant in the image domain. Therefore, the sky is partitioned into facets, with the facet center determined by the brightest source or the approximate center of a source group. At each step, the calibration solutions are obtained for the facet center, which is then applied over the whole facet. It is worth mentioning that, to the best of our knowledge, this method has only been applied to point sources models of the sky, which at most include few extended sources. Also note that, when this DDE calibration method is combined with an imaging method, the obtained global reconstruction algorithm does not benefit from any global convergence guarantee.

In this paper, we develop a method estimating jointly the unknown intensity image and the DDEs, reconstructing not only point sources but also sophisticated extended sources, and considering a smooth DDE screen applied to all the sources across the image, preventing the need to select calibrator directions. Importantly, it is to our knowledge the first joint imaging and DDE calibration algorithm with convergence guarantees. Moreover, it is important to emphasize that the proposed method can be applied to joint DIE calibration and imaging problem as a particular case, providing a robust convergent self-calibration (selfcal) method. Our approach is inspired by both the imaging techniques using optimization and compressive sensing theories, and the alternating calibration method StEFCal. It aims to minimize a regularized non-linear LS criterion, with respect to both the image and the DDEs. To regularize the problem, in addition to an $\ell_{1}$ regularization term applied to the image (to promote the sparsity of the image, or its sparsity in a given dictionary), we make use of the assumption that the brightest sources of the image are known. Moreover, concerning the DDEs, we assume that they are smooth functions across the field of view, i.e. spatially band-limited. Therefore, we design an algorithm based on recent non-convex optimization techniques, which benefits from convergence guarantees deduced from Bolte et al. (2014); Frankel et al. (2015); Chouzenoux et al. (2016). It involves alternating between the estimation of the image and the DDEs, relying on an iterative structure based on the same optimization toolbox for both image and DDEs estimation, consisting of a gradient step followed by a projection step. In addition to the convergence guarantees, the key point of our algorithm is that it works globally on the whole image and in an automatic manner, without any requirement of partitioning the sky into different directions. Furthermore, it is adapted not only to images of point sources, but also to reconstruct images with extended and complex structures. Finally, the results of our simulations suggest that using our method to jointly estimate the DDEs and the image leads to large improvement of the dynamic range, which is orders of magnitude higher when compared to accounting for DIEs only.

The remainder of the paper is organised as follows. A simple description of our method is provided in Section 2, followed by a complete description given in Sections 3, 4 and 5. More precisely, in Section 3, we describe the imaging inverse problem and the associated minimization problem appearing in the case of a perfectly calibrated telescope array. In Section 4, we introduce the calibration problem to estimate the DDEs from a known sky intensity image. The joint imaging and DDE calibration problem and the proposed alternating method are given in Section 5. Moreover, an alternative DDE calibration method based on the StEFCal algorithm is described in Section 6. The simulations performed and the results obtained thereby using the proposed approach are discussed in Section 7. Finally, in Section 8, we briefly summarize the main contributions presented in this paper and we describe envisaged future work directions. Note that even though the main formal explanations concerning the proposed approach are given in Sections 3, 4 and 5, the reader should be able to understand the main results of Section 7 after the intuitive description of our method given in Section 2.

2 Block-coordinate approach in a nutshell

The joint imaging and DDE calibration method proposed in this article is based on a block-coordinate technique. The mathematical details of this approach are given in Sections 3, 4 and 5, while Section 6 gives details of an alternative StEFCal-based method. Nevertheless, the current section is dedicated to an intuitive description that should enable readers to jump to Section 7 and understand our results and conclusions. Moreover, tables of notations are provided in Appendix A.

The proposed approach is in the spirit of the current imaging techniques based on optimization and compressive sensing theories, and can be seen as a generalization of an existing DIE calibration method, namely StEFCal (Salvini & Wijnholds, 2014).

2.1 Imaging problem

The imaging problem for RI corresponds to an inverse problem (Rau et al., 2009; Carrillo et al., 2014; Onose et al., 2016). In this context, the objective is to estimate an original unknown intensity image ${\boldsymbol{x}}$ from the complex and noisy visibilities $\boldsymbol{y}$ given by

[TABLE]

where $\boldsymbol{\mathsf{F}}$ is the Fourier matrix, ${\boldsymbol{\mathsf{G}}}$ is the matrix containing on each row the antenna-based gain related to each antenna pair acquiring the complex visibilities, and $\boldsymbol{b}$ is a realization of a complex i.i.d. Gaussian additive noise. In the case of a perfectly calibrated antenna array (i.e. $\boldsymbol{\mathsf{G}}$ is known), problem (1) is linear and can be solved efficiently using convex optimization methods (Onose et al., 2016). Basically, these approaches define the estimate of the image as a solution to a regularized LS minimization problem (see Section 3 for more details). In particular, compressive sensing theory ensuring that combining an LS data-fidelity with a regularization term promoting sparsity leads to good reconstruction results (Candès, 2006), $\ell_{1}$-based regularization terms have been extensively used in the last years for RI imaging. Therefore, the underlying minimization problem consists of a sum of an $\ell_{2}$ term and an $\ell_{1}$ term. A simple optimization algorithm to solve this problem is the forward-backward algorithm (Combettes & Wajs, 2005). Basically, it is an iterative method, alternating at each iteration between a gradient step (forward step) on the $\ell_{2}$ function, and a proximity step (backward step) on the $\ell_{1}$ function. Note that in the context of the $\ell_{1}$ norm, the backward step corresponds to a soft-thresholding operation, obliging the small values of the signal to be set to zero, thus promoting sparsity of this signal. Furthermore, as pointed out by Rau et al. (2009); Carrillo et al. (2014); Onose et al. (2016), the forward-backward algorithm is very similar to the CLEAN methods. Indeed, the forward step computes the residual dirty image, while the backward step is used to promote sparsity. Thus, similarities can be seen between these two steps and the major and the minor cycles of CLEAN methods (see e.g. Onose et al. (2016) for details). Nevertheless, it is worth emphasizing that the forward-backward algorithm, and in general optimization algorithms (Combettes & Pesquet, 2010; Komodakis & Pesquet, 2015), can be used to introduce more sophisticated regularization terms. Moreover, it has been shown in the last years that using optimization and compressive sensing theories leads to very competitive results with respect to traditional radio interferometric methods such as CLEAN (Onose et al., 2016; Onose et al., 2017), yet comparisons have mostly been performed on smaller images, and the computational costs so far have not been competitive.

2.2 Calibration problem

In practice, the antenna-based gains contained in ${\boldsymbol{\mathsf{G}}}$ are also unknown and have to be estimated. In the case when they reduce to DIEs and the image is known, Salvini & Wijnholds (2014) have proposed an efficient method solving for the calibration problem, namely the StEFCal method. In the Fourier domain, DIEs consist of only one complex non-zero coefficient centred at the zero spatial frequency. In this context, problem (1) can be rewritten using a matrix formulation:

[TABLE]

where $\boldsymbol{\mathsf{C}}$ is the diagonal matrix containing the non-zero elements of the direction-independent Fourier coefficients, $(.)^{*}$ denotes the complex conjugate of its argument, $\boldsymbol{\mathsf{X}}$ is the matrix containing the Fourier coefficients of the image ${\boldsymbol{x}}$ associated with the frequencies acquired by the antenna pairs, and $\boldsymbol{\mathsf{B}}$ is the matrix formulation of the additive noise. In order to have a bi-linear inverse problem, Salvini & Wijnholds (2014) propose to introduce in eq. (2) matrices $\boldsymbol{\mathsf{C}}_{1}=\boldsymbol{\mathsf{C}}_{2}=\boldsymbol{\mathsf{C}}$, and to solve the non-linear LS problem associated with eq. (2) using an iterative method based on an ADI algorithm. This method alternates at each iteration between the estimation of ${\boldsymbol{\mathsf{C}}}_{1}$ and ${\boldsymbol{\mathsf{C}}}_{2}$, assuming they are independent. More precisely, the update of ${\boldsymbol{\mathsf{C}}}_{1}$ is taken to be the exact minimizer of the LS objective function, i.e.

[TABLE]

while the update of ${\boldsymbol{\mathsf{C}}}_{2}$ is taken to be exactly equal to $\boldsymbol{\mathsf{C}}_{1}$, i.e. $\boldsymbol{\mathsf{C}}_{2}^{+}=\boldsymbol{\mathsf{C}}_{1}^{+}$. As explained by Salvini & Wijnholds (2014), this method converges when the original image is known exactly. However, it has been shown that in practice, it works when only the bright sources of the original image are known. In this case, Salvini & Wijnholds (2014) suggest to combine the StEFCal calibration method with an imaging algorithm in order to estimate more accurately the image once we have an estimation of the DIEs, and to iterate this process. However, this combined DIEs calibration and imaging approach does not benefit from the StEFCal convergence guarantees, and is not adapted to DDE calibration.

2.3 Proposed joint DDE calibration and imaging method

Taking advantage of the above-mentioned methods, we propose a joint imaging and DDE calibration algorithm. In this context, we generalize the inverse problem (2) in order to take into account the DDEs in ${\boldsymbol{\mathsf{C}}}_{1}$ and ${\boldsymbol{\mathsf{C}}}_{2}$. To this aim, we assume that the DDEs are smooth functions across the field of view, i.e. they are spatially band-limited (see Figure 2 in Section 7). Therefore, we propose to reduce drastically the dimensionality of the problem by estimating only the non-zero Fourier coefficients of the DDEs in ${\boldsymbol{\mathsf{C}}}_{1}$ and ${\boldsymbol{\mathsf{C}}}_{2}$, now represented by matrices $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ (see Section 4 for more details). Moreover, we make use of the prior information on the bright sources of the original image. Note that this assumption is common in the context of DDE calibration methods (Yatawatta et al., 2009; van Weeren et al., 2016), and is useful to reduce the ambiguity problems appearing between the image and the DDEs (Smirnov, 2015). More precisely, we assume that the original image $\boldsymbol{x}$ can be split as a sum of two images $\boldsymbol{x}_{o}$ and $\boldsymbol{\epsilon}$, where $\boldsymbol{x}_{o}$ is assumed to be known exactly, while $\boldsymbol{\epsilon}$ has to be estimated (see Figure 3 in Section 7 for an example in the case of point sources). Then, instead of using two different algorithms to estimate the DDEs and the image, respectively, we design a joint framework. In essence, this joint method involves alternating between the estimation of $\boldsymbol{\mathsf{U}}_{1}$, $\boldsymbol{\mathsf{U}}_{2}$ and the faint sources $\boldsymbol{\epsilon}$ contained in ${\boldsymbol{x}}$, using the same algorithmic structure, based on the forward-backward iterations. To do so, we propose to

[TABLE]

where $h$ is the data fidelity term corresponding to a least squares criterion and depending on both the image and the DDEs, $r$ is the regularization function for the image, and $p$ is the regularization function for the DDEs. In particular, $r$ is chosen to constrain the image to be positive and to promote sparsity either directly in the image domain or in a given dictionary. Concerning the DDEs, $p$ is chosen to incorporate constraints on the direction-dependent Fourier coefficients, and to control the similarities between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ (see Section 4 for more details). Note that problem (4) is non-convex with respect to the concatenation of the variables $(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$, but is convex with respect to each of them. In other words, keeping $(\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$ (resp. $(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{2})$ and $(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{1})$) fixed, problem (4) is convex with respect to the variable $\boldsymbol{\epsilon}$ (resp. $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$).

The proposed iterative method is based on a block-coordinate forward-backward approach, solving the regularized LS problem (4), and benefits from the convergence guarantees in Chouzenoux et al. (2016). More precisely, at each iteration, we first estimate approximately the DDEs computing a fixed number of forward-backward iterations, and then estimate approximately the image again using forward-backward iterations. It is important to emphasize that, as a particular case, the proposed approach can be applied to solve the joint DIE calibration and imaging problem.

It is interesting to note that the global structure of the proposed algorithm is very similar to the traditional selfcal method, since they both aim to alternate between the estimation of the gains and the estimation of the image. To illustrate the similarities and differences between the two methods, diagrams depicting the relevant steps of the proposed method and the traditional selfcal method are given in Figure 1 (left) and (right), respectively333The selfcal method presented in this diagram consists in alternating between the StEFCal algorithm described in Section 2.2 and an imaging method, namely CLEAN. However, note that StEFCal can be coupled with other imaging methods, as described in Section 2.1.. One can observe that the global structure of the two diagrams are very similar. However, the structure of the inner-loops are different. In particular, in our method, $J_{\text{cyc}}-1$ inner-loops are performed to estimate approximately the DDEs $(\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$. Within each of these inner-loops, $J_{\boldsymbol{\mathsf{U}}_{1}}$ forward-backward steps are performed to estimate $\boldsymbol{\mathsf{U}}_{1}$, followed by $J_{\boldsymbol{\mathsf{U}}_{2}}$ forward-backward steps to estimate $\boldsymbol{\mathsf{U}}_{2}$. Note that $J_{\boldsymbol{\mathsf{U}}_{1}}$ and $J_{\boldsymbol{\mathsf{U}}_{2}}$ are finite so as to have approximate estimations of $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$, respectively. Then, to complete one global iteration (i.e. a cycle), $J_{\boldsymbol{\epsilon}}$ forward-backward steps are performed to estimate approximately the image $\boldsymbol{\epsilon}$. Note that computing only approximated estimates is very important in practice to ensure the convergence of the overall algorithm. Intuitively, one can understand this point as follows: when the algorithm is initialized with a very poor estimation of the image, it is obvious that estimating completely the DDEs from this incorrect image can be inefficient. Therefore, it is important to control the accuracy of the estimates at each iteration in order to make the overall algorithm converge step by step for the image and the DDEs together. This adopted methodology for approximate estimates is different from the traditional selfcal method, where the problem of estimating the gains (restricted to DIEs) is solved completely before solving the imaging problem, and vice versa. As explained above, this strategy can lead to poor reconstruction results. Furthermore, while the imaging update of our method is very similar to the forward-backward approach described in Section 2.1, the calibration part is different from the StEFCal approach given in Section 2.2. Firstly, using a forward-backward based algorithm allows us to introduce constraints on the DDEs. Moreover, StEFCal is an implicit method, stating directly that the update of $\boldsymbol{\mathsf{C}}_{2}$ is equal to $\boldsymbol{\mathsf{C}}_{1}$, whereas our method updates independently the variables $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$. Thus, in order to constrain them to be equal at convergence, the distance between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ is explicitly controlled in the minimization problem. Finally, StEFCal is only designed for DIE calibration, while our method jointly corrects for DDEs and estimates the image.

Therefore our method can be seen as a generalization of selfcal method, with theoretical convergence guarantees. Note that traditional selfcal method does not benefit from the convergence guarantees of the proposed method since both rely on different algorithmic structures. More details of the proposed algorithm are given in Section 5.

For the sake of completeness, a naive generalization of the StEFCal method for DDE calibration is also proposed in Section 6. Note that in contrast with the StEFCal method for DIE calibration, in the DDE case, it is not enough to consider the LS data-fidelity term, and the minimization problem needs to be regularized. As for the basic StEFCal method for DIE calibration, we propose to combine this naive generalization with an imaging technique. The obtained method consists in alternating between a DDE calibration algorithm (i.e. the naive generalization of StEFCal) and an imaging algorithm (not necessarily based on the CLEAN method). Nevertheless, the obtained method is not guaranteed to converge, fully justifying the development of a more sophisticated approach, described in Sections 3, 4 and 5.

In Section 7, we describe the simulations performed considering both point source images and an image of M31. In both cases, we show the results obtained using our method, reconstructing the DDEs and the image $\boldsymbol{\epsilon}$, and the results obtained with the StEFCal algorithm solving only for the DIEs, and combined with an imaging method based on the forward-backward iterations. In particular, in Section 7.2 we focus on point source images, where the original image $\boldsymbol{x}$ is the sum of three point source images: the known image $\boldsymbol{x}_{o}$, and the unknown images ${\boldsymbol{\epsilon}}_{1}$ and ${\boldsymbol{\epsilon}}_{2}$ (i.e. $\boldsymbol{\epsilon}=\boldsymbol{\epsilon}_{1}+\boldsymbol{\epsilon}_{2}$, see Figure 3). In our simulations, ${\boldsymbol{\epsilon}}_{2}$ is considered to be an astrophysical noise of very low intensity, which can be seen as the confusion noise induced by the resolution limitation of the radio telescope. Thus, the main objective is to estimate ${\boldsymbol{\epsilon}}_{1}$. Our experiments in this context aim to study the behaviour of the proposed approach varying the dynamic range of the image (see Figures 4 and 5), the size of the DDEs (see Figures 6 and 7), the number of antennas considered, and the time dependency of the DDEs (see Figures 9 and 10). In Section 7.3, we describe the results obtained by considering an image of M31 (see Figure 11). In this case, we assume that the brightest part of the image $\boldsymbol{x}_{o}$ is known, and the objective is to estimate the unknown sources ${\boldsymbol{\epsilon}}$. In particular, we perform three experiments, changing the prior information on the bright sources. Figures 12, 13 and 14 show the results obtained for these simulations, giving in the first row the bright sources contained in $\boldsymbol{x}_{o}$ on the left, and the unknown image ${\boldsymbol{\epsilon}}$ on the right. It is worth noting that all the results presented in Section 7 suggest that using our method to jointly estimate the DDEs and the image enables to improve the dynamic range by orders of magnitude compared to accounting for DIEs only. Moreover, in order to present a study as complete as possible, we compared our approach with the naive StEFCal method for DDE calibration presented in Section 6, and we show that our method leads to better reconstruction results.

The following three sections give more details on the proposed method. In particular, Section 3 is dedicated to the imaging part, Section 4 to the calibration part, and the complete algorithm is explained in detail in Section 5. However, at this stage the reader might directly be interested in the simulations and results given in Section 7.

3 Imaging problem formulation

3.1 Observation model

A radio interferometer consists of an array of $n_{a}$ antennas, measuring the radio emission from a given area of the sky. The complex data, namely the visibilities, acquired at instant $t\in\{1,\ldots,T\}$, are determined by the relative position between each antenna pair indexed by $(\alpha,\beta)\in\{1,\ldots,n_{a}\}^{2}$, with $\alpha<\beta$. The baseline components, denoted by \big{(}u_{t,\alpha,\beta},v_{t,\alpha,\beta},w_{t,\alpha,\beta}\big{)}, are measured in units of the observation wavelength. On the one hand, the components in the plane perpendicular to the line of sight \boldsymbol{v}_{t,\alpha,\beta}=\big{(}u_{t,\alpha,\beta},v_{t,\alpha,\beta}\big{)} specify the coordinates of the projected the baseline, $\boldsymbol{v}_{t,\alpha,\beta}=\boldsymbol{v}_{t,\alpha}-\boldsymbol{v}_{t,\beta}$, $\boldsymbol{v}_{t,\alpha}$ (resp. $\boldsymbol{v}_{t,\beta}$) being the position of antenna $\alpha$ (resp. $\beta$) in units of the observation wavelength. On the other hand, the third component $w_{t,\alpha,\beta}$ identifies the coordinate in the line of sight. The sky brightness distribution $\overline{x}$ is described in the same coordinate system, with components $l$, $m$, $n$ and with $\boldsymbol{l}=(l,m)$ and $n(\boldsymbol{l})=\sqrt{1-l^{2}-m^{2}}$, $l^{2}+m^{2}\leqslant 1$. Then, the general measurement equation of non-polarized monochromatic RI imaging, for the antenna pair $(\alpha,\beta)$ at instant $t$, is given by

[TABLE]

where $\overline{\boldsymbol{d}}_{t,\alpha}$ represents a DDE in the image domain related to antenna $\alpha$ at instant $t$.

To recover the original unknown image from incomplete visibility measurements, a discretized version of problem (5) is considered, which corresponds to a sampled version of the continuous sky brightness distribution $\overline{x}$. We denote this sampled original intensity signal by the vector $\overline{\boldsymbol{x}}=(\overline{x}(n))_{-N/2\leqslant n\leqslant N/2-1}\in\mathbb{R}_{+}^{N}$, and the complex visibilities are represented by the vector $\boldsymbol{y}\in\mathbb{C}^{TM}$. More precisely, at each instant $t\in\{1,\ldots,T\}$, an interferometer with $n_{a}$ antennas acquires $M=n_{a}(n_{a}-1)/2$ different measurements in the Fourier domain of the image of interest. Therefore, the visibility $y_{t,\alpha,\beta}\in\mathbb{C}$ measured by the antenna pair $(\alpha,\beta)$ at instant $t$ at the discrete spatial frequency $k_{t,\alpha,\beta}=k_{t,\alpha}-k_{t,\beta}$ can be modelled as the following inverse problem:

[TABLE]

where $\overline{\boldsymbol{d}}_{t,\alpha}=(\overline{d}_{t,\alpha}(n))_{-N/2\leqslant n\leqslant N/2-1}\in\mathbb{C}^{N}$ is the sampled DDE related to antenna $\alpha$, and $\boldsymbol{b}=(b_{t,\alpha,\beta})_{\underset{1\leqslant\alpha<\beta\leqslant n_{a}}{1\leqslant t\leqslant T}}\in\mathbb{C}^{TM}$ is a realization of a complex i.i.d. Gaussian additive noise. Note that, using these notations, the DIEs can be seen as a special case of the DDEs where $\overline{\boldsymbol{d}}_{t,\alpha}=\delta_{t,\alpha}\boldsymbol{1}_{N}$, with $\delta_{t,\alpha}\in\mathbb{C}$ and $\boldsymbol{1}_{N}$ being the unitary vector of dimension $N$.

When the antenna array is perfectly calibrated, i.e. when the DDEs/DIEs are perfectly known, the observations  described in eq. (6) can be rewritten using a matrix formulation

[TABLE]

where $\boldsymbol{\mathsf{F}}\in\mathbb{C}^{N\times N}$ represents the Fourier matrix, and $\overline{\boldsymbol{\mathsf{G}}}\in\mathbb{C}^{TM\times N}$ is the matrix containing on each row the antenna-based gain for the pair $(\alpha,\beta)$. More precisely, each row of $\overline{\boldsymbol{\mathsf{G}}}$, indexed by \big{(}t,\alpha,\beta\big{)}, corresponds to the convolution of the Fourier transforms $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ and $\widehat{\overline{\boldsymbol{d}}^{*}}_{t,\beta}$ of the gains $\overline{\boldsymbol{d}}_{t,\alpha}$ and $\overline{\boldsymbol{d}}^{*}_{t,\beta}$ respectively, centred at the spatial frequency $k_{t,\alpha,\beta}$.

Note that the relevant variables and dimensions associated with the RI inverse problem are summarized in Table 3, given in Appendix A.

3.2 Minimization problem

Considering the inverse problem (7), new methods based on both convex optimization and compressive sensing theory have been developed recently to find an estimate $\boldsymbol{x}^{\star}$ of $\overline{\boldsymbol{x}}$ from the observations $\boldsymbol{y}$ (Wiaux et al., 2009a). In this context, the estimated image is defined as the minimizer of a sum of two functions: the data fidelity term related to the observation model, and the regularization term incorporating a priori information on the target image. In the context of RI imaging, the additive noise being assumed i.i.d. Gaussian, a usual choice for the data fidelity term is the least squares criterion:

[TABLE]

where $\|\cdot\|_{2}$ denotes the standard Euclidean norm. Then, the estimated image can be defined as a solution to:

[TABLE]

where $\widetilde{r}\colon\mathbb{R}^{N}\to]-\infty,+\infty]$ is the regularization function, and $\eta>0$ is the regularization parameter chosen to balance the data fidelity and the regularization terms.

Note that in previous works (Wiaux et al., 2009a; Carrillo et al., 2014; Onose et al., 2016), the authors have chosen to solve a constrained version of problem (9):

[TABLE]

where $\theta>0$ is a parameter to be fixed by the user, associated with the upper bound of the norm of the additive noise in eq. (7). Historically, this formulation has been preferred over problem (9) since in theory, $\theta$ can be determined as a function of the additive noise in eq. (7), whereas $\eta$ is a free parameter. However, due to technical assumptions related to the algorithm proposed in the current work, we will only focus on problem (9) for the remaining of this article (Chouzenoux et al., 2016). Moreover, the constrained problem given in eq. (10) might be unstable in the context of joint imaging and calibration, since, in this case, the bound $\theta$ cannot be properly fixed.

Compressive sensing theory makes use of the assumption that the image $\overline{\boldsymbol{x}}$ has a sparse representation $\boldsymbol{\Psi}^{\dagger}\overline{\boldsymbol{x}}\in\mathbb{C}^{D}$ in a given dictionary $\boldsymbol{\Psi}\in\mathbb{C}^{N\times D}$ (Fornasier & Rauhut, 2011), where $\dagger$ is the transpose conjugate operator. This framework has shown to give good reconstruction results in several application fields such as astronomical remote sensing (Bobin et al., 2008), optical interferometric imaging (Auria et al., 2014; Birdi et al., 2017), and RI imaging (Wiaux et al., 2009a; Li et al., 2011; Garsden et al., 2015; Dabbech et al., 2015). In this context, the regularization function $\widetilde{r}$ in eq. (9) is chosen in order to promote sparsity. Intuitively, the $\ell_{0}$ pseudo-norm, counting the non-zero elements of a signal, is the best choice to promote sparsity (Donoho et al., 1995). More generally $\ell_{p}$ norms, with $0\leqslant p<1$, can be efficiently used (Bouman & Sauer, 1996). However, due to their non-convexity and non-differentiability, these functions are difficult to optimize in practice, and often an $\ell_{1}$ convex relaxation is used (Bect et al., 2004; Donoho, 2006; Figueiredo et al., 2007). Note that (re)weighted $\ell_{1}$ regularization (Candès et al., 2008; Carrillo et al., 2012) can also be used, in order to approximate the $\ell_{0}$ regularization function. The weighted $\ell_{1}$ regularization is defined by

[TABLE]

where $\boldsymbol{\mathsf{W}}\in\mathbb{R}^{D\times D}$ is a weighting matrix, often chosen to be diagonal (if $\boldsymbol{\mathsf{W}}$ is the identity matrix, the usual $\ell_{1}$ regularization term is recovered). Concerning $\boldsymbol{\Psi}$, different dictionaries can be chosen, depending on the nature of the target image. On the one hand, if the image is assumed to have only few non-zero coefficients (e.g. point sources), one can choose $\boldsymbol{\Psi}=\boldsymbol{\mathsf{I}}_{N}$ to be the identity operator in $\mathbb{R}^{N}$. On the other hand, natural images are not necessarily sparse, but can have a sparse representation in other domains. For instance, piece-wise constant images are sparse on the gradient domain, and thus total variation based regularizations can be used for such images (Rudin et al., 1992; Gilboa & Osher, 2008; Wiaux et al., 2010). Another choice for $\boldsymbol{\Psi}$ is to promote sparsity in the wavelet domain, using, e.g., isotropic undecimated wavelet (IUW) transforms (Starck & Murtagh, 1994; Li et al., 2011; Dabbech et al., 2015; Garsden et al., 2015), or a concatenation of wavelet transforms (Vannier et al., 2010; Carrillo et al., 2012).

Note that constraints can also be taken into account in the formulation given in eq. (9). To this aim, we introduce the indicator function of a closed non-empty subset $\mathbb{E}$ of $\mathbb{R}^{N}$, penalizing all the coefficients of the image which are not in $\mathbb{E}$. More formally, we define the indicator function as follows:

[TABLE]

In problem (9), introducing this indicator function in the regularization term is equivalent to constrain $\boldsymbol{x}$ to belong to $\mathbb{E}$, since the objective function cannot be minimized as long as $\boldsymbol{x}\not\in\mathbb{E}$. Constraints, and subsequently indicator functions, can be used to impose positivity of the image by defining $\mathbb{E}=[0,+\infty[^{N}$, or more generally to constrain the amplitude of a signal to stay in a box. They can be used as well to impose sparsity in the case when the sparsity degree $\kappa\in\mathbb{N}^{*}$ is known, by choosing \mathbb{E}=\big{\{}{\boldsymbol{x}\in\mathbb{R}^{N}}~{}\big{|}~{}{\|\boldsymbol{x}\|_{1}\leqslant\kappa}\big{\}}.

4 Calibration problem formulation

4.1 Observation model

In practice, in the context of RI, for every $t\in\{1,\ldots,T\}$, the antenna-based gains $(\overline{\boldsymbol{d}}_{t,\alpha})_{1\leqslant\alpha\leqslant n_{a}}$ described in eq. (6) are unknown and have to be calibrated. To this aim, the visibilities $\boldsymbol{y}$ can be reordered, and a matrix formulation can be used. Let $\overline{\boldsymbol{\mathsf{D}}}_{t,1}\in\mathbb{C}^{n_{a}\times N}$ be the matrix containing on each row $\alpha$ the reordered kernel \big{(}\widehat{\overline{d}}_{t,\alpha}(-k)\big{)}_{k} centred at $k_{t,\alpha}$, and $\overline{\boldsymbol{\mathsf{D}}}_{t,2}\in\mathbb{C}^{n_{a}\times N}$ be the matrix containing on each row $\alpha$ the conjugate $\widehat{\overline{\boldsymbol{d}}}^{\,*}_{t,\alpha}$ of $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ centred at $-k_{t,\alpha}$. Moreover, let \boldsymbol{\mathsf{Y}}_{t}=\big{(}{\mathsf{Y}}_{t,(\alpha,\beta)}\big{)}_{1\leqslant\alpha,\beta\leqslant n_{a}}\in\mathbb{C}^{n_{a}\times n_{a}} be the matrix containing the visibilities, with zero coefficients on its diagonal and such that ${\mathsf{Y}}_{t,(\alpha,\beta)}=y_{t,\alpha,\beta}$ if $\alpha<\beta$, and ${\mathsf{Y}}_{t,(\alpha,\beta)}=y_{t,\beta,\alpha}^{*}$ if $\alpha>\beta$. Then, the inverse problem described in eq. (6) can be equivalently reformulated as follows

[TABLE]

where $(.)^{\top}$ denotes the transpose operation of its argument, $\mathcal{S}\colon\mathbb{C}^{n_{a}\times n_{a}}\to\mathbb{C}^{n_{a}\times n_{a}}$ is the operator setting to zero the coefficients on the diagonal of its argument, $\boldsymbol{\mathsf{B}}_{t}$ is the matrix formulation associated with the realization of the additive noise in eq. (6), and $\overline{\boldsymbol{\mathsf{X}}}\in\mathbb{C}^{N\times N}$ contains on each row/column a shifted version of the Fourier transform of the original image $\widehat{\overline{\boldsymbol{x}}}$. More formally, for every $(k,k^{\prime})\in\{1,\ldots,N\}$, we have $\overline{{\mathsf{X}}}(k,k^{\prime})=\widehat{\overline{x}}(k+k^{\prime}-1)$ if $k+k^{\prime}-1\leqslant N$, and $\overline{{\mathsf{X}}}(k,k^{\prime})=0$ otherwise. This matrix is used to model the convolution operation in the matrix formulation to define $\boldsymbol{\mathsf{Y}}_{t}$. Note that the matrix formulation given by eq. (13) symmetrizes the data, contrary to the inverse problem given in eq. (7). Thus, each measurement appears twice in $\boldsymbol{\mathsf{Y}}_{t}$ (itself and its conjugate). However, we use this formulation in order to simplify the notation in the remaining of the section. Another remark concerns the fact that, in our formulation, we assume that the DDEs are time-dependent. Therefore, when the image is known, $T$ inverse problems of the form of eq. (13) have to be independently solved.

In the particular case when only DIEs are considered, the inverse problem given by eq. (13) can be simplified, and the formulation given by Salvini & Wijnholds (2014) can be used. Indeed, in this case, the Fourier transform $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ of $\overline{\boldsymbol{d}}_{t,\alpha}$ has only one non-zero coefficient centred at the zero spatial frequency, denoted by $\widehat{\overline{d}}_{t,\alpha}(0)$. Thus, each row of $\overline{\boldsymbol{\mathsf{D}}}_{t,1}$ and $\overline{\boldsymbol{\mathsf{D}}}_{t,2}$ contains only one non-zero element, centred in $k_{t,\alpha}$ and $-k_{t,\alpha}$, respectively, and eq. (13) reduces to

[TABLE]

where $\overline{\boldsymbol{\mathsf{C}}}_{t}\in\mathbb{C}^{n_{a}\times n_{a}}$ is the diagonal matrix containing the elements \big{(}\widehat{\overline{d}}_{t,\alpha}(0)\big{)}_{1\leqslant\alpha\leqslant n_{a}} on its diagonal. In this case, $\widetilde{\boldsymbol{\mathsf{X}}}$ is a matrix of size $n_{a}\times n_{a}$ with entries corresponding to the coefficients of $\overline{\boldsymbol{\mathsf{X}}}$ associated with the frequencies selected by the antenna pairs (i.e. the rows associated with the spatial frequency $k_{t,\alpha}$, and the columns associated with the spatial frequency $-k_{t,\beta}$, for every $(\alpha,\beta)\in\{1,\ldots,n_{a}\}^{2}$, with $\beta\neq\alpha$).

4.2 Minimization problem

As described in Section 3.2, for the imaging problem, we assume that the additive noise in the model described by eq. (13) is an i.i.d. Gaussian noise. In this context, several works have proposed to choose the data fidelity term associated with the DIEs or DDEs to be the LS criterion (Salvini & Wijnholds, 2014; Smirnov & Tasse, 2015). In particular, in the case when only DIEs are considered, according to eq. (14), the estimate of $\overline{\boldsymbol{\mathsf{C}}}_{t}$ can be defined as a solution to

[TABLE]

Note that this minimization problem is non-convex. Then, to solve it, Salvini & Wijnholds (2014) proposed to introduce two intermediary variables $\boldsymbol{\mathsf{C}}_{t,1}$ and $\boldsymbol{\mathsf{C}}_{t,2}$ such that $\boldsymbol{\mathsf{C}}_{t}=\boldsymbol{\mathsf{C}}_{t,1}=\boldsymbol{\mathsf{C}}_{t,2}$. In this case, the overall minimization problem remains non-convex, but it becomes convex with respect to $\boldsymbol{\mathsf{C}}_{t,1}$ when $\boldsymbol{\mathsf{C}}_{t,2}$ is fixed, and convex with respect to $\boldsymbol{\mathsf{C}}_{t,2}$ when $\boldsymbol{\mathsf{C}}_{t,1}$ is fixed.

More generally, as described by Smirnov & Tasse (2015), the calibration problem for DDEs can also be written as a LS minimization problem. In particular, by analogy with eq. (15), we can define estimates of the DDEs, stored in the matrices $(\overline{\boldsymbol{\mathsf{D}}}_{t,1},\overline{\boldsymbol{\mathsf{D}}}_{t,2})$, as a solution to a LS minimization problem. However, only the non-zero coefficients of the Fourier transforms of the DIEs appear in eq. (14), while $\overline{\boldsymbol{\mathsf{D}}}_{t,1}$ and $\overline{\boldsymbol{\mathsf{D}}}_{t,2}$ both contain $n_{a}\times N$ unknown coefficients. Therefore, the DDE calibration minimization problem is not only non-convex, but also highly under-determined. Furthermore, the non-convexity becomes even more problematic when the image is also unknown, and the associated minimization problem tends to be more difficult to solve. In order to overcome this difficulty, we make use of the assumption that DDEs are smooth functions of the sky, i.e. spatially band-limited, and we propose to regularize the LS minimization problem accordingly. More precisely, we assume that, for every $t\in\{1,\ldots,T\}$ and $\alpha\in\{1,\ldots,n_{a}\}$, the Fourier transform $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ of ${\overline{\boldsymbol{d}}}_{t,\alpha}$ has a known bounded support, we denote by $\mathbb{S}$ the largest of these supports (i.e. for every $k\not\in\mathbb{S}$, $\widehat{\overline{d}}_{t,\alpha}(k)=0$), and by $S$ its cardinal. In addition, we assume that the amplitude of the coefficients of $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ are bounded, and the amplitude of the central coefficient $\widehat{\overline{d}}_{t,\alpha}(0)$ is much larger than the amplitude of the other coefficients. More formally, we assume that there exist $\vartheta_{1}\gg\vartheta_{2}>0$ such that

[TABLE]

and

[TABLE]

Under this assumption, the DDE calibration problem can be reformulated in order to estimate only the non-zero coefficients of $\overline{\boldsymbol{\mathsf{D}}}_{t,1}$ and $\overline{\boldsymbol{\mathsf{D}}}_{t,2}$. To this aim, we introduce the matrix $\overline{\boldsymbol{\mathsf{U}}}_{t,1}\in\mathbb{C}^{n_{a}\times S}$ (resp. $\overline{\boldsymbol{\mathsf{U}}}_{t,2}\in\mathbb{C}^{n_{a}\times S}$) containing on each row $\alpha$ only the non-zero coefficients of $\overline{\boldsymbol{\mathsf{D}}}_{t,1}$ (resp. $\overline{\boldsymbol{\mathsf{D}}}_{t,2}$). More formally, the $\alpha$-th row of $\overline{\boldsymbol{\mathsf{U}}}_{t,1}$, denoted by $\overline{\boldsymbol{u}}_{t,\alpha,1}\in\mathbb{C}^{S}$, contains \big{(}\widehat{\overline{d}}_{t,\alpha}(-k)\big{)}_{-k\in\mathbb{S}}^{\top}\in\mathbb{C}^{1\times S}, and the $\alpha$-th row of $\overline{\boldsymbol{\mathsf{U}}}_{t,2}$, denoted by $\overline{\boldsymbol{u}}_{t,\alpha,2}\in\mathbb{C}^{S}$, contains \big{(}\widehat{\overline{d}}_{t,\alpha}(k)^{*}\big{)}_{k\in\mathbb{S}}^{\top}\in\mathbb{C}^{1\times S}. Moreover, let $\mathcal{D}_{t,1}\colon\mathbb{C}^{n_{a}\times S}\to\mathbb{C}^{n_{a}\times N}$ and $\mathcal{D}_{t,2}\colon\mathbb{C}^{n_{a}\times S}\to\mathbb{C}^{n_{a}\times N}$ be the functions building respectively matrices $\overline{\boldsymbol{\mathsf{D}}}_{t,1}$ and $\overline{\boldsymbol{\mathsf{D}}}_{t,2}$ from the non-zero coefficient matrices $\overline{\boldsymbol{\mathsf{U}}}_{t,1}$ and $\overline{\boldsymbol{\mathsf{U}}}_{t,2}$. Thus, we have \overline{\boldsymbol{\mathsf{D}}}_{t,1}=\mathcal{D}_{t,1}\big{(}\overline{\boldsymbol{\mathsf{U}}}_{t,1}\big{)} and \overline{\boldsymbol{\mathsf{D}}}_{t,2}=\mathcal{D}_{t,2}\big{(}\overline{\boldsymbol{\mathsf{U}}}_{t,2}\big{)}. Then, we propose to reformulate the minimization problem associated with the DDE calibration as

[TABLE]

where $\boldsymbol{\mathsf{U}}_{1}$ (resp. $\boldsymbol{\mathsf{U}}_{2}$), for all $t\in\{1,\ldots,T\}$, is the concatenation of matrices $\boldsymbol{\mathsf{U}}_{t,1}$ (resp. $\boldsymbol{\mathsf{U}}_{t,2}$), the function $\overline{h}\colon\mathbb{C}^{Tn_{a}\times S}\times\mathbb{C}^{Tn_{a}\times S}\to\mathbb{R}$ is the data fidelity term chosen to be the least squares criterion, and ${p}\colon\mathbb{C}^{Tn_{a}\times S}\times\mathbb{C}^{Tn_{a}\times S}\to]-\infty,+\infty]$ is the regularization term enforcing similarity between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$, and constraining the direction-dependent Fourier coefficients to satisfy conditions (16)-(17).

In addition to the dimensionality reduction, working only on the non-zero coefficients $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ allows to handle independently the DDE associated with each of the antennas. In particular, in this context the data-fidelity term becomes additively separable in terms of antenna specific terms and each time instant. This separability will be a key feature in Section 5 to design a highly parallelizable algorithm. More precisely, we have

[TABLE]

where \boldsymbol{\mathsf{Y}}_{t,\alpha}=\big{(}{\mathsf{Y}}_{t,(\alpha,\beta)}\big{)}_{\underset{\footnotesize\beta\neq\alpha}{\footnotesize\beta\in\{1\ldots,n_{a}\}}}, and, in eq. (19), \boldsymbol{\mathsf{H}}_{t,\alpha,1}=\mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)}\mathcal{L}_{t,\alpha,1}\big{(}\overline{\boldsymbol{\mathsf{X}}}\big{)} with \mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)}\in\mathbb{C}^{(n_{a}-1)\times N} corresponds to the matrix \boldsymbol{\mathsf{D}}_{t,2}=\mathcal{D}_{t,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)} without the row $\alpha$, and $\mathcal{L}_{t,\alpha,1}\colon\mathbb{C}^{N\times N}\to\mathbb{C}^{N\times S}$ is the operator selecting only the columns of $\overline{\boldsymbol{\mathsf{X}}}$ associated with the frequencies corresponding to the support $\mathbb{S}$ centred in $k_{t,\alpha}$. Similarly, in eq. (20), \boldsymbol{\mathsf{H}}_{t,\alpha,1}=\mathcal{D}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{U}}_{t,1}\big{)}\mathcal{L}_{t,\alpha,2}\big{(}\overline{\boldsymbol{\mathsf{X}}}\big{)} where \mathcal{D}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{U}}_{t,1}\big{)}\in\mathbb{C}^{(n_{a}-1)\times N} corresponding to the matrix \boldsymbol{\mathsf{D}}_{t,1}=\mathcal{D}_{t,1}\big{(}\boldsymbol{\mathsf{U}}_{t,1}\big{)} without the row $\alpha$, and $\mathcal{L}_{t,\alpha,2}\colon\mathbb{C}^{N\times N}\to\mathbb{C}^{N\times S}$ being the operator selecting only the columns of $\overline{\boldsymbol{\mathsf{X}}}$ associated with the frequencies $k\in\{-N/2,\ldots,N/2-1\}$ corresponding to the support $\mathbb{S}$ centred in $-k_{t,\alpha}$. Note that the data fidelity term $\overline{h}$ defined by eq. (19)-(20) is non-convex in $(\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$. However, (19) is convex with respect to $\boldsymbol{\mathsf{U}}_{1}$ when $\boldsymbol{\mathsf{U}}_{2}$ is fixed, and (20) is convex with respect to $\boldsymbol{\mathsf{U}}_{2}$ when $\boldsymbol{\mathsf{U}}_{1}$ is fixed.

On the other hand, the regularization function $p$ in eq. (18) is given by

[TABLE]

where $\sigma\colon\mathbb{C}^{S}\to\mathbb{C}^{S}$ is the bijection defined as \sigma(\boldsymbol{u}_{t,\alpha,1})=\big{(}u_{t,\alpha,1}(-s)^{*}\big{)}_{1\leqslant s\leqslant S}, flipping the elements of $\boldsymbol{u}_{t,\alpha,1}$ to obtain the vector \big{(}u_{t,\alpha,1}(-s)\big{)}_{1\leqslant s\leqslant S}, and then taking the complex conjugate of each of its element $u_{t,\alpha,1}(-s)^{*}$, for every $s\in\{1,\ldots,S\}$. Note that we have $\sigma=\sigma^{-1}$. This bijection is used in order to impose the kernels contained in $\boldsymbol{\mathsf{U}}_{1}$ to be equal to the complex conjugate of the kernels contained in $\boldsymbol{\mathsf{U}}_{2}$, up to a reorganization of its coefficients. The set $\mathbb{D}\subset\mathbb{R}^{2S}$ is chosen in order to constrain the DDEs amplitudes to satisfy conditions (16)-(17). Formally, $\mathbb{D}$ is defined by

[TABLE]

5 Proposed joint imaging and DDE calibration

5.1 Proposed blind approach

In this work, we develop a new method for the joint estimation of the original image $\overline{\boldsymbol{x}}$, and the DDEs $\overline{\boldsymbol{d}}_{t,\alpha}$ when both are unknown. To this aim, we consider a global minimization problem, mixing eq. (9) and eq. (18). Moreover, we assume that the original image $\overline{\boldsymbol{x}}$ is the sum of two images, $\boldsymbol{x}_{o}\in\mathbb{R}^{N}$ and $\overline{\boldsymbol{\epsilon}}\in\mathbb{R}^{N}$. The first image contains the brightest sources of the original image, and is assumed to be known as prior information, while the second image contains the other unknown sources and has to be estimated. Note that this assumption is a usual assumption in the context of DDE calibration methods (Yatawatta et al., 2009; van Weeren et al., 2016), and is useful to reduce the ambiguity problems. In particular, the joint calibration and imaging problem being a blind deconvolution problem, any arbitrary modulation of the intensity image can be absorbed in the DDEs without affecting the data in eq. (7). The data are also blind to a common translation of the DDE kernels in the Fourier domain. These two points are mitigated by the assumption that both the intensity and the position of the brightest sources in the sky are known. Then, the global minimization problem can be written as follows:

[TABLE]

where $h\colon\mathbb{R}^{N}\times\mathbb{C}^{Tn_{a}\times S}\times\mathbb{C}^{Tn_{a}\times S}\to\mathbb{R}$ is the data fidelity term corresponding to a least squares criterion and depending on both the image and the DDEs, $r$ is the regularization function for the image (see Section 3.2) and $p$ is the regularization function for the DDEs given by eq. (21).

According to Sections 3 and 4, the data fidelity term has three different formulations, respectively to the image, to the DDEs contained in $\boldsymbol{\mathsf{U}}_{1}$, and those contained in $\boldsymbol{\mathsf{U}}_{2}$. More precisely, we have

[TABLE]

In the formulation (24), $\mathcal{G}\colon\mathbb{C}^{Tn_{a}\times S}\times\mathbb{C}^{Tn_{a}\times S}\to\mathbb{C}^{TM\times N}$ is the function defined such that $\mathcal{G}(\overline{\boldsymbol{\mathsf{U}}}_{1},\overline{\boldsymbol{\mathsf{U}}}_{2})=\overline{\boldsymbol{\mathsf{G}}}$, $\overline{\boldsymbol{\mathsf{G}}}$ being the matrix described in eq. (7), where each row indexed by \big{(}t,\alpha,\beta\big{)} relates to the convolution between the kernels $\boldsymbol{u}_{t,\alpha,1}$ and $\boldsymbol{u}_{t,\beta,2}$. This data fidelity term is used to estimate the image, and we can notice that eq. (24) is convex with respect to $\boldsymbol{\epsilon}$ when $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ are fixed. In the formulation given by eq. (25) (resp. eq. (26)), $\mathcal{X}_{t,\alpha,1}\colon\mathbb{C}^{N}\to\mathbb{C}^{N\times S}$ (resp. $\mathcal{X}_{t,\alpha,2}\colon\mathbb{C}^{N}\to\mathbb{C}^{N\times S}$) is the function defined such that \mathcal{X}_{t,\alpha,1}(\boldsymbol{\mathsf{F}}\overline{\boldsymbol{x}})=\mathcal{L}_{t,\alpha,1}\big{(}\overline{\boldsymbol{\mathsf{X}}}\big{)} (resp. \mathcal{X}_{t,\alpha,2}(\boldsymbol{\mathsf{F}}\overline{\boldsymbol{x}})=\mathcal{L}_{t,\alpha,2}\big{(}\overline{\boldsymbol{\mathsf{X}}}\big{)}), where $\overline{\boldsymbol{\mathsf{X}}}$ is the matrix defined in eq. (13). The data fidelity term given in eq. (25) (resp. eq. (26)) is convex with respect to $\boldsymbol{\mathsf{U}}_{1}$ (resp. $\boldsymbol{\mathsf{U}}_{2}$) when the other variables are fixed. Therefore, even if the function $h$ is non-convex with respect to $(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$, it is convex when each of the three variables is taken independently and the two others are fixed.

Note that the relevant variables and dimensions associated with the joint calibration and imaging minimization problem are summarized in Table 4, and the variables used in the proposed algorithm are summarized in Table 5, given in Appendix A.

5.2 Optimization tools

In order to describe our optimization method to solve the joint calibration and imaging problem given in Section 5.1, we first present in this section the basics of optimization. Note that the definitions given below stand for non-necessarily convex functions. We refer the reader to Rockafellar (1970); Bauschke & Combettes (2011) for an overview in convex optimization, and to Rockafellar & Wets (1997); Mordukhovich (2006) for non-convex optimization.

A function $\psi\colon\mathbb{R}^{Q}\to]-\infty,+\infty]$ is proper if its domain, denoted by $\text{\rm dom}\,\psi$, is non-empty.

Let $\psi$ be a proper, lower-semicontinuous function, bounded from below by an affine function. Its proximity operator (Moreau, 1965) defined at $\boldsymbol{z}\in\mathbb{R}^{Q}$, denoted by $\text{\rm prox}_{\psi}(\boldsymbol{z})$, corresponds to the set of minimizers of $\psi+\frac{1}{2}\|\cdot-\boldsymbol{z}\|_{2}^{2}$, i.e.

[TABLE]

Note that, in the case when $\psi$ is convex, $\text{\rm prox}_{\psi}(\boldsymbol{z})$ reduces to a singleton. Proximity operators are used in optimization algorithms to deal with non-smooth functions, whereas differentiable functions are usually handled by computing gradient steps.

A particular case of proximity operator is the projection operator onto a closed non-empty set $\mathbb{E}$, denoted by $\Pi_{\mathbb{E}}$, obtained when $\psi$ is the indicator function of $\mathbb{E}$. Another well known proximity operator, often used in signal processing, is the soft-thresholding operator (Chaux et al., 2007) corresponding to the case when $\psi$ is the $\ell_{1}$ norm. The $n$-th component of this operator, when applied to a signal $\boldsymbol{z}$, is given by

[TABLE]

where $\eta>0$ is the threshold parameter. Basically, this operator imposes small values of a signal $\boldsymbol{z}$ to be equal to [math] (more precisely, values $z(n)$ such that $|z(n)|\leqslant\eta$), while the other values are reduced accordingly to $\eta$. Note that this operator is a convex relaxation of the hard-thresholding operator (Blumensath & Davies, 2008) corresponding to the case when $\psi$ is the $\ell_{0}$ pseudo-norm.

5.3 Proposed alternated forward-backward algorithm

To solve the joint calibration and imaging problem described in Section 5.1, we propose to exploit the block-variable structure of problem (23). To this aim, we design an iterative alternated minimization algorithm, alternating between the estimation of the unknown image $\overline{\boldsymbol{\epsilon}}$ and the estimation of the DDEs represented by the matrices $\overline{\boldsymbol{\mathsf{U}}}_{1}$ and $\overline{\boldsymbol{\mathsf{U}}}_{2}$. Furthermore, as explained in Sections 3 and 4, due to the constraints and the $\ell_{1}$ regularization term, the global objective function is non-smooth. More precisely, it corresponds to a sum of smooth and non-smooth functions, and we propose to solve it using an alternated forward-backward algorithm (Chouzenoux et al., 2016). This approach combines, at each iteration, a gradient step (forward step) on the Lipschitz-differentiable functions with a proximity step (backward step) on the non-smooth functions.

In the context of problem (23), on the one hand, for the image estimation, at each iteration, a gradient step is performed on the differentiable data fidelity term given by eq. (24) while a proximity step is computed to deal with the non-smooth regularization function $r$. On the other hand, for the DDEs estimation, at each iteration the forward step involves not only the data fidelity term given by eq. (25) (resp. eq. (26)) but also the differentiable part of the function $p$ given by eq. (21), i.e. the squared norm of the difference between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$. Then a projection step is required to ensure that the updates of the DDEs satisfy conditions (16)-(17).

The iterations of this method to solve problem (23) are given in Algorithm 1, and the details are given below.

5.3.1 Number of sub-iterations

At each global iteration $i\in\mathbb{N}$ of the algorithm, for both variables $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$, a given number of sub-iterations $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$ (resp. $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$) are performed in order to obtain an update $\boldsymbol{\mathsf{U}}_{1}^{(i+1)}$ (resp. $\boldsymbol{\mathsf{U}}_{2}^{(i+1)}$) as accurate as desired by the user. Similarly, at iteration $i$, a fixed number of sub-iterations $J_{\boldsymbol{\epsilon}}^{(i)}$ are performed in order to obtain the update $\boldsymbol{\epsilon}^{(i+1)}$ of the image.

Basically the user can choose any values for $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$, $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$ and $J_{\boldsymbol{\epsilon}}^{(i)}$ provided they satisfy the conditions given in Section 5.4 (i.e. they are finite values, and each of the variables $\boldsymbol{\mathsf{U}}_{1}$, $\boldsymbol{\mathsf{U}}_{2}$ and $\boldsymbol{\epsilon}$ need to be estimated at least once during the iteration process). However, to simplify the explanation of our method, we propose to estimate at each iteration $i$ either the DDEs or the image. To this purpose, note that for a given iteration $i$, when $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}=J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}=0$ and $J_{\boldsymbol{\epsilon}}^{(i)}>0$, only the image is updated while the DDEs are kept fixed, i.e. $\boldsymbol{\mathsf{U}}_{1}^{(i+1)}=\boldsymbol{\mathsf{U}}_{1}^{(i)}$ and $\boldsymbol{\mathsf{U}}_{2}^{(i+1)}=\boldsymbol{\mathsf{U}}_{2}^{(i)}$. In contrast, if we choose $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)},J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}>0$ and $J_{\boldsymbol{\epsilon}}^{(i)}=0$, then only the DDEs are updated and the image is kept fixed at iteration $i$, i.e. $\boldsymbol{\epsilon}^{(i+1)}=\boldsymbol{\epsilon}^{(i)}$.

As mentioned in Section 2, Figure 1 (left) shows a diagram giving the relevant steps of the proposed method. It resumes the global structure of the algorithm and describes the number of sub-iterations performed on $\boldsymbol{\mathsf{U}}_{1}$, $\boldsymbol{\mathsf{U}}_{2}$ and $\boldsymbol{\epsilon}$, respectively. In particular, in our method, we define a cycle to be the minimum iterations of Algorithm 1 needed to estimate at least once the DDEs and once the image, and we denote by $J_{\text{cyc}}\geqslant 2$ the global number of iterations of Algorithm 1 to perform a cycle. More precisely, to estimate the DDEs, $(J_{\text{cyc}}-1)$ iterations of Algorithm 1 are performed to estimate $(\boldsymbol{\mathsf{U}}_{1},\boldsymbol{\mathsf{U}}_{2})$. For each of these iterations, $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$ forward-backward steps are computed to estimate $\boldsymbol{\mathsf{U}}_{1}$, followed by $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$ forward-backward steps to estimate $\boldsymbol{\mathsf{U}}_{2}$. Then, to complete one global cycle, $J_{\boldsymbol{\epsilon}}^{(i)}$ forward-backward steps are performed to estimate the image $\boldsymbol{\epsilon}$.

The parameters $(J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)},J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)},J_{\boldsymbol{\epsilon}}^{(i)},J_{\text{cyc}})$ are used in order to alternate more than once between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ before updating the image $\boldsymbol{\epsilon}$. This flexibility is important in practice since it can accelerate the convergence speed of the global algorithm. Furthermore, the reconstruction quality is not highly dependent on the numbers of sub-iterations, provided that they are finite to ensure the convergence of the algorithm (see Section 5.4). In fact, if these parameters are chosen too large, the number of iterations performed on the DDEs and on the image may not be balanced, and the algorithm may reach local minima.

Note that to simplify the readability of Algorithm 1, the parameter $J_{\text{cyc}}$ does not appear explicitly. Then to take it into account, we propose to choose $(J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)},J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)})$ such that

For every (i+1)\neq 0\Big{[}\text{mod }J_{\text{cyc}}\Big{]}, we have $J_{\boldsymbol{\epsilon}}^{(i)}=0$ and $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)},J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}>0$, 2. 2.

For every (i+1)=0\Big{[}\text{mod }J_{\text{cyc}}\Big{]}, we have $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}=J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}=0$ and $J_{\boldsymbol{\epsilon}}^{(i)}>0$.

For instance, if $J_{\text{cyc}}=2$, at iterations $i=0,2,4,\ldots$ of the algorithm, only the DDEs $\boldsymbol{\mathsf{U}}_{1}^{(i)}$ and $\boldsymbol{\mathsf{U}}_{2}^{(i)}$ are updated, computing respectively $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}>0$ and $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}>0$ forward-backward steps, while at iterations $i=1,3,5,\ldots$ of the algorithm only the image $\boldsymbol{\epsilon}^{(i)}$ is updated, computing $J_{\boldsymbol{\epsilon}}^{(i)}>0$ forward-backward steps. If we want to update 2 times more the DDEs, we choose $J_{\text{cyc}}=3$. Then, at iterations $i=0,1,3,4,6,7,\ldots$ of the algorithm, only the DDEs are updated, while the image is only updated at iterations $i=2,5,8,\ldots$. Consequently, increasing $J_{\text{cyc}}$ allows to alternate more often between $\boldsymbol{\mathsf{U}}_{1}^{(i)}$ and $\boldsymbol{\mathsf{U}}_{2}^{(i)}$ before estimating $\boldsymbol{\epsilon}^{(i)}$. However, this number has to be finite in accordance with the convergence conditions given in Section 5.4.

5.3.2 Estimation of DDEs

Algorithm 1 alternates between the estimation of the non-zero coefficients of the Fourier transform of the DDEs (steps 5 to 26) and the estimation of the image (steps 27 to 33). Concerning the estimation of the DDEs, the proposed algorithm updates alternately $\boldsymbol{\mathsf{U}}_{1}$ in steps 5-15 and $\boldsymbol{\mathsf{U}}_{2}$ in steps 16-26. It is worth noting that these steps update in parallel the DDEs for all the antennas and at all time instants. In the following, several additional remarks are made on these steps.

Gradient steps

Steps 10 and 21 correspond to the gradient steps performed on the variables $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$, respectively. These steps being symmetric, we focus on the description of step 10.

In problem (23), there are two differentiable functions corresponding to the data fidelity term given by eq. (24) and the smooth regularization function controlling the distance between $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$. These two functions are complex valued, and as described in (Smirnov & Tasse, 2015), to deal directly with these functions, adapted tools from complex analysis have to be used. However, one can note that, for every $\boldsymbol{u}\in\mathbb{C}^{S}$ and for every $\boldsymbol{\mathsf{H}}\in\mathbb{C}^{(n_{a}-1)\times S}$, we have $\boldsymbol{\mathsf{H}}\boldsymbol{u}=[\boldsymbol{\mathsf{H}}\,|\,{\mathsf{i}}\,\boldsymbol{\mathsf{H}}]\left[\begin{matrix}\real(\boldsymbol{u})\\ \text{\rm Im}(\boldsymbol{u})\end{matrix}\right]$. This relation implies that working with complex valued variables of dimension $S$ is equivalent to deal with real valued variables of dimension $2S$. Thus, to leverage recent results from non-convex optimization, the gradient given in step 10 corresponds to the gradient of the smooth functions involved in eq. (23), rewritten as real-valued functions.

Projection steps

Steps 11 and 22 correspond to projection steps on $\mathbb{D}$, where $\mathbb{D}$ is defined by eq. (22), for each row of variables $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$, respectively. Similarly to the gradient steps, these steps are symmetric. Thus we focus on the description of step 11.

Let \boldsymbol{u}=\big{(}u(s)\big{)}_{-S/2\leqslant s\leqslant S/2-1}\in\mathbb{C}^{S} be a complex vector which can be decomposed as $\boldsymbol{u}=\real(\boldsymbol{u})+{\mathsf{i}}\,\text{\rm Im}(\boldsymbol{u})$. Then the projection of $\boldsymbol{u}$ in $\mathbb{D}$, is given by \Pi_{\mathbb{D}}(\boldsymbol{u})=\big{(}v(s)\big{)}_{-S/2\leqslant s\leqslant S/2-1}, where

[TABLE]

and, for every $s\neq 0$,

[TABLE]

This operation corresponds to a shrinkage of the highest/lowest amplitude values of the Fourier coefficients of the DDEs, in order to satisfy conditions (16)-(17).

5.3.3 Estimation of the image

The proposed Algorithm 1 has the same structure for the estimation of the DDEs and the image, in the sense that forward-backward iterations are used in both cases.

Gradient step

Step 30 corresponds to the gradient step performed on the image. For the imaging problem, this step involves only the least-squares criterion given by eq. (24). Note that, the gradient of this function, computed at the current sub-iterate $\widetilde{\boldsymbol{\epsilon}}^{(i,j)}$ relies on computing the residual image associated with the full current image $(\boldsymbol{x}_{o}+\widetilde{\boldsymbol{\epsilon}}^{(i,j)})$.

Proximity step

Step 31 involves the computation of the proximity operator of the function $r$, corresponding to the regularization term associated with the image. According to remarks made in Section 3.2, we propose to promote sparsity of the global image $\boldsymbol{x}_{o}+\boldsymbol{\epsilon}$ using the $\ell_{1}$ norm and to penalize its negative coefficients thanks to the indicator function of the positive orthant $[0,+\infty[^{N}$.

In the case when the original image $\overline{\boldsymbol{x}}$ is assumed to be sparse, we choose

[TABLE]

where $\mathbb{E}$ is an $N$ dimensional box. More precisely, we know that the original global image $\overline{\boldsymbol{x}}$ corresponds to the sum of $\boldsymbol{x}_{o}$ and $\overline{\boldsymbol{\epsilon}}$. Let $\mathbb{I}\subset\{-N/2,\ldots,N/2-1\}$ be the set of indices giving the exact position of the known bright sources contained in the image $\boldsymbol{x}_{o}$. In other words, $\mathbb{I}$ is the support of $\boldsymbol{x}_{o}$. Then, we have, for every $n\in\mathbb{I}$, $x_{o}(n)=\overline{x}(n)$. In this context, we oblige the coefficients $\epsilon(n)$ such that $n\in\mathbb{I}$ to be equal to [math], while we constrain the other coefficients to be positive. More formally, we define

[TABLE]

In this context, the proximity operator $\text{\rm prox}_{\gamma_{\boldsymbol{\epsilon}}^{(i)}\eta r}$ of $\gamma_{\boldsymbol{\epsilon}}^{(i)}\eta r$ at $\boldsymbol{z}\in\mathbb{R}^{N}$, for $\gamma_{\boldsymbol{\epsilon}}^{(i)}>0$ and $\eta>0$, corresponds to the projection onto $\mathbb{E}$ composed with a soft-thresholding operator (see eq. (28)), i.e.

[TABLE]

As described in Section 3.2, when the unknown original image $\overline{\boldsymbol{x}}$ is not sparse in its domain, we can decompose it in a given dictionary $\boldsymbol{\Psi}$ in order to promote sparsity of its coefficients $\boldsymbol{\Psi}^{\dagger}\overline{\boldsymbol{x}}$. In this case, the regularization term is chosen to be

[TABLE]

The proximity operator of the above function does not have an explicit formula, and it has to be computed approximately with sub-iterations (e.g. using the Dual forward-backward algorithm (Combettes et al., 2011)).

More generally, note that the regularization function $r$ can be non-convex but, for convergence purposes, needs to be semi-algebraic444Semi-algebraicity is a property satisfied by a wide class of functions, which means that their graph is a finite union of sets defined by a finite number of polynomial inequalities. In particular, it is satisfied for the functions mentioned in eq. (31) and eq. (34)..

5.4 Convergence results

The proposed method is based on a block-coordinate forward-backward algorithm developed by Chouzenoux et al. (2016) (other versions can be found in Frankel et al. (2015); Bolte et al. (2014)). In particular, we can deduce from this paper the convergence guarantees of Algorithm 1 given below.555Note that the algorithms, and the associated convergence results, presented in Chouzenoux et al. (2016); Frankel et al. (2015); Bolte et al. (2014) are designed to minimize real valued functions, and thus are not directly applicable to solve problem (23). However, to design Algorithm 1, problem (23) has been reformulated to deal with real valued functions.

Assume that, for each iteration $i\in\mathbb{N}$ of Algorithm 1, the step-size $\gamma_{\boldsymbol{\epsilon}}^{(i)}$, associated with the image update, and, for every $(t,\alpha)\in\{1,\ldots,T\}\times\{1,\ldots,n_{a}\}$ the step-sizes $\gamma_{\boldsymbol{u}_{t,\alpha,1}}^{(i)}$, $\gamma_{\boldsymbol{u}_{t,\alpha,2}}^{(i)}$, associated with the DDEs updates, satisfy

[TABLE]

where $\|.\|_{S}$ denotes the spectral norm, and $\boldsymbol{\mathsf{H}}_{t,\alpha,1}^{(i)}$, $\boldsymbol{\mathsf{H}}_{t,\alpha,2}^{(i)}$, and $\boldsymbol{\mathsf{G}}^{(i)}$ are the matrices built in steps 8, 19 and 28, respectively. Moreover, assume that the number of sub-iterations determined by $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$, $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$, and $J_{\boldsymbol{\epsilon}}^{(i)}$ are finite and that there exists $(i^{\prime},i^{\prime\prime},i^{\prime\prime\prime})\in\mathbb{N}^{3}$ such that $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i^{\prime})}>0$, $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i^{\prime\prime})}>0$ and $J_{\boldsymbol{\epsilon}}^{(i^{\prime\prime\prime})}>0$. This latter condition ensures that all the different variables (i.e. the DDEs and the image) will be updated at least once during the iterative process. Then, the sequence of iterates \Big{(}\boldsymbol{\epsilon}^{(i)},\boldsymbol{\mathsf{U}}_{1}^{(i)},\boldsymbol{\mathsf{U}}_{2}^{(i)}\Big{)}_{i\in\mathbb{N}} generated by Algorithm 1 converges to a critical point \Big{(}\boldsymbol{\epsilon}^{\star},\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star}\Big{)} of the objective function in eq. (23). Moreover, the objective function value is decreasing along the iterations.

Finally, Chouzenoux et al. (2016, Cor. 3.1) ensures that Algorithm 1 converges globally to a solution to problem (23) if it is initialized not too far from this solution. This result shows the importance of the initialization of algorithms solving non-convex problems. Therefore, the initialization of Algorithm 1 is explained in detail in the following section.

5.5 Initialization

As explained above, problem (23) being non-convex, the final solution of the DDEs and the image obtained by Algorithm 1 is likely to be sensitive to the initialization used for these unknowns. Therefore, we propose an empirical approach to initialize the proposed algorithm, leading in practice to a stable global method, as attested by the simulations presented in Section 7. This initialization is composed of two steps described below.

Step 1.

The initialization of the image and the DDEs is done during the same process. It relies on the assumption that the DDEs are band-limited, and that the zero spatial frequency coefficients (the DIEs) of the DDEs $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ are much larger than the other coefficients. Thus, we propose to first jointly estimate the image $\overline{\boldsymbol{\epsilon}}$ and only the DIEs $\widehat{\overline{d}}_{t,\alpha}(0)$ using Algorithm 1. This involves solving problem (23) considering a support of size $1\times 1$ for the DDEs, ignoring the other Fourier coefficients \big{(}\widehat{\overline{d}}_{t,\alpha}(k)\big{)}_{k\in\mathcal{S}\setminus\{0\}} modelled in the inverse problem (6). In other words, even if the inverse problem of interest includes DDEs, the initialization step we propose ignores this fact and solves the problem taking into account only DIEs in the algorithm.

For this first run of Algorithm 1, DIEs $\widehat{\overline{d}}_{t,\alpha}(0)$ are initialized randomly, and the image $\overline{\boldsymbol{\epsilon}}$ is chosen to be the zero image.

At the end of this first step, we obtain an estimation of the DIEs, denoted by $\check{\boldsymbol{\mathsf{U}}}_{1}$ and $\check{\boldsymbol{\mathsf{U}}}_{2}$, and an estimation of the unknown image $\overline{\boldsymbol{\epsilon}}$, denoted by $\check{\boldsymbol{\epsilon}}$.

Step 2.

Once we have a first estimation of the DIEs from Step 1, we use it to initialize the variables $\widehat{{\boldsymbol{d}}}_{t,\alpha}^{(0)}$ to estimate the full DDEs. More precisely, in Algorithm 1, $\widehat{{d}}_{t,\alpha}^{(0)}(0)$ is taken to be the DIE estimated in Step 1, while the other coefficients are chosen randomly.

Concerning the initialization of the image, preliminary simulations have shown that the first estimation obtained from Step 1 is not necessarily good. In some cases, the estimation only contains artefacts and will slow down the convergence of the Algorithm, but in worse cases, it can even lead to poor reconstruction results. This can be explained by the fact that in Step 1 only the DIEs are estimated and not the DDEs, which is not enough to reconstruct high dynamic range images.

Therefore we examine the information content of the reconstructed image from Step 1 by comparing it with the considered prior image $\boldsymbol{x}_{o}$.

Intuitively, if the amplitudes of the unknown sources belonging to $\overline{\boldsymbol{\epsilon}}$ are close to the amplitudes of the known sources in $\boldsymbol{x}_{o}$, considering DIEs may lead to a first good estimation of $\overline{\boldsymbol{\epsilon}}$. However, in the case when the amplitudes of the unknown sources are much lower than the amplitudes of the known sources, DIEs will not lead to an accurate estimation. For this purpose, we design an automatic selective method, based on the comparison of the energy contained in both images, in the Fourier domains.

Let $\check{\boldsymbol{\mathsf{G}}}=\mathcal{G}(\check{\boldsymbol{\mathsf{U}}}_{1},\check{\boldsymbol{\mathsf{U}}}_{2})$ be the convolution matrix associated with the DIEs estimated in Step 1 ($\mathcal{G}$ being defined in eq. (24)). Then, we compare the energy, i.e. the $\ell_{2}$ norm, of $\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}}$ with that of $\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o}$, corresponding respectively to the energies of images $\check{\boldsymbol{\epsilon}}$ and $\boldsymbol{x}_{o}$ taking into account only the measured frequencies from their Fourier coefficients. In the context of radio interferometry imaging, $u-v$ coverages tend to measure more frequently low frequencies in the Fourier domain, containing mainly the information on the approximations of the image. Then, comparing the energy of $\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}}$ with the energy of $\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o}$ gives a weighted information, without taking into account details of the image which are known to be not recovered by DIEs.

Empirically, we propose to distinguish three different cases:

If $\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}})<\underline{\tau}\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o})$, where $\underline{\tau}\geqslant 0$, we assume that the unknown sources have an amplitude much smaller than the amplitude of the known bright sources, and thus cannot be reconstructed correctly only by considering DIEs. Then we consider that the image estimated with Step 1 cannot be trusted. In this case, we discard this estimate, and we initialize Algorithm 1 for joint imaging with DDE calibration considering $\boldsymbol{\epsilon}^{(0)}$ to be a zero image. 2. 2.

If $\underline{\tau}\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o})\leqslant\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}})\leqslant\overline{\tau}\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o})$, with $\overline{\tau}\geqslant\underline{\tau}$, we assume that the unknown sources and the known bright sources have similar amplitudes, and thus $\overline{\boldsymbol{\epsilon}}$ could be reconstructed approximately only considering DIEs. Then we consider that the image estimated in Step 1 can be trusted. In this case, we use $\check{\boldsymbol{\epsilon}}$ to initialize $\boldsymbol{\epsilon}^{(0)}$ in Algorithm 1 for joint imaging with DDE calibration. 3. 3.

If $\overline{\tau}\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o})\leqslant\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}})$, we consider that the image estimated with Step 1 can be trusted but may contain some artefacts due to DIEs. Thus, we threshold the estimate $\check{\boldsymbol{\epsilon}}$, until having $\overline{\tau}\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}_{o})\geqslant\ell_{2}(\check{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\check{\boldsymbol{\epsilon}})$. Then, $\boldsymbol{\epsilon}^{(0)}$ in Algorithm 1 for joint imaging with DDE calibration is taken to be the thresholded image.

Remarks.

Note that the proposed initialization process is totally automatic. Moreover, we observed during preliminary simulations that it leads to stable results, in the sense that considering, for the same experiment, different random initializations both for DIEs in Step 1 and DDEs in Step 2 leads to similar reconstruction results. In addition, it is interesting to emphasize that Step 1 of the proposed initialization is very similar to perform the StEFCal method developed by Salvini & Wijnholds (2014), coupled with an imaging method based on the forward-backward algorithm (Combettes & Wajs, 2005). In particular, we observed numerically that these two methods give very similar results. However, from an algorithmic view point, the calibration part of these methods are different since Step 1 relies on an alternated forward-backward iterations, while StEFCal is based on an alternating direction implicit method. Moreover, unlike Step 1 which can be seen as a joint DIE calibration and imaging method, StEFCal is not an imaging method. Finally, note that if only DIEs are considered in the observation model, then Step 1 can be seen as a selfcal method, with convergence guarantees (and is enough to jointly estimate the DIEs and the image).

5.6 Additional stopping criteria

In order to avoid useless gradient and proximal steps computation, we introduce additional stopping criteria in Algorithm 1. Firstly, as explained earlier, the algorithm performs $J_{\text{cyc}}-1$ updates of the DDEs for a fixed image, where $J_{\text{cyc}}$ is chosen in order to avoid a complete convergence of the DDEs when the image is only known approximately. However, in practice, when the global algorithm has already performed a large number of iterations, the DDE updates can stabilise in less than $J_{\text{cyc}}-1$ iterations. Therefore, we introduce parameter $\xi_{\boldsymbol{\mathsf{U}}_{\text{tot}}}>0$, corresponding to an additional stopping criterion such that, for every $i\in\mathbb{N}$ satisfying (i+1)\neq 0\;\Big{[}\text{mod }J_{\text{cyc}}\Big{]},

[TABLE]

Then, if this criterion is satisfied, the algorithm stops estimating the DDEs during the current cycle, and estimates directly the image.

Similarly, we introduce parameter $\xi_{\boldsymbol{\epsilon}}>0$, corresponding to an additional stopping criterion to avoid useless updates of the image in the inner-loop given by steps 29-32 in Algorithm 1. Then, for every $i\in\mathbb{N}$ satisfying (i+1)=0\;\Big{[}\text{mod }J_{\text{cyc}}\Big{]}, we consider

[TABLE]

Then, if this criterion is satisfied, the algorithm stops estimating the image and starts a new cycle, estimating again the DDEs.

Finally, to stop the global algorithm at convergence, we use a stopping criterion depending on the objective function. More precisely, we introduce $\zeta>0$, and we define, for every $i\in\mathbb{N}$,

[TABLE]

where $\varphi^{(i+1)}$ is the value of the objective function to be minimized in eq. (23) evaluated at the current update $(\boldsymbol{\epsilon}^{(i+1)},\boldsymbol{\mathsf{U}}_{1}^{(i+1)},\boldsymbol{\mathsf{U}}_{2}^{(i+1)})$.

Note that these stopping criteria are used in the algorithm to avoid useless calculus when the iterates are almost constant along the iterations. They can be seen as additional optional parameters to $J_{\text{cyc}}$, $J_{\boldsymbol{\mathsf{U}}_{1}}$, $J_{\boldsymbol{\mathsf{U}}_{2}}$, and $J_{\boldsymbol{\epsilon}}$ to ensure the global convergence of Algorithm 1.

5.7 Computational cost

The objective of this section is to provide information on the computational cost of the proposed algorithm. As described in Section 5.3.1, at each iteration $i\in\mathbb{N}$, either the DDEs or the image are updated, and one global cycle of the algorithm performs $J_{\text{cyc}}-1$ iterations updating only the DDEs and one iteration updating only the image. Therefore, we describe below the behaviour of the first cycle of Algorithm 1, which has to be repeated until convergence, in order to provide rough information on the computational cost of the proposed method.

The first $J_{\text{cyc}}-1$ iterations of the cycle are used to update the DDEs. To this aim, a discrete Fourier transform of the global image $(\boldsymbol{x}_{o}+\boldsymbol{\epsilon})$ is performed in step 4. Then, for every $i\in\{0,\ldots,(J_{\text{cyc}}-1)-1\}$, $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ are updated sequentially. Firstly, $\boldsymbol{\mathsf{U}}_{1}$ is updated by implementing $Tn_{a}$ parallel steps consisting in (i) the construction of the observation matrices of dimension $(n_{a}-1)\times S$ (described in step 8), and (ii) $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$ forward-backward steps on a variable of dimension $S$ (described in steps 10-11). Then $\boldsymbol{\mathsf{U}}_{2}$ is updated similarly (using steps 19, 21, and 22). Therefore, in total for the update of the DDEs, $Tn_{a}$ parallel steps are performed, each computing $2(J_{\text{cyc}}-1)$ times (i), and $\sum_{i=0}^{J_{\text{cyc}}-2}(J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}+J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)})$ times (ii).

The last iteration $i=J_{\text{cyc}}-1$ of the first cycle is dedicated to the update of the image. In this case, the convolution matrix is built in step 28, followed by $J_{\boldsymbol{\epsilon}}^{(i)}$ forward-backward steps on a variable of dimension $N$. These forward-backward steps are described in steps 30-31. Note that they are the heaviest steps of the algorithm. Indeed, the forward step consists in performing a discrete Fourier transform on the global image, followed by an inverse discrete Fourier transform. Concerning the backward step, if the regularisation described in eq. (31) is considered, then step 31 reduces to a soft-thresholding operation. However, in the case when eq. (34) is considered, then step 31 is computed using to sub-iterations.

Note that the description of the algorithm given here does not take into account the stopping criteria described in Section 5.6. In particular, the condition given in eq. (38) is used to perform at maximum $J_{\text{cyc}}-1$ iterations estimating the DDEs, while the condition given in eq. (39) is used to perform at maximum $J_{\boldsymbol{\epsilon}}^{(i)}$ sub-iterations estimating the image. Examples of numbers of sub-iterations $J_{\boldsymbol{\epsilon}}^{(i)}$, $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$, $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$ and $J_{\text{cyc}}$ and of stopping criteria values are given in the simulation sections, in Tables 1 and 2.

Finally, it is important to emphasize that for the sake of simplicity the current implementation uses fast Fourier transform (FFT) for the Fourier transforms and do not take into account the gridding and degridding steps needed to perform non-uniform Fourier transforms (Fessler & Sutton, 2003). These steps need to be incorporated in the convolution matrices built in steps 8-19 (for the DDEs estimation) and 28 (for the image estimation).

6 StEFCal for DDE calibration

For the sake of completeness, we propose in this section a naive generalization of the StEFCal method developed by Salvini & Wijnholds (2014) for DDE calibration. However, the obtained method is not guaranteed to converge, fully justifying the development of a more sophisticated approach, described in the previous sections.

In the particular case when only DIEs appear in the inverse problem (7) (or its equivalent matrix formulation), the calibration problem can be solved efficiently using the StEFCal algorithm. This method consists in estimating the zero spatial frequency elements \big{(}\widehat{\overline{d}}_{t,\alpha}(0)\big{)}_{1\leqslant\alpha\leqslant n_{a}} of the DIEs, by solving problem (15) using an alternating direction implicit (ADI) method. Note that this approach is a DIE calibration method, which can be combined with an imaging algorithm, and is not adapted for the DDE calibration problem.

The StEFCal approach used in the previous sections can be naively generalized to solve the DDE calibration problem described in eq. (18), with $\nu=0$. In this case, the ADI algorithm defines at each iteration the update $\boldsymbol{\mathsf{U}}_{1}^{+}$ of $\boldsymbol{\mathsf{U}}_{1}$ by setting, for every $t\in\{1,\ldots,T\}$ and $\alpha\in\{1,\ldots,n_{a}\}$,

[TABLE]

where \mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)} and $\boldsymbol{\mathsf{Y}}_{t,\alpha}$ are defined in eq. (19), $\mathcal{X}_{t,\alpha,1}$ is the function defined in eq. (25), and $\boldsymbol{\epsilon}^{\star}$ is the current estimation of the image. Then, $\boldsymbol{\mathsf{U}}_{2}$ is updated by taking \boldsymbol{u}_{t,\alpha,2}^{+}=\sigma\big{(}\boldsymbol{u}_{t,\alpha,1}^{+}\big{)}, where $\sigma$ is the bijection defined in eq. (21). Note that this approach requires the inversion of the matrix \Big{(}\mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)}\mathcal{X}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{F}}(\boldsymbol{x}_{o}+\boldsymbol{\epsilon}^{\star})\big{)}\Big{)}^{\dagger}\mathcal{D}_{t,\alpha,2}\big{(}\boldsymbol{\mathsf{U}}_{t,2}\big{)}\mathcal{X}_{t,\alpha,1}\big{(}\boldsymbol{\mathsf{F}}(\boldsymbol{x}_{o}+\boldsymbol{\epsilon}^{\star})\big{)} at each iteration. Unfortunately, this matrix is not necessarily invertible, and the StEFCal method cannot be directly applied to the DDE calibration problem.

We propose to modify the minimization problem solved by Salvini & Wijnholds (2014) by introducing a regularization term. More precisely, we propose to solve problem (18) where $p$ is given by:

[TABLE]

In this context, the update of $\boldsymbol{\mathsf{U}}_{1}$ is given by, for every $t\in\{1,\ldots,T\}$ and $\alpha\in\{1,\ldots,n_{a}\}$,

[TABLE]

which admits an explicit solution. In eq. (43), the parameter $\nu$ not only is a regularization parameter, but also allows to have a more stable minimization problem to be solved at each iteration. Moreover, as mentioned earlier, to estimate both the DDEs and the image, one can combine a calibration method to an imaging algorithm. Therefore, the proposed generalized StEFCal based method repeats the following two steps until convergence:

Estimate the DDEs from the current estimate $\boldsymbol{\epsilon}^{\star}$ of the original image: Update iteratively $\boldsymbol{\mathsf{U}}_{1}$ using eq. (43), from the current estimation of the image, and $\boldsymbol{\mathsf{U}}_{2}$ such that, for every $t\in\{1,\ldots,T\}$ and $\alpha\in\{1,\ldots,n_{a}\}$, \boldsymbol{u}_{t,\alpha,2}^{+}=\sigma\big{(}\boldsymbol{u}_{t,\alpha,1}^{+}\big{)}. Stop at convergence. 2. 2.

Estimate the unknown faint sources $\overline{\boldsymbol{\epsilon}}$ of the image from the current estimate $(\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star})$ of the DDEs : Minimize $h(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star})+\eta r(\epsilon)$, where $h$ is given by eq. (24), and $r$ is defined by eq. (31) or eq. (34).

In step 2, a simple method to minimize $h(\boldsymbol{\epsilon},\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star})+\eta r(\epsilon)$ is to use the forward-backward algorithm (Combettes & Wajs, 2005). It corresponds to the steps 30-31 in Algorithm 1, when $J_{\boldsymbol{\epsilon}}^{(i)}\to+\infty$.

The main differences between this StEFCal-based approach and our method described in Section 5 are the following:

Our method allows to take into account constraints in the amplitude of the Fourier coefficients of the DDEs, thanks to the projection steps 11 and 22 in Algorithm 1. Note that these constraints could be taken into account also using the ADI algorithm by using sub-iterations, which may slow down drastically the algorithm. 2. 2.

The ADI method is an implicit algorithm, meaning that it takes advantage from the symmetric formulation of eq. (18) by obliging the second variable $\boldsymbol{\mathsf{U}}_{2}$ to be exactly equal to the first variable $\boldsymbol{\mathsf{U}}_{1}$, up to a transformation. This symmetry is not directly used in the structure of our algorithm, since in Algorithm 1 both $\boldsymbol{\mathsf{U}}_{1}$ and $\boldsymbol{\mathsf{U}}_{2}$ are updated alternately using forward-backward iterations. 3. 3.

While the StEFCal algorithm is initially only a calibration method, which can be combined with an imaging algorithm, the proposed method described in Section 5 is directly designed to estimate jointly the image and the DDEs. Therefore, at each iteration, our method does not require a prefect estimation of both DDEs and the image. The accuracy of these estimates is controlled, at each iteration $i\in\mathbb{N}$ of Algorithm 1, by the parameters $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$, $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$, $J_{\text{cyc}}$, and $J_{\boldsymbol{\epsilon}}^{(i)}$. Thus, in our method, the global estimations are obtained at the convergence of the algorithm, as described in Section 5.4, while the proposed naive generalization of StEFCal does not have any convergence guarantee. In particular, we will show in Section 7.2.4 that the results obtained with the StEFCal based method are less accurate than those obtained using our algorithm.

7 Simulations and results

7.1 Simulation setting

To show the performance of the proposed method, we computed numerical experiments, implemented in MATLAB666The MATLAB code is available at https://basp-group.github.io/joint-img-dde-calibration/.. We consider $n_{a}$ randomly distributed antennas, where each antenna pair acquires $T$ measurements, $n_{a}$ and $T$ being specified for the different simulations performed in this section. Note that Earth rotation is incorporated to track the $u-v$ positions of each resulting baseline by considering a time interval of 10 hours777The $u-v$ tracks are simulated using the code available at http://www.astro.umd.edu/$\sim$cychen/MATLAB/ASTR410/uvAnd Beams.html. Moreover, in order to simplify the experiments, we use discrete versions of the associated $u-v$ coverages. This is done by considering the nearest discrete $u-v$ position of each antenna. This approximation is adopted in order to avoid to introduce gridding and degridding steps to model the non-uniform FFT in this first presentation of the proposed DDE selfcal algorithm.

For every $\alpha\in\{1,\ldots,n_{a}\}$, the DDEs kernel $\widehat{\overline{\boldsymbol{d}}}_{t,\alpha}$ associated with antenna $\alpha$ at instant $t\in\{1,\ldots,T\}$ is simulated randomly in the Fourier domain and the support $\mathbb{S}$ of the direction-dependent Fourier coefficients is assumed to be known exactly. In our approach, we do not assume that calibration transfer has been performed but that both DIEs and DDEs belong to some complex valued intervals centred in [math]. In practice, in the proposed algorithm, these intervals are handled by considering the set $\mathbb{D}$ defined by eq. (22) and parametrized by $(\vartheta_{1},\vartheta_{2})$. In particular, we choose $\vartheta_{1}=1.5+\theta$ and $\vartheta_{2}=\theta$, where $\theta>0$ is a parameter depending on the standard deviation $\upsilon>0$ used to generate the DDEs. In other words, we impose to our algorithm to estimate DDEs such that the real and imaginary parts of their zero spatial frequency coefficients (i.e. the DIEs) belong to $[-1.5-\theta,+1.5+\theta]$, while the other coefficients belong to $[-\theta,+\theta]$. Note that it is reasonable to consider small values for the higher order spatial frequencies since they represent direction-dependent variations in the gain across the field of view with respect to the mean gain. Therefore, in our simulations, we generate the DIEs such that both their real and imaginary parts have values around $\pm 1$, while we generate the other coefficients of the DDEs such that their real and imaginary parts have values in the neighbourhood of [math] (in both cases with standard deviation $\upsilon$). In the context when calibrator transfer has been performed, normalized DIEs are obtained. Then, in our method we can introduce stronger constraints obliging the reconstructed DIEs to belong to a complex neighbourhood of $1+0\mathsf{i}$. More formally, these constraints can be taken into account by our algorithm by imposing that the real part of the DIEs belongs to $[1-\theta,1+\theta]$ and their imaginary part belongs to $[-\theta,+\theta]$.

An example of one considered DDE is shown in Figure 2, in image and Fourier spaces. For instance, in the context of the SPAM method (Intema et al., 2009), such DDEs can be seen as ionospheric phase screen.

An extensive study has been performed by considering a wide variety of cases, both varying parameters for the calibration and the imaging part. In Section 7.2, we consider a simulated sky of point sources, and we investigate the impact of the range of the image, the size of the support of the Fourier direction-dependent kernels along with the associated standard deviation $\upsilon$, and the time dependency. Finally, in Section 7.3 we show that our method can be used as well to reconstruct complex structured images such as the image of M31. For all these cases, results will be compared in terms of different metrics. In particular, we will consider the signal to noise ratio (SNR) of the reconstructed image $\boldsymbol{\epsilon}^{\star}$ with respect to the original $\overline{\boldsymbol{\epsilon}}$, which is expressed as

[TABLE]

and the $\ell_{2}$ norm of the residual images, obtained either considering the estimated DDEs or the true DDEs. More precisely, we will consider the weighted $\ell_{2}$ norm of the residual images \|\boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N}, where $\boldsymbol{x}^{\star}=\boldsymbol{x}_{o}+\boldsymbol{\epsilon}^{\star}$ corresponds to the global estimated image and $\boldsymbol{\mathsf{G}}^{\star}$ corresponds to the estimated DDEs, i.e. $\boldsymbol{\mathsf{G}}^{\star}=\mathcal{G}(\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star})$ ($\mathcal{G}$ being defined in eq. (24)). Similarly, we will consider the weighted $\ell_{2}$ norm of the residual images \|\boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N}, where the matrix $\overline{\boldsymbol{\mathsf{G}}}$ introduced in eq. (7) corresponds to the original DDEs. Note that the latter metric is used in order to point out the ambiguity problems which can appear between $\boldsymbol{x}^{\star}$ and $\boldsymbol{\mathsf{G}}^{\star}$, i.e. due to the inverse problem formulation given in eq. (7), parts of $\boldsymbol{x}^{\star}$ can be absorbed in $\boldsymbol{\mathsf{G}}^{\star}$, and vice versa.

In the two simulation settings, we show results obtained using both the proposed method and the StEFCal algorithm (Salvini & Wijnholds, 2014). Note that the latter has been designed to solve only for DIEs, not for DDEs, and without estimating the image. Thus, the results shown in our experiments are obtained combining the StEFCal algorithm, correcting for DIEs, with an imaging algorithm based on the forward-backward algorithm (Combettes & Wajs, 2005), considering the same regularization terms as used in our method, and the same associated parameters. We alternate between the DIEs estimation and the image estimation until stabilization of the overall process. In the following, we will refer to this method using the notation StEFCal-FB. It is interesting to note that in our experiments, the StEFCal-FB method gives results of the same quality as those obtained after the first step of the initialization of our algorithm, where only the DIEs are considered, as described in Section 5.5.

Finally, as proposed in Section 6, we perform simulations on images of point sources using the proposed naive generalization of the StEFCal algorithm to correct for the DDEs. Discarding for one instant the important problem of the lack of convergence guarantee, we show that this generalization of StEFCal combined with a forward-backward algorithm already gives good reconstruction results, although they are not as good as those obtained using our method.

7.2 Sparse images with point sources

In this first part, we consider simulated sky images of size $128\times 128$, consisting of point sources, where each source corresponds to a small 2D Gaussian kernel of size $3\times 3$. As described in Section 5.1, we consider that the original image can be decomposed as $\overline{\boldsymbol{x}}=\boldsymbol{x}_{o}+\overline{\boldsymbol{\epsilon}}$, where $\boldsymbol{x}_{o}$ is known, while $\overline{\boldsymbol{\epsilon}}$ is unknown and has to be estimated jointly with the DDEs. More precisely, in all the experiments of this section, we assume that $\boldsymbol{x}_{o}$ corresponds to 10 bright sources generated randomly such that their total energy888The energy ${\mathsf{E}}(\boldsymbol{x})$ of a vector $\boldsymbol{x}\in\mathbb{C}^{N}$ is defined as

${\mathsf{E}}(\boldsymbol{x})=(\sum_{n=1}^{N}|x(n)|^{2})^{1/2}$. is ${\mathsf{E}}(\boldsymbol{x}_{o})=10$. An example of $\boldsymbol{x}_{o}$ is shown in the first image in Figure 3. Moreover, the sources in $\overline{\boldsymbol{\epsilon}}$ are generated randomly with two intensity levels, i.e. $\overline{\boldsymbol{\epsilon}}=\overline{\boldsymbol{\epsilon}}_{1}+\overline{\boldsymbol{\epsilon}}_{2}$. In all the experiments of this section, we consider $\overline{\boldsymbol{\epsilon}}_{2}$ having an energy ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{2})=10^{-5}$ and to be composed of $200$ sources (see the third image of Figure 3 for an example), while, depending on the experiments, we vary the energy and the number of sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$. In particular, an example of $\overline{\boldsymbol{\epsilon}}_{1}$ containing $10$ sources with ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$ is shown in the second image of Figure 3, and the global image $\overline{\boldsymbol{x}}=\boldsymbol{x}_{o}+\overline{\boldsymbol{\epsilon}}_{1}+\overline{\boldsymbol{\epsilon}}_{2}$ in this case is given in the right image of Figure 3 (in log scale). In this context, $\overline{\boldsymbol{\epsilon}}_{2}$ has sources with intensity $\sim 10^{7}$ weaker than the sources belonging to $\boldsymbol{x}_{o}$. Thus, $\overline{\boldsymbol{\epsilon}}_{2}$ is considered in our simulations as astrophysical noise, and we do not focus on its estimation. Therefore, our main objective is to find an estimate of $\overline{\boldsymbol{\epsilon}}_{1}$, reconstructing its faint sources.

In this section, the original image being sparse in its domain, we use eq. (31) as regularization term for the image, and we fix the regularization parameters $\eta=10^{-3}$ for the image, $\nu=2000$ for the initialization estimating the zero spatial frequency coefficients of the DDEs (i.e. the DIEs), and $\nu=1000$ for the main algorithm estimating DDEs. Moreover, the different parameters for Algorithm 1 are given in Table 1. As described in Section 5.3, on the one hand, at each iteration $i\in\mathbb{N}$ such that (i+1)\neq 0\,\Big{[}\text{mod }J_{\text{cyc}}\Big{]}, the algorithm alternates between the estimation of $\boldsymbol{\mathsf{U}}_{1}^{(i)}$ and $\boldsymbol{\mathsf{U}}_{2}^{(i)}$ computing $J_{\boldsymbol{\mathsf{U}}_{1}}^{(i)}$ (resp. $J_{\boldsymbol{\mathsf{U}}_{2}}^{(i)}$) forward-backward sub-iterations, while the image is not updated. On the other hand, if (i+1)=0\,\Big{[}\text{mod }J_{\text{cyc}}\Big{]}, then the algorithm does not update the DDEs but computes $J_{\boldsymbol{\epsilon}}^{(i)}$ sub-iterations to update the image $\boldsymbol{\epsilon}^{(i)}$. Moreover, $\xi_{\boldsymbol{\mathsf{U}}_{\text{tot}}}$, $\xi_{\boldsymbol{\epsilon}}$, and $\zeta$ are the stopping criteria defined in Section 5.6.

Finally, in this section the considered images consisting of point sources, in addition to the SNR and the $\ell_{2}$ norm of the residual images, we will consider a success rate counting the positions of the sources in $\overline{\boldsymbol{\epsilon}}_{1}$ which are recovered successfully. More precisely, we compare the positions of the brightest sources found in the estimate ${\boldsymbol{\epsilon}}^{\star}$ with the true positions of sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$. Then, we deduce from this comparison a percentage of success for the position estimation of the faint sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$.

7.2.1 Dynamic range of the image

Simulation settings.

One of the main advantages of taking DDEs into account, and not only DIEs, is to be able to reconstruct more accurately images with high dynamic range (Smirnov, 2015; Smirnov & Tasse, 2015). To illustrate this assertion, and in order to show the performance and the limits of our method, in this first test, we vary the energy of $\overline{\boldsymbol{\epsilon}}_{1}$ and the number of sources it contains, while all other parameters are kept fixed. More precisely, we consider the measurements made by $n_{a}=200$ antennas for a single time interval $T=1$, and considering DDEs Fourier kernels with support size $S=7\times 7$, and standard deviation set to be $\upsilon=0.05$ (see Section 7.1). Then, we perform simulations for 10, 50 and 90 sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$, and for each of these values, we consider ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})\in\{10,1,10^{-1},10^{-2},10^{-3}\}$. Note that the only case when the intensity of sources belonging to $\boldsymbol{x}_{o}$ is of the same order as those belonging to $\overline{\boldsymbol{\epsilon}}_{1}$ appears if $\overline{\boldsymbol{\epsilon}}_{1}$ has 10 sources with global energy ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=10$.

Simulation results.

Figure 4 gives the results, as a function of ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})$ using different metrics, obtained for an average over 10 realizations varying the antenna distribution, the random images, and the DDEs coefficients. The first graph shows the SNR of the reconstructed image $\boldsymbol{\epsilon}^{\star}$ with respect to the original image $\overline{\boldsymbol{\epsilon}}$. In the second graph, we consider a success rate counting the positions of the sources in $\overline{\boldsymbol{\epsilon}}_{1}$ which are recovered successfully. These first two graphs show that, whatever the number of considered sources in $\overline{\boldsymbol{\epsilon}}_{1}$ and the value of ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})$, our method outperforms the StEFCal-FB method, showing the importance of reconstructing the DDEs. More precisely, on the one hand, our method is able to recover $100\%$ of sources’ positions for ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})\in\{0.1,1,10\}$, with respective SNR values of $10.7$, $16.9$ and $13.4$ dBs, independently of the number of sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$. It is important to emphasize that, in these three cases, the reconstruction quality depends mainly on the energy of the image and not on its intensity level. Indeed, if $\overline{\boldsymbol{\epsilon}}_{1}$ contains $90$ sources with global energy equals to ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=10$, the intensity of sources belonging to $\overline{\boldsymbol{\epsilon}}_{1}$ is almost two times lower than the intensity of the sources in $\boldsymbol{x}_{o}$. On the other hand, the StEFCal-FB obtains approximately $100\%$ only in the case ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=10$ with SNR value of $11$ dB (resp. $9.1$ dB and $8.6$ dB) when 10 (resp. 50 and 90) sources belong to $\overline{\boldsymbol{\epsilon}}_{1}$. Note that this latter corresponds to the case when the sources in $\overline{\boldsymbol{\epsilon}}_{1}$ are almost of same intensity as the sources in $\boldsymbol{x}_{o}$. Concerning our method, the only case when the number of sources considered in $\overline{\boldsymbol{\epsilon}}_{1}$ gives different reconstruction results is when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=10^{-2}$. Therefore, this case can be seen as a limit case, since for ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=10^{-3}$ our method has a success rate for position recovery of $0\%$. It suggests that using the proposed method we are able to improve the dynamic range by at least three orders of magnitude compared to accounting for DIEs only. The third graph of Figure 4 shows the weighted $\ell_{2}$ norm of the residual images \|\boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N}. Here, $\boldsymbol{x}^{\star}=\boldsymbol{x}_{o}+\boldsymbol{\epsilon}^{\star}$ corresponds to the global estimated image (obtained either using our method or StEFCal-FB), and $\boldsymbol{\mathsf{G}}^{\star}$ corresponds either to the DDEs estimated with our method, i.e. $\boldsymbol{\mathsf{G}}^{\star}=\mathcal{G}(\boldsymbol{\mathsf{U}}_{1}^{\star},\boldsymbol{\mathsf{U}}_{2}^{\star})$ ($\mathcal{G}$ being defined in eq. (24)), or to the DIEs obtained using the StEFCal-FB method. Similarly, the last graph in Figure 4 shows the weighted $\ell_{2}$ norm of the residual images \|\boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N}, where the matrix $\overline{\boldsymbol{\mathsf{G}}}$ introduced in eq. (7) corresponds to the original DDEs. These two last graphs show again the advantage of reconstructing the full DDEs instead of considering only DIEs. In particular, the third graph shows that our method gives a smaller norm of the residual images (of order $10^{-3}$ vs. $10^{-1}$ with StEFCal-FB), taking into account both the estimated image and the estimated DDEs. However, the fourth graph suggests that there is an ambiguity error between $\boldsymbol{\mathsf{G}}^{\star}$ and $\boldsymbol{x}^{\star}$, mainly when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})\in\{1,10\}$ (i.e. corresponding to the cases when the sources in $\overline{\boldsymbol{\epsilon}}$ are of similar range as the sources in $\boldsymbol{x}_{o}$) leading to important errors in the residual images when the true DDEs are considered.

The above global observations can be observed as well by Figure 5 showing images obtained considering 50 sources in $\overline{\boldsymbol{\epsilon}}_{1}$, with energy ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})\in\{0.01,0.1,1\}$. The first two rows consider the case when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=1$, showing the original image $\overline{\boldsymbol{\epsilon}}$ in the first column. The second column gives the reconstructed images obtained using the proposed method (first row) and the StEFCal-FB method (second row). For this realization, our method finds $100\%$ of the positions of the sources in $\overline{\boldsymbol{\epsilon}}_{1}$ and the estimate has an SNR equal to $18.9$ dB, while the StEFCal-FB method finds $90\%$ of the sources with a reconstructed image of SNR $=-3.22$ dB. This difference of SNR in the reconstruction can be understood by noticing that the background of the estimated image with StEFCal-FB is very noisy due to the non-estimation of the DDEs, which is not the case using our method. The third and the fourth columns correspond respectively to the residual images |\boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}| and |\boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}| obtained with our method (first row) and the StEFCal-FB method (second row). In both cases, one can notice the errors of the residual images obtained using our method and considering $\boldsymbol{\mathsf{G}}^{\star}$ (resp. $\overline{\boldsymbol{\mathsf{G}}}$) are $100$ (resp. $10$) times smaller than using the StEFCal-FB method. However, as already noticed in Figure 4, there is obviously an ambiguity error between $\boldsymbol{\mathsf{G}}^{\star}$ and $\boldsymbol{x}^{\star}$ leading to larger errors when the true observation matrix $\overline{\boldsymbol{\mathsf{G}}}$ is considered. The third and fourth rows show similar results, in the cases when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$ and ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.01$, respectively. However, since in these cases the StEFCal-FB method has a success rate of $0\%$ for finding the positions of faint sources, only the results obtained using our method are shown. In the case ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$ presented in the third row, the proposed method recovers all the positions of the faint sources, and the reconstructed image has an SNR of $12.32$ dB. Similarly, the reconstructed image presented in the fourth row, in the case ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.01$, our method recovers $90\%$ of the faint sources positions, and the SNR of the estimate is equal to $6.5$ dB. Thus, when the SNR is low, one can still observe visually that our method leads to good reconstruction results. Furthermore, as observed in Figure 4, the $\ell_{2}$ norm of the residual images decreases with the energy of $\overline{\boldsymbol{\epsilon}}_{1}$.

Ambiguity.

As mentioned in the previous paragraph, Figure 4 and Figure 5 emphasize the ambiguity in the reconstruction of the DDEs and the reconstruction of the image. This ambiguity can be observed by comparing the $\ell_{2}$ errors of the third and fourth graphs of Figure 4, and comparing the error images shown in the third and fourth graphs of Figure 5. Indeed, one can observe that the residual images obtained with both the estimated image and DDEs (i.e. \boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}) have much smaller $\ell_{2}$ norm and amplitude than the residual images obtained considering the true DDEs (i.e. \boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}). Note that in particular, in Figure 5, the most significant ambiguities are associated with the brightest sources positions. This behaviour has already been observed in previous works, even in the DIE calibration case, and is known as “ghost sources” (see e.g. Grobler et al. (2014)).

7.2.2 Size of direction-dependent Fourier kernels

Simulation settings.

A second set of simulations have been performed considering images with point sources, investigating the size $S$ of the support of the Fourier direction-dependent kernels, and the standard deviation $\upsilon$ defined in  Section 7.1. In this section, we consider the measurements obtained with $n_{a}=200$ antennas for a single time interval $T=1$. Moreover, we generate images $\overline{\boldsymbol{\epsilon}}_{1}$ such that they contain $10$ sources, and we fix ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$. We perform simulations considering direction-dependent Fourier kernels with support size $S\in\{3\times 3,7\times 7,11\times 11,15\times 15\}$, and, for each of these values, we consider the standard deviation to be $\upsilon\in\{0.005,0.01,0.05,0.1\}$ (see Section 7.1). The objective of these tests is to analyse the estimation of the unknown image by the proposed method, depending on the importance of the direction-dependent Fourier kernels. Typically, for a given $\upsilon$, larger is the support size $S$, the more is the number of unknowns for DDEs. Similarly, for a given $S$, the DDEs will affect much more the observed sky image for larger values of $\upsilon$, leading to a more important need for DDE calibration. Hence, for both these cases, accurate estimation of both the image and the DDEs becomes difficult.

Simulation results.

A comparison of the results obtained using the proposed approach and the StEFCal-FB method is presented in Figure 6. As in Figure 4, the four graphs represent respectively: the SNR of the estimated images $\boldsymbol{\epsilon}^{\star}$, the success rate determining the percentage of recovered sources positions in $\overline{\boldsymbol{\epsilon}}_{1}$, and the $\ell_{2}$ norm of the residual images obtained considering either the estimated DDEs (third graph) or the true DDEs (fourth graph). As expected, the difficulty of reconstructing accurately the image and the DDEs increases with $S$ and $\upsilon$. In particular, we can observe that the only case when StEFCal-FB estimates correctly $100\%$ of the sources’ positions is for $S=3\times 3$ and $\upsilon=0.005$. However, even in this case, our method leads to a better estimation of $\overline{\boldsymbol{\epsilon}}$ as attested by the SNR graph. Indeed, in this case the images reconstructed with our method have an averaged SNR of $26.9$ dB, while it is equal to $5.7$ dB with StEFCal-FB. In general, we can observe that for a support of Fourier direction-dependent kernels of sizes $S=3\times 3$ and $S=7\times 7$, our method recovers $100\%$ of the sources in $\overline{\boldsymbol{\epsilon}}_{1}$. Moreover, for $S=3\times 3$ the reconstructed images have a stable average SNR, reaching $20.2$ dB in the worst case when $\upsilon=0.1$, and the residual images have small errors for all $\upsilon$ considered. This leads to the conclusion that even when direction-dependent Fourier kernels have very small amplitude, it is not enough to estimate only the DIEs, and our method improves significantly the quality of the reconstructed image. In the cases when $S=11\times 11$ and $S=15\times 15$, our method is able to detect correctly $100\%$ of the positions of sources in $\overline{\boldsymbol{\epsilon}}_{1}$ until $\upsilon=0.01$ and $\upsilon=0.05$, respectively. This shows that the proposed method is very stable even considering large Fourier direction-dependent kernels, impacting significantly the observations. However, it seems that there is a breaking point when $S=15\times 15$ and $\upsilon=0.1$. In this case our algorithm is unstable and can lead to very poor reconstructions as attested by the four graphs.

In addition, images associated with these simulations are provided in Figure 7. The first two rows show the images associated with the simpler case when $S=3\times 3$ and $\upsilon=0.005$, showing the original image $\overline{\boldsymbol{\epsilon}}$ in the first column. The other images in these rows consists of the results obtained using our method (first row) and the StEFCal-FB algorithm (second row). Visually, these three images are very similar and both methods recover the exact positions of the 10 sources. However, the quality of the image recovered with our method (first row, second column) is significantly better, since it has an SNR equals to $28.1$ dB, against $6.7$ dB for StEFCal-FB (second row, second column). This shows that even with very small DDEs, the details of the image are estimated more accurately with the proposed method, estimating the full DDEs. The third and fourth columns show the residual images obtained considering the estimated DDEs and the true DDEs, respectively. Note that in both cases, our method leads to residual images with two orders of magnitude lower than the residual images obtained using StEFCal-FB. Similarly, the third and fourth rows show the results obtained with our method for the cases $S=11\times 11$ and $\upsilon=0.05$, and $S=15\times 15$ and $\upsilon=0.01$, respectively. The third row considers a case with very large DDEs for both the support size and the standard deviation. For this simulation, the SNR is equal to $8.26$ dB, and we can observe few small artefacts in the reconstructed image. The fourth row gives results considering a very large support size $S$ but with a small standard deviation. In this case, the SNR of the reconstructed image is equal to $12.35$ dB, and we can observe that there are errors in the amplitude estimation but the sources positions are recovered exactly. Note that in both cases presented in third and fourth rows, the errors in the residual images are very small. Moreover, accordingly to the third and fourth graphs of Figure 6, it seems that the ambiguity error is decreasing when both the standard deviation $\upsilon$ and the support size $S$ are decreasing.

7.2.3 Time dependency and antennas distribution

Simulation settings.

As explained in Section 4.1, we consider a general model for DDEs where they can be time dependent. In the simulations presented in Sections 7.2.1 and 7.2.2, we have considered snapshot imaging, i.e. the measurements are acquired at a single instant corresponding to $T=1$. However, in practice, the interferometric measurements can be acquired over a period of time at several instants, thereby increasing the sampling of the Fourier plane, and the number of different DDEs to be recovered. Moreover, at each instant, the sampling of the Fourier plane is determined by the number of antennas $n_{a}$ considered. Therefore, we implement our algorithm for different values of $T$ and $n_{a}$, fixing all the other parameters. In this experiment, we consider images with $10$ sources in $\overline{\boldsymbol{\epsilon}}_{1}$ and global energy ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$, while for the direction-dependent Fourier kernels we fix $S=7\times 7$ with $\upsilon=0.05$. Then, we perform simulations with $n_{a}\in\{50,100,200\}$, considering $T\in\{1,10,20\}$ for $n_{a}\in\{100,200\}$, and $T\in\{10,20\}$ for $n_{a}=50$. Note that we have not considered the case when $T=1$ and $n_{a}=50$ since this configuration leads to a highly under-sampled $u-v$ coverage giving very poor results. Figure 8 shows two examples of discrete versions of $u-v$ coverages obtained in the cases when (left) $n_{a}=200$ and $T=20$, and (right) $n_{a}=50$ and $T=20$. They corresponds to 2D masks selecting the measured frequencies in the Fourier domain of the image of interest.

Simulation results.

A comparison of the results obtained using the proposed approach and the StEFCal-FB method is presented in Figure 9. As in the previous sections, the first graph represents the SNR of the estimated image $\boldsymbol{\epsilon}^{\star}$, the second graph gives the percentage of recovered sources positions in $\overline{\boldsymbol{\epsilon}}_{1}$, and the third and fourth graphs show the errors of the residual images obtained considering respectively the estimated DDEs and the true DDEs. Considering the reconstruction quality of $\boldsymbol{\epsilon}^{\star}$, we can observe that our method detects correctly $100\%$ of the sources positions for all the considered cases, which is not the case for StEFCal-FB. Note that increasing the number of integrations per antenna pair leads to better reconstruction results in terms of SNR of the image. However, it also increases the number of DDE unknowns which are uncorrelated due to time dependency. As a consequence, the error in the residual images increases with $T$.

These results are corroborated by the images shown in Figure 10. The first row corresponds to the case when $n_{a}=200$ and $T=20$, with images recovered considering the $u-v$ coverage associated with Figure 8 (left). It gives, from left to right, the original unknown image $\overline{\boldsymbol{\epsilon}}$, its reconstruction obtained using our method with corresponding SNR equal to $25.58$ dB, and the residual images obtained using either the estimated DDEs, or using the true DDEs. The second row gives similar results for the case when $n_{a}=50$ and $T=20$, corresponding to the $u-v$ coverage associated with Figure 8 (right). In this case, the SNR of the reconstructed image is equal to $10$ dB.

7.2.4 StEFCal for DDE calibration

As described in Section 6, the StEFCal-FB method can be generalized to solve for the DDE calibration problem. However, this empirical method does not benefit from any convergence guarantees.

We have implemented the naive StEFCal-based method for DDE calibration described in Section 6, considering the measurements made by $n_{a}=200$ antennas for a single time interval $T=1$, and considering DDEs Fourier kernels with support size $S=7\times 7$, with standard deviation to be $\upsilon=0.05$ (see Section 7.1). Moreover, we consider images with second level containing $10$ sources such that ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}}_{1})=0.1$.

As already shown in the previous experiments, our method recover the exact position of $100\%$ of the sources in $\overline{\boldsymbol{\epsilon}}$ and the average SNR of the recovered image (over 10 experiments) is equal to $10.67$ dB. In comparison, the proposed generalized StEFCal based method recovers the position of $98\%$ of the faint sources on average, with an average SNR equal to $8$ dB. Note also that our method is on average 4 times faster since, as explained above, it does not require the DDEs estimation and image estimation steps to converge completely.

The $\ell_{2}$ norm of the residual image \|\boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N} is equal to $9.75\times 10^{-4}$ (resp. $1.88\times 10^{-4}$) with $\boldsymbol{\mathsf{G}}^{\star}$ and $\boldsymbol{x}^{\star}$ obtained using our method (resp. using the proposed StEFCal based method). Moreover, the $\ell_{2}$ norm of the residual image \|\boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}\|_{2}/\sqrt{N} is equal to $1.31\times 10^{-3}$ (resp. $5.78\times 10^{-3}$) with $\boldsymbol{x}^{\star}$ obtained using our method (resp. using the proposed StEFCal based method). These results show that the $\ell_{2}$ norm of the residual image with the estimated DDEs is smaller using the StEFCal based method than using our method. However, it suggests that the true $\ell_{2}$ norm of the residual image is smaller using our method.

Therefore, according to the comments given above, we can conclude that, compared to the proposed naive StEFCal based method, our method based on the alternated forward-backward algorithm is faster, more stable and less sensitive to ambiguity problems.

7.3 Reconstruction of images with an extended source

In this second part, we consider an image of M31 of size $128\times 128$ given in Figure 11 (left). The objective of the simulations presented in this section is to show the behaviour of the proposed algorithm in the context of images involving extended sources. Therefore we consider only one simulation setting, consisting of measurements acquired by $n_{a}=100$ antennas at $T=10$ instants. The corresponding discrete $u-v$ coverage is shown in Figure 11 (right). Moreover, we fix the size of the support of the Fourier direction-dependent kernels to be $S=9\times 9$, with associated standard deviation, defined in Section 7.1, set to be $\upsilon=0.05$. As described in Section 5.1, we consider that the original image of M31 is decomposed as $\overline{\boldsymbol{x}}=\boldsymbol{x}_{o}+\overline{\boldsymbol{\epsilon}}$, where $\boldsymbol{x}_{o}$ is known, while $\overline{\boldsymbol{\epsilon}}$ is unknown and has to be estimated jointly with the DDEs. In this case, we construct $\boldsymbol{x}_{o}$ and $\overline{\boldsymbol{\epsilon}}$ such that

[TABLE]

where $\kappa>0$. In the performed simulations, we investigate the reconstruction of the image $\overline{\boldsymbol{\epsilon}}$ for the cases when $\kappa\in\{0.1,0.5,1\}$.

In this section, the original image is not sparse in its domain. Therefore, we propose to use the function given by eq. (34) as regularization term for the image, considering $\boldsymbol{\Psi}$ to be the SARA collection of wavelets (Carrillo et al., 2012). More precisely, $\boldsymbol{\Psi}$ is chosen to be the concatenation of a Dirac basis with the first eight Daubechies wavelets (Mallat, 2009). Moreover, we fix the regularization parameters in eq. (23) to be $\eta=10^{-2}$ and $\nu=2000$ for the initialization estimating the zero spatial frequency coefficients of DDEs, and $\nu=1000$ for the main algorithm estimating DDEs. Finally, the different parameters for Algorithm 1 are given in Table 2, including the stopping criteria defined in Section 5.6.

The results obtained considering $\kappa\in\{0.1,0.5,1\}$ are displayed in Figures 12, 13, and 14, respectively. For each of these figures, the first row shows in the left image the known approximation $\boldsymbol{x}_{o}$, and in the right image the unknown background $\overline{\boldsymbol{\epsilon}}$ to be estimated. The second row shows the estimated image $\boldsymbol{\epsilon}^{\star}$ of $\overline{\boldsymbol{\epsilon}}$, obtained (left) using the StEFCal-FB method solving only for DIEs, and (right) after the complete estimation of $\overline{\boldsymbol{\epsilon}}$ and the DDEs using the proposed method. Note that, as mentioned before, the StEFCal-FB methods gives similar reconstruction results as those obtained after the initialization of our method considering only the zero spatial frequency coefficients of the DDEs described in Section 5.5 (Step 1). The SNR between the true $\overline{\boldsymbol{\epsilon}}$ and its estimate $\boldsymbol{\epsilon}^{\star}$ obtained using StEFCal-FB is equal to $0.65$ dB (resp. $5.52$ dB and $3.64$ dB) for $\kappa=0.1$ (resp. $\kappa=0.5$ and $\kappa=1$). Similarly, the estimate $\boldsymbol{\epsilon}^{\star}$ obtained after the complete estimation of the DDEs using our method has an SNR equal to $16.83$ dB (resp. $17.47$ dB and $13.54$ dB) for $\kappa=0.1$ (resp. $\kappa=0.5$ and $\kappa=1$). These results suggest that our method is efficient not only for the reconstruction of point sources, but also to estimate the background of extended sources. Similarly to the simulations presented in Section 7.2.1, the accuracy of the reconstruction for extended sources depends on the total energy of the known approximation $\boldsymbol{x}_{o}$ with respect to the total energy of the unknown image $\overline{\boldsymbol{\epsilon}}$. As observed previously, if ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}})\ll{\mathsf{E}}(\boldsymbol{x}_{o})$, our method reconstructs accurately the unknown sources. However, in the case when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}})$ is of the same order as ${\mathsf{E}}(\boldsymbol{x}_{o})$, the reconstruction is more difficult. Nonetheless, in the worst case considered when ${\mathsf{E}}(\overline{\boldsymbol{\epsilon}})={\mathsf{E}}(\boldsymbol{x}_{o})$ (i.e. $\kappa=1$) presented in Figure 14, visually it can be observed that our method gives a good estimate $\epsilon^{\star}$ of $\overline{\boldsymbol{\epsilon}}$. The two last rows of Figures 12, 13 and 14 are dedicated to the residual images obtained considering $\kappa\in\{0.1,0.5,1\}$. More precisely, the third rows show the residual images |\boldsymbol{\mathsf{F}}^{\dagger}\boldsymbol{\mathsf{G}}^{\star\,\dagger}\big{(}\boldsymbol{\mathsf{G}}^{\star}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}| considering the estimated $\boldsymbol{\mathsf{G}}^{\star}$ and $\boldsymbol{x}^{\star}$ obtained (left) using the StEFCal-FB method solving only for DIEs, and (right) after the complete estimation of the image and the DDEs using our method. Similarly, the last rows show the residual images |\boldsymbol{\mathsf{F}}^{\dagger}\overline{\boldsymbol{\mathsf{G}}}^{\dagger}\big{(}\overline{\boldsymbol{\mathsf{G}}}\boldsymbol{\mathsf{F}}\boldsymbol{x}^{\star}-\boldsymbol{y}\big{)}| where $\overline{\boldsymbol{\mathsf{G}}}$ corresponds to the true unknown DDEs. As expected, it can be observed that estimating the full DDEs with respect to estimating only the DIEs leads to residual images with smaller amplitudes for all the cases presented. Moreover, for the final results with DDEs estimation, the residual images have smaller amplitude considering the estimated DDEs (third rows) than the original DDEs (fourth rows). This observation suggests that there may be ambiguity errors between the image and the direction-dependent Fourier kernels which need to be corrected.

8 Conclusions

We propose a non-convex optimization algorithm to jointly calibrate DDEs and estimate the sky intensity. The algorithm is based on a block-coordinate forward-backward approach, making use of suitable priors on both the image and the DDEs. More precisely, we use an $\ell_{1}$ based regularization term to promote sparsity on the image and the bright sources are assumed to be known, while we model the DDEs as smooth functions of the sky, i.e. spatially band-limited. Our method presents several advantages. Firstly, it benefits from convergence guarantees for both the image and the DDEs. Secondly, in contrast with DDE calibration methods developed recently, our method does not require a selection of calibrator directions since it constructs a smooth DDE screen applied to all the sources across the image. Finally, the proposed method is very general and can be easily adapted to the nature of the considered image. Indeed, we present two cases: the first case corresponds to images with point sources where the regularization term is chosen to promote sparsity of the image itself, while the second case deals with sophisticated extended sources where sparsity is promoted in a given dictionary. It is important to emphasize as well that, in the proposed algorithm, the DDEs are updated in parallel for all the antennas and at all time instants.

We study the performance of the proposed method considering a large class of simulations, varying the images, their dynamic range, the size and the variations of the direction-dependent Fourier kernels, the antennas distribution and taking into account the time dependency. In all the presented simulations, we show that our method leads to better reconstruction quality than obtained by only estimating the DIEs. Moreover, our simulations suggest that using our method to jointly estimate the DDEs and the image result in significant improvement of the the dynamic range, which is orders of magnitude higher when compared to accounting for DIEs only. Remark that all the parameters of the proposed method are fixed for all the different simulations presented in this paper. Therefore, it is worth noting that our method leads to very good reconstructions using a unique general framework converging in a complete automatic manner. The MATLAB code of the proposed method is available on GitHub (https://basp-group.github.io/joint-img-dde-calibration/).

Even though the proposed method is very promising, we intend to improve some points in future works. Firstly, in the presented simulations the bright sources in the images are assumed to be known exactly. Thus, we intend to extend our method to be able to take into account possible errors in the amplitudes and/or the positions of the bright sources. Secondly, we plan to investigate further regularization terms in order to reduce the ambiguity error appearing between the DDEs and the image. Moreover we aim to perform simulations on realistic data simulations, incorporating the gridding steps in the algorithm to use the non-uniform FFT, considering unknown size support for the direction-dependent Fourier kernels, and modelling time dependence in order to investigate the robustness of the proposed method. Finally, we plan to compare our approach to other state of the art algorithms.

Acknowledgements

This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC, grant EP/M008843/1). We would like to thank Oleg Smirnov for insightful discussions.

Appendix A Tables of notations

In this paper, we use italic bold letters to denote vectors, italic letters to denote scalar variables, and capital bold letters for matrices. For a given variable, e.g. $\boldsymbol{x}$, $x(n)$ denotes its $n$-th coefficient, $\widehat{\boldsymbol{x}}$ denotes its Fourier transform, $\overline{\boldsymbol{x}}$ denotes the original unknown associated variable, and $\boldsymbol{x}^{\star}$ denotes the estimated associated variable.

Table 3 gives the notations used for the main variables describing the radio interferometric inverse problem, Table 4 gives the notations used for the imaging and calibration minimization problem, and Table 5 gives the notations used for the proposed algorithm.