Super-resolution Full Polarimetric Imaging for Radio Interferometry with Sparse Modeling

TL;DR

This paper introduces a sparse modeling technique for radio interferometry that achieves super-resolution full polarization imaging, surpassing traditional methods in resolution and fidelity, demonstrated through simulations of M87 with the EHT.

Contribution

The paper presents a novel sparse modeling approach using combined $ extit{l}_1$ and TV regularization for super-resolution full-polarimetric imaging in radio interferometry.

Findings

Achieves 25-30% of diffraction limit resolution.

Outperforms CLEAN and single regularization methods.

Provides higher image fidelity across multiple scales.

Abstract

We propose a new technique for radio interferometry to obtain super-resolution full polarization images in all four Stokes parameters using sparse modeling. The proposed technique reconstructs the image in each Stokes parameter from the corresponding full-complex Stokes visibilities by utilizing two regularization functions: the $\ell _1$-norm and total variation (TV) of the brightness distribution. As an application of this technique, we present simulated linear polarization observations of two physically motivated models of M87 with the Event Horizon Telescope (EHT). We confirm that $\ell _1$+TV regularization can achieve an optimal resolution of $\sim 25-30$\% of the diffraction limit $\lambda/D_{\rm max}$, which is the nominal spatial resolution of a radio interferometer for both the total intensity (i.e. Stokes $I$) and linear polarizations (i.e. Stokes $Q$ and $U$). This optimal resolution is better than that obtained from the widely used Cotton-Schwab CLEAN algorithm or from using $\ell _1$ or TV regularizations alone. Furthermore, we find that $\ell _1$+TV regularization can achieve much better image fidelity in linear polarization than other techniques over a wide range of spatial scales, not only in the super-resolution regime, but also on scales larger than the diffraction limit. Our results clearly demonstrate that sparse reconstruction is a useful choice for high-fidelity full-polarimetric interferometric imaging.

Equations

$$\underset{\mathbf{x}}{\min }||\mathbf{y}-\mathbf{Ax}|{|}_2^2\mathrm{subject}\mathrm{to}||\mathbf{x}|{|}_0\le T,$$

$$||\mathbf{x}|{|}_p={\left(\sum _i|x_i{|}^p\right)}^{\frac{1}{p}}$$

$$\underset{\mathbf{x}}{\min }||\mathbf{y}-\mathbf{Ax}|{|}_2^2\mathrm{subject}\mathrm{to}||\mathbf{x}|{|}_1\le T,$$

$$\underset{\mathbf{x}}{\min }\left(||\mathbf{y}-\mathbf{A}\mathbf{x}|{|}_2^2+{\Lambda }_{\ell }||\mathbf{x}|{|}_1\right).$$

$$\tilde{S}(u,v)=\int dxdyS(x,y)\exp (-i2\pi (ux+vy)),$$

$$\tilde{\mathbf{S}}=\mathbf{F}\mathbf{S}(\mathrm{for}\mathbf{S}=\mathbf{I},\mathbf{Q},\mathbf{U},\mathbf{V}).$$

$$\mathbf{S}={\mathrm{argmin}}_{\mathbf{S}}\left(||\tilde{\mathbf{S}}-\mathbf{F}\mathbf{S}|{|}_2^2+{\Lambda }_{\ell }||\mathbf{S}|{|}_1+{\Lambda }_t||\mathbf{S}|{|}_{\mathrm{tv}}\right)$$

$${\tilde{\Lambda }}_{\ell }\equiv {\Lambda }_{\ell }\underset{i}{\max }|\widetilde{I_i}|/N,$$

$$||\mathbf{x}|{|}_{\mathrm{tv}}=\sum _i\sum _j\sqrt{|x_{i+1,j}-x_{i,j}{|}^2+|x_{i,j+1}-x_{i,j}{|}^2}.$$

$${\tilde{\Lambda }}_t\equiv 4{\Lambda }_t\underset{i}{\max }|\widetilde{I_i}|/N,$$

$$\mathrm{NRMSE}(\mathbf{I},\mathbf{K})=\sqrt{\frac{{\sum }_i|I_i-K_i{|}^2}{{\sum }_i|K_i{|}^2}},$$

$$||\mathbf{x}|{|}_{\mathrm{t}{\mathrm{v}}^2}=\sum _i\sum _j\left(|x_{i+1,j}-x_{i,j}{|}^2+|x_{i,j+1}-x_{i,j}{|}^2\right).$$

$$\underset{\mathbf{x}\in C}{\min }\{F(\mathbf{x})\equiv f(\mathbf{x})+g(\mathbf{x})\},$$

$$Q(\bm{x},\bm{y};L)=f(\bm{y})+\bigl{\langle}\bm{x}-\bm{y},\nabla f(\bm{y})\bigr{\rangle}+\frac{L}{2}\bigl{\|}\bm{x}-\bm{y}\bigr{\|}_{2}^{2}+g(\bm{x}),$$

$$p_{C}(\bm{y};L)=\operatorname*{argmin}_{\bm{x}\in C}Q(\bm{x},\bm{y};L)=\operatorname*{argmin}_{\bm{x}\in C}\Biggl{\{}\frac{L}{2}\Bigl{\|}\bm{x}-\bigl{(}\bm{y}-\frac{1}{L}\nabla f(\bm{y})\bigr{)}\Bigr{\|}_{2}^{2}+g(\bm{x})\Biggr{\}}.$$

$$F({\mathbf{x}}_k)-F({\mathbf{x}}^{\ast })\le \frac{2\alpha L(f)\Vert {\mathbf{x}}_0-{\mathbf{x}}^{\ast }{\Vert }_2^2}{(k+1{)}^2},\forall {\mathbf{x}}^{\ast }\in C.$$

$$f(\mathbf{I})=\Vert \tilde{\mathbf{I}}-\mathbf{FI}{\Vert }_2^2+{\Lambda }_{\ell }\sum _i{I}_i,g(\mathbf{I})={\Lambda }_t\Vert \mathbf{I}{\Vert }_{\mathrm{tv}},C=\{\mathbf{I}|I_i\ge 0,\text{mbox}for\forall i\}.$$

$$f(\mathbf{S})=\Vert \tilde{\mathbf{S}}-\mathbf{FS}{\Vert }_2^2,g(\mathbf{S})={\Lambda }_1\Vert \mathbf{S}{\Vert }_1+{\Lambda }_t\Vert \mathbf{S}{\Vert }_{\mathrm{tv}},\mathbf{S}=\mathbf{Q},\mathbf{U}\text{mbox}or\mathbf{V}.$$

Full text

techniques: high angular resolution — techniques: image processing — techniques: interferometric — techniques: polarimetric — polarization

††slugcomment: Accepted for publication in the Astronomical Journal

1 Introduction

Polarization is a unique tracer of the magnetic field and magnetized plasma distribution in the universe. The distribution of magnetic field lines can be inferred from linear polarization in a variety of sources, including synchrotron emission from non-thermal or relativistic thermal electrons in high-energy objects (e.g. Pacholczyk, 1970), maser emission from star-forming regions or evolved stars (e.g. Fish & Reid, 2006), and thermal emission partially absorbed by aligned dust grains (e.g. Girart et al., 2006, 2009). Polarized emission also contains information about the magnetized plasma along the line of sight (e.g. Faraday Tomography; Burn, 1966; Brentjens & de Bruyn, 2005) via Faraday rotation of linear polarization or Faraday conversion from linear to circular polarization (Legg & Westfold, 1968; Pacholczyk, 1970; Jones & Odell, 1977). Recent theoretical and observational studies suggest that linear polarization can be a unique tracer of the dust evolution in proto-planetary disks (Kataoka et al., 2015, 2016).

High-resolution imaging of polarized emission is therefore a fundamental part of the modern observational toolkit in astronomy. The angular resolution of a telescope (often referred to as “beam size” in radio astronomy and “diffraction limit” in optical astronomy) is given by $\theta\approx\lambda/D$, where $\lambda$ and $D$ are the observing wavelength and the diameter of the telescope, respectively. At radio wavelengths, interferometry is the most effective approach to obtain high angular resolution. The nominal resolution of an interferometer is given by $\theta\approx\lambda/D_{\rm max}$, where $D_{\rm max}$ is the maximum length of the baseline between two telescopes, projected in the plane normal to the direction of observation. Of all observational techniques across the electromagnetic spectrum, radio interferometry provides the imaging capability at the finest angular resolution (e.g. Thompson et al., 2001). In particular, very long baseline interferometry (VLBI), which utilizes intercontinental baselines (or even baselines to space), has achieved the highest angular resolution in the history of astronomy.

The Event Horizon Telescope (EHT; Doeleman et al., 2009a) is a ground-based VLBI array with an angular resolution of a few tens of microarcseconds at short/sub-millimeter wavelengths ($\lambda\lesssim 1.3$ mm, $\nu\gtrsim 230$ GHz) (e.g. Doeleman et al., 2008, 2012; Fish et al., 2011, 2016; Lu et al., 2012, 2013; Akiyama et al., 2015; Johnson et al., 2015). The EHT resolves compact structures of a few Schwarzschild radii ($R_{s}$) in the vicinity of the supermassive black holes in the Galactic Center source Sgr A* (Doeleman et al., 2008; Fish et al., 2011, 2016; Johnson et al., 2015) and the nucleus of M87 (Doeleman et al., 2012; Akiyama et al., 2015). Direct imaging of these scales will be accessible in the next few years with technical developments and the addition of new (sub)millimeter telescopes such as the Atacama Large Submillimeter/millimeter Array (ALMA) to the EHT (e.g. Fish et al., 2013). Polarimetric imaging with the EHT will be especially transformational, opening a new field to study the properties of the magenetic field distribution and magnetized plasma in the regions of strong gravitation. Indeed, early EHT observations of Sgr A* resolves ordered and time-variable magnetic fields on $R_{s}$ scales (Johnson et al., 2015). High-fidelity images of the linearly polarized emission will be crucial for understanding processes of black hole accretion and jet formation.

The imaging problem of interferometry is formulated as an underdetermined linear problem (see §2.1) of reconstructing an image from complex visibilities that represent Fourier components of the source image. The CLEAN algorithm (Högbom, 1974) and its variants (e.g. Clark, 1980; Schwab, 1984) have been the most successful and widely used algorithms in radio interferometry. CLEAN was independently rediscovered as the Matching Pursuit algorithm (MP; Mallat & Zhang, 1993) and has been widely used in many other fields to derive a sparse solution ${\bf x}$ of an underdetermined linear problem ${\bf y}={\bf Ax}$, where ${\bf y}$ and ${\bf A}$ are observational data and observation matrix, respectively. For real data with noise, this can be mathematically described by

[TABLE]

where $||{\bf x}||_{p}$ is the $\ell_{p}$-norm of the vector ${\bf x}$ given by

[TABLE]

for $p>0$, and defined as the number of non-zero components for $p=0$. The term to be minimized is the traditional $\chi^{2}$ term, and $T$ is a threshold for the $\ell_{0}$-norm representing the sparsity of the solution. Thus, the solution is equivalent to minimizing the $\chi^{2}$ term within a given sparsity. A direct approach to solve this equation is to try all possible combinations of zero components of ${\bf x}$ one-by-one. However, the computational cost of this exhaustive search is so large that it is intractable for large dimensional $\bm{x}$. CLEAN, MP and their variants are select a non-zero element one-by-one and incrementally in a greedy manner in order to solve this problem efficiently.

A popular relaxation of sparse reconstruction comes from replacing the $\ell_{0}$-norm with the $\ell_{1}$-norm as

[TABLE]

which is known as LASSO (least absolute shrinkage and selection operator; Tibshirani, 1996). This is a convex relaxation of Eq. (1), and can be transformed in the Lagrange form,

[TABLE]

Many efficient algorithms have been proposed to solve LASSO (e.g., the fast iterative shrinkage-thresholding algorithm (FISTA); Beck & Teboulle, 2009b). This method has become popular in many fields such as medical imaging, particularly after the appearance of compressed sensing (also known as compressive sensing; Donoho, 2006; Candes & Tao, 2006) techniques, which have shown that LASSO can solve many ill-posed linear problems accurately if the solution vector is sparse — the number of elements with nonzero value is much small compared to its dimension. We have applied LASSO to Stokes $I$ imaging with radio interferometry (Honma et al., 2014; Ikeda et al., 2016; Akiyama et al., 2017) and found that LASSO can potentially reconstruct structure on scales $\sim 4$ times finer than $\lambda/D_{\rm max}$ (Honma et al., 2014). Techniques of compressed sensing are beginning to be used in other fields of radio interferometry (see Garsden et al., 2015, and references therein).

A critical assumption in techniques with $\ell_{1}$ regularization is that the solution (i.e., the true image) is sparse. If the number of pixels with nonzero brightness is not small compared to the number of data points, simple $\ell_{1}$-norm regularization may reconstruct an image that is too sparse. This situation can arise when reconstructing an extended source or even for a compact source if the imaging pixel size is set to be much smaller than the size of the emission structure. A promising approach to overcome this issue is to change the basis of the image to a more sparse one. Pioneering work in this area has made use of transforms to wavelet or curvelet bases, in which the image can be represented sparsely (e.g. Li et al., 2011; Carrillo et al., 2014; Garsden et al., 2015; Dabbech et al., 2015). We have taken another approach by adding total variation (TV) regularization (e.g. Wiaux et al., 2010; McEwen & Wiaux, 2011; Uemura et al., 2015; Chael et al., 2016), which produces an image that is sparse in its gradient domain. TV regularization has been shown to be effective for imaging with visibility amplitudes and closure phases (e.g. Akiyama et al., 2017) in the super-resolution regime finer than the diffraction limit.

In interferometric imaging, another class of widely-used imaging techniques is the Maximum Entropy Methods (MEM), utilizing different functions (named as “entropy terms”) to regularize images (see Narayan & Nityananda, 1986, for a review). Image reconstruction with MEM has been practically extended to polarimetry (Holdaway & Wardle 1990; Sault et al. 1999; and see Chael et al. 2016 and Coughlan & Gabuzda 2016 for a review of polarimetric MEM techniques).

Sparse modeling techniques utilizing $\ell_{1}$ and TV terms have heretofore been applied only to Stokes $I$ image reconstruction. In this paper, we extend the framework of sparse imaging techniques for radio interferometry with $\ell_{1}$ and TV regularizations to full-polarization imaging for the first time. As an example, we apply our new technique to simulated EHT data of the accretion and jet launching region immediately around the black hole of M87.

2 The proposed method

2.1 Polarimetric imaging with radio intererometry

The intensity distribution of the emission from the sky can be described with four Stokes parameters, $I$, $Q$, $U$ and $V$, which are all real. Stokes $I$ represents the total intensity of the emission, which is generally non-negative for astronomical images. On the other hand, $Q$ and $U$, which represent linear polarization, and $V$, which represents circular polarization, can take on negative values. Stokes $Q$ and $U$ are often combined into the complex quantity $P\equiv Q+iU$, where $|P|$ and $\chi=\arg(P)/2$ are the linear polarization intensity and the electric vector polarization angle (EVPA), respectively.

A radio interferometer samples Fourier components of each Stokes parameter, known as the Stokes visibilities $\tilde{I}$, $\tilde{Q}$, $\tilde{U}$ and $\tilde{V}$ defined by

[TABLE]

where $S$ and $\tilde{S}$ represent a Stokes parameter and corresponding Stokes visibility (i.e. $S=I,Q,U,V$). Here, the spatial frequency $(u,v)$ corresponds to the baseline vector (in units of the observing wavelength $\lambda$) between two antennas (or receivers) projected to the tangent plane of the celestial sphere at the phase-tracking center.

Observed visibilities are discrete quantities, and the sky image can be approximated by a pixellated version where the pixel size is much smaller than the nominal resolution of the interferometer. The Stokes parameters can therefore be represented as discrete vectors ${\bf I}$, ${\bf Q}$, ${\bf U}$ and ${\bf V}$, related to the observed Stokes visibilities $\tilde{{\bf I}}$, $\tilde{{\bf Q}}$, $\tilde{{\bf U}}$ and $\tilde{{\bf V}}$ by a discrete Fourier transform ${\bf F}$:

[TABLE]

The sampling of Stokes visibilities is almost always incomplete. Since the number of visibility samples $\tilde{{\bf S}}$ is smaller than the number of pixels in the image, solving the above equation for the image ${\bf S}$ is an ill-posed problem. One or more regularizers must therefore be added to find a unique solution to equation (6).

2.2 The Proposed Methods

A natural extension of our previous work (Honma et al., 2014; Ikeda et al., 2016; Akiyama et al., 2017) to full polarimetric imaging is given by

[TABLE]

for each Stokes parameter (i.e. ${\bf S}={\bf I},~{}{\bf Q},~{}{\bf U},~{}{\bf V}$) and corresponding Stokes visibility (i.e. $\tilde{{\bf S}}=\tilde{{\bf I}},~{}\tilde{{\bf Q}},~{}\tilde{{\bf U}},~{}\tilde{{\bf V}})$. This equation consists of the traditional $\chi^{2}$ term, which represents deviations between the model image and observed visibilities, and two terms consisting of a regularizer and a regularization parameter.

One of the additional terms is $\ell_{1}$-regularization (e.g. Honma et al., 2014). $\Lambda_{\ell}$ is its regularization parameter, adjusting the degree of sparsity by changing the weight of the $\ell_{1}$-norm penalty. In general, a large $\Lambda_{\ell}$ prefers a solution with very few nonzero components, while a small $\Lambda_{\ell}$ imposes less sparsity. In this paper, we use the normalized regularization parameter $\tilde{\Lambda}_{\ell}$ defined by

[TABLE]

which is less affected by the number of visibilities $N$ and the total flux density of the target source that should be close to the maximum value of the visibility amplitudes at Stokes $I$ (i.e. $\max_{i}|\tilde{I}_{i}|$), following Akiyama et al. (2017).

The other additional term is total variation (TV) regularization with a regularization parameter $\Lambda_{t}$. A large $\Lambda_{t}$ leads to a sparse solution in the gradient domain – a piecewise smooth solution. In this paper, we adopt the isotropic TV expression (Rudin et al., 1992), a typical form for two-dimensional images, defined by

[TABLE]

We have examined the effects of TV regularization on Stokes $I$ imaging in our previous work (Ikeda et al., 2016; Akiyama et al., 2017), and TV regularization is also used in other applications, such as Doppler tomography (Uemura et al., 2015). As with the $\ell_{1}$-norm, we use a normalized regularization parameter $\tilde{\Lambda}_{t}$ defined by

[TABLE]

again following Akiyama et al. (2017).

The Stokes $I$ image is solved with a non-negative condition (i.e. ${\bf I}\geq 0$). The linear polarization image (henceforth $P$ image) is calculated from reconstructed $Q$ and $U$ images by $P=Q+iU$. In other words, Stokes $Q$ and $U$ images are solved independently. Since the Stokes $Q$ and $U$ images can be negative, we solve these images without the non-negative condition.

The optimization problem, equation (7), is convex, and therefore its solution is uniquely determined regardless of initial conditions. Many algorithms have been proposed to solve this problem. We adopt the fast iterative shrinking thresholding algorithm (FISTA), originally proposed by Beck & Teboulle (2009b) for $\ell_{1}$ regularization and by Beck & Teboulle (2009a) for TV regularization. We use a monotonic FISTA algorithm (MFISTA) designed for $\ell_{1}$+TV regularization (see Appendix A for details).

3 Imaging Simulations

3.1 Physically Motivated Models

In this paper, we adopt two physically motivated models of the 1.3 mm emission from M87 on event-horizon scales (Figure 1). In this paper, we focus on imaging the total intensity $I$ and linear polarization $Q$ and $U$ emission

The first model is a simple force-free jet model (hereafter, forward-jet model) in the magnetically dominated regime Broderick & Loeb (2009); Lu et al. (2014). We adopt the model image presented in Lu et al. (2014), which is based on model parameters fitted to the spectral energy distribution of M87 and results of EHT observations at 1.3 mm (Doeleman et al., 2012). The approaching jet is the dominant feature in this model.

The second model (henceforth, counter-jet model) is based on results of GRMHD simulations (Dexter et al., 2012) and full polarimetric general relativistic radiation transfer calculations (Dexter 2016; Mościbrodzka, Dexter, Davelaar et al., in prep.). The dominant emission region is the counter jet illuminating the last photon orbit.

3.2 Simulated Observations

We simulate observations of model M87 images with the EHT at 1.3 mm (230 GHz) using the MIT Array Performance Simulator (MAPS)111http://www.haystack.mit.edu/ast/arrays/maps/. In most aspects, the data generation parameters are identical to those used in Akiyama et al. (2017), except that here we use an integration time of 10 sec. We simulate data for a six-station array with a band width of 3.5 GHz at each polarization, system temperatures described in the proposer’s guid of 1-mm VLBI observations in ALMA Cycle 4, and a correlation efficiency of 0.7 that includes a quantization efficiency of 2-bit sampling and other potential losses. Observations are performed with an observational efficiency of 25% in time, during a GST (Greenwich sidereal time) range of 13-0 hour. This GST range corresponds to a time range when M87 can be observed by either of two anchor stations of the EHT, the Atacama Large Millimeter/submillimeter Array (ALMA) or the Large Millimeter Telescope (LMT), at an elevation greater than 20∘. Here, we consider only thermal errors. See Akiyama et al. (2017) for more details about the conditions of simulated observations. Fig. 2 shows the resultant $uv$-coverage of simulated observations. Note that the maximum baseline length of observations is 7.2 G$\lambda$, corresponding to $\lambda/D_{\rm max}=28.5$ $\rm\mu$as.

3.3 Imaging

We reconstruct Stokes $I$, $Q$ and $U$ images from simulated data-sets based on the method described in §2.2. In addition, we attempt to solve images with the Cotton-Schwab CLEAN algorithm (henceforth, CS-CLEAN; Schwab, 1984) using uniform weighting to evaluate the relative performance of our techniques in the exactly same way with Akiyama et al. (2017). We use an implementation of CS-CLEAN in the Common Astronomy Software Applications (CASA) package222https://casa.nrao.edu/. We adopt a field of view (FOV) of 200 $\mu$as, gridded into 100 pixels in each of right ascension and declination for both models. The resulting pixel size of $\sim 2$ ${\rm\mu}$as corresponds to a physical scale of $\sim 0.21$ $R_{s}$.

The proposed method has two regularization parameters $\tilde{\Lambda}_{\ell}$ and $\tilde{\Lambda}_{t}$. Images at each Stokes parameter were reconstructed at 5 regularization parameters for both $\tilde{\Lambda}_{\ell}$ and $\tilde{\Lambda}_{t}$, equally spaced in logarithmic steps in the range $10^{-1},...,10^{+3}$. In addition to employing $\ell_{1}$+TV regularization, we also explore pure $\ell_{1}$ and pure TV regularization separately (i.e., $\tilde{\Lambda}_{t}=0$ or $\tilde{\Lambda}_{\ell}=0$, respectively).

We evaluate the goodness-of-fit for each image and then select the best-fit images with 10-fold Cross Validation (henceforth CV; Akiyama et al., 2017). The observational data (i.e. sampled visibilities) are randomly partitioned into 10 equal-sized subsamples. 9 of 10 subsamples are used in the image reconstruction as the training set, and we obtain the trained image. The remaining single subsample is used as the validation set for testing the model using $\chi^{2}$. The $\chi^{2}$ between the validation set and the image from the training set, which is so-called the validation error, is a good indicator of goodness-of-fit. For too small regularization parameters, the trained image would be over-fitted and too complicated, resulting in a large deviation between the trained image and the validation set (i.e. large validation error). On the other hand, for too large regularization parameters, the trained image would be too simple and not well-fitted to the training set, also resulting in a large validation error. Thus, reasonable parameters can be estimated by finding a parameter set that minimizes the validation error. We repeat the procedure by changing the subsample for validation data 10 times, until all subsamples are used for both training and validation. As a result, we get 10 validation errors. The validation errors are averaged and then used to determine optimal regularization parameters at each Stokes parameter that minimize the averaged validation error. Note that the final image is reconstructed by full sample of the observed visibilities.

To reduce the computational cost, we grid the observed visibilities with the classic cell-averaging method (see Thompson et al., 2001) prior to imaging. We adopt a FOV size of 2 mas for gridding, corresponding to a grid size of $\sim 0.1$ G$\lambda$ in $uv$-space.

3.4 Evaluation of the image fidelity

We evaluate the quality of reconstructed images in two ways. First, we employ the normalized root mean square error (NRMSE) metric following Chael et al. (2016) and Akiyama et al. (2017), defined as,

[TABLE]

where ${\bf I}$ and ${\bf K}$ are the image to be evaluated and the reference image, respectively. For linear polarization images, we use the complex linear-polarization intensity (i.e. ${\bf P}={\bf Q}+i{\bf U}$) to evaluate NRMSEs. Since both model images have finer resolutions than is recoverable using the EHT, we adjust the pixel size of the reconstructed image to that of the model image with bi-cubic spline interpolation. Second, we measure structural dissimilarity (Wang et al., 2004) between the model and reconstructed images using the DSSIM metric adopted in previous work (Lu et al., 2014; Fish et al., 2014). Since both metrics show a similar trend, we show only the behavior of the NRMSE in the figures that follow.

Of potential interest for future EHT observations is to detect hypothesized blob-like flaring structure(s) in the accretion flow or jet due to partially heated or overdense plasma (e.g. Broderick & Loeb, 2006; Doeleman et al., 2009b). However, image reconstructions can generate artifacts that mimic bright components, making it difficult to identify such signatures accurately. A useful evaluation tool for imaging techniques is to identify how many bright blobs appear in the reconstructed image. The input model images do not contain flaring structures, so reconstructed images that show more than one cluster of pixels falsely recover blob-like features. We therefore also perform a cluster analysis on each image using Density-Based Spatial Clustering of Applications with Noise (DBSCAN; Ester et al., 1996) to identify these false features. We configure DBSCAN as follows. The images contain a range of pixel brightness values; therefore, we cluster the pixels that have a brightness $>50\%$ of its peak intensity, of which separations are larger than 2 pixels $\rm\mu$as. Then we cluster the reconstructed images with the same parameters to find if false blob-like features (clusters) are obtained.

4 Results

4.1 Stokes $I$ images

The results for Stokes $I$ images of $\ell_{1}$, TV and $\ell_{1}$+TV regularizations are shown along with the model and CS-CLEAN images in Figure 3. We also plot the NRMSE metric for reconstructed images as well in the bottom panel, along with the full width at half maximum (FWHM) size of a convolving circular Gaussian beam. The best-case scenario — the differences from the original input due solely to a loss of resolution, not to errors in reconstructing the image — is shown by the black curve labeled “Model” following previous work (Chael et al., 2016; Akiyama et al., 2017). This is calculated by taking the NRMSE between the model image convolved with a circular Gaussian beam with a FWHM and the original (unconvolved) model image. The NRMSE of each of the reconstructed images, convolved with circular Gaussian beams, is shown in the bottom panel.

All techniques reconstruct Stokes $I$ images equally well on scales comparable to or greater than the diffraction limit. The NRMSEs of the reconstructed images only start to deviate from the model images in the super-resolution regime — namely on scales finer than the diffraction limit. In this regime, the NRMSEs are different by techniques. $\ell_{1}$-regularization and CS-CLEAN have a common trend for both models. The minimum errors are achieved at a resolution of $\sim 30-50$% of the diffraction limit, and then the NRMSEs show a rapid increase in errors at smaller scales, broadly consistent with previous studies on different model images and data products (Chael et al., 2016; Akiyama et al., 2017). This clearly shows that, on such small scales, the image is no longer sparse and breaks the underlying assumption of both techniques. In contrast, TV and $\ell_{1}$+TV regularizations show much more modest variations in the super-resolution regime. The minimum errors are achieved at a resolution of $\sim 25-30$% of the diffraction limit, smaller than $\ell_{1}$-regularization and CS-CLEAN. In addition, the NRMSEs show only a slight increase in smaller scales. Both the TV and $\ell_{1}$+TV reconstructions produce images that have a smooth distribution similar to the model images, resulting in smaller errors than $\ell_{1}$-regularization and CS-CLEAN, even if the TV and $\ell_{1}$+TV are not convolved with a restoring beam.

A clustering analysis with DBSCAN shows that images with smoother regularizations ($\ell_{1}$+TV and TV) have only one cluster of bright imaging pixels regardless of resolution. The other two sparse techniques (CS-CLEAN and $\ell_{1}$) show more than one cluster in smaller resolutions, as clearly seen in Figure 3 (a) for both models. Thus, even though all techniques have similar optimal resolutions and minimum NRMSEs for the counter-jet model, the bright emission has more than one clusters at optimal resolutions for CS-CLEAN and $\ell_{1}$. This indicates that sparse reconstructions with $\ell_{1}$-regularization and CS-CLEAN are more likely to misidentify flaring substructures. We also note that, simultaneously, this indicates that the NRMSE and DSSIM image fidelity metrics may not always be an appropriate indicator for goodness of feature reconstruction.

4.2 Linear polarization images

The results for linear polarization images (Stokes $Q$ and $U$) are shown in Figure 4. Similar to Stokes $I$ (Figure 3), we show the model/reconstructed images at the optimal resolution of $\ell_{1}$+TV regularization in panel (a) and the NRMSEs in panel (b). The NRMSEs show different trends for polarization than they do for Stokes $I$ §4.1 because Stokes $Q$ and $U$ can be negative and because the linearly polarized flux is significantly smaller than in total intensity, leading to a smaller signal-to-noise ratio in these simulations.

CS-CLEAN shows the highest NRMSE over almost the entire range of spatial scales, and its NRMSEs do not have a global minimum at resolutions smaller than 1.5$\lambda/D_{\rm max}$. Indeed, as shown in Figure 4, CS-CLEAN can recover only a tiny fraction of linearly polarized emission, and the reconstructed EVPA distribution is inaccurate.

TV regularization, which shows good performances for Stokes $I$ imaging, is insufficient on its own for polarimetric imaging. The much lower signal-to-noise ratios in Stokes $Q$ and $U$ visibilities require a regularization parameter $\sim$10 times larger than for Stokes $I$ to minimize the validating error of 10-fold CV, resulting in image distributions that are much blurrier than the model images. The TV-regularized $P$ images reconstruct bright emissions better than CS-CLEAN, but there are a lot of artificial diffuse emissions that dominate the NRMSEs, as seen by the flat curves in Figure 4 (b).

The $\ell_{1}$ regularization exhibits better NRMSEs than CS-CLEAN for both models on most spatial scales. Since $\ell_{1}$ regularization suppresses the artificial diffuse emission seen in TV regularization, the achieved NRMSE is better than TV until at $\gtrsim 20$% of the diffraction limit. However, as with Stokes $I$, the images become too sparse on scales smaller than $\sim 30-50$% of the diffraction limit, causing a rapid rising in NRMSEs.

$\ell_{1}$+TV regularization provides reasonable linear polarization images with the most reasonable sparseness and smoothness, stably showing good performance across the entire range of spatial scales. The optimal resolution of $\sim 20-25$% is comparable to Stokes $I$ and is the best among the four techniques. These results clearly shows that $\ell_{1}$+TV can achieve the best image fidelity among techniques presented in this paper not only for super-resolution imaging but also more general imaging on scales larger than diffraction limit.

5 Discussions and Summary

We have presented a new technique for full polarimetric imaging with radio interferometry using sparse modeling. As shown in §3 and §4, $\ell_{1}$+TV regularization stably shows better performance than either $\ell_{1}$ or TV regularization alone, and than the most widely-used Cotton-Schwab CLEAN. This applies regardless of Stokes parameters. In particular, the superiority of the combined $\ell_{1}$+TV regularization is significant for linear polarizations on a wide range of spatial scales from super-resolution regimes to scales larger than the diffraction limit.

Our technique can be applied to most existing radio interferometers whose data products are full complex visibilities in all four Stokes parameters. Although we did not image the circular polarization (Stokes $V$) in this work, our results suggest that $\ell_{1}$+TV regularization would likely achieve a better performance than the Cotton-Schwab CLEAN for circular polarization too, since it is mathematically similar to Stokes $Q$ and $U$ imaging.

We note future prospects for the application of our techniques to VLBI, including future EHT observations. In VLBI observations, the absolute phase information generally can not be obtained due to non-synchronized local oscillators and quite different atmospheric phase delays at different sites (see Thompson et al., 2001). In addition, at short-mm/sub-mm wavelengths, even the source visibility phase cannot be measured due to the rapidly varying atmospheric delays. In VLBI, the visibility phase is traditionally calibrated based on phase closure, using the self-calibration technique with hybrid/differential mapping (e.g. Walker, 1995). The systematic phase errors derived using Stokes $I$ data can be applied to Stokes $Q$ and $U$, since station-based systematic errors should be the same among the Stokes parameters. Our techniques can be applied to the VLBI data after self-calibrating Stokes $I$ data with the traditional hybrid/differential mapping or the Stokes $I$ image obtained with new state-of-art imaging techniques based on closure quantities such as Maximum Entropy Methods (e.g. Buscher, 1994; Chael et al., 2016), a patch prior (CHIRP; Bouman et al., 2015), sparse modeling (Akiyama et al., 2017) and PRECL (Ikeda et al., 2016), which have been developed for optical interferometers and/or the EHT. In a forthcoming paper, we will evaluate the performance of the proposed technique for data with station-based systematic phase errors, which are common in VLBI.

We also note that there is a new method for Stokes I and linear polarization, very recently proposed in Chael et al. (2016), which is designed for VLBI. This method solves first for the Stokes-$I$ images from visibility amplitudes and closure phases at Stokes $I$. Then, the linear polarization images are solved using the reconstructed Stokes-$I$ images and complex polarimetric ratios (ratios of the Stokes $Q$, $U$ visibilities to the Stokes $I$ visibilities). All of these VLBI observables are robust against station-based phase errors. Through the above processes, Stokes $I$ visibility phases are recovered from imaging with visibility amplitudes and closure phases, and the visibility phases at Stokes $Q$ and $U$ are phase-referenced from Stokes $I$ through the polarimetric ratio. The Stokes-$I$ and polarization intensity (i.e. $|P|$) images are regularized by the entropy term of MEM, while the EVPA distribution (i.e. $\arg(P)/2$) is regularized independently by a smooth regularization term such as TV. Chael et al. (2016) demonstrate that this method can also achieve a better fidelity and superior optimal resolution than the Cotton-Schwab CLEAN. An advantage of this technique is that it can simultaneously reconstruct $I$, $Q$ and $U$ images from robust VLBI observables. In addition, the reconstructed images strictly satisfy $|I|>|P|$, which can suppress artifacts in $P$ images in regions where $I$ is not bright. The disadvantage of this technique is that the optimization problem is highly non-linear and non-convex, and that the solution can potentially be initial-condition dependent not only in Stokes $I$ but also in Stokes $Q$ and $U$. Furthermore, $Q$ and $U$ images are reconstructed from the polarimetric ratio that can have larger uncertainties than the Stokes visibilities, particularly at long baselines, limiting the dynamic range, image sensitivity and optimal spatial resolution. An alternative, mathematically equivalent way — phase-referencing with self-calibration techniques — will avoid such disadvantages in polarimetric imaging.

Future work for the techniques proposed in this paper will include other sparse regularizations for multi-resolution imaging, such as $\ell_{1}$+wavelet/curvelet transformation (e.g. Li et al., 2011; Carrillo et al., 2014; Garsden et al., 2015; Dabbech et al., 2015). In addition, the application of, and experimentation with other forms of TV would be important. We have been using the most widely-used isotropic TV (Rudin et al., 1992) for TV regularization, which preserves sharp edges in the image. This would be useful for optically thick objects like stars, but might not be optimal for emission from optically thin objects that have smoother edges in general. An alternative form of sparse regularization in the gradient domain that favors smoother edges is given by, for instance,

[TABLE]

This function is also convex, similar to the TV term adopted in this work, and can be an alternative option. The performance of these sparse regularizers has not yet been fully evaluated for super-resolution imaging of compact objects with complicated structures on scales comparable to the diffraction limit, as is expected for black hole shadow imaging with the EHT. We will study this issue both for imaging with the full-complex visibility and closure quantities as an extension of this work and our previous work (Honma et al., 2014; Ikeda et al., 2016; Akiyama et al., 2017).

A relevant issue of our proposed methods is the computational time, most of which is spent in determining the optimum parameters for the regularization terms. Since we adopt 10-fold CV for determining regularization parameters, we need to do image reconstruction 11 times for each set of regularization parameters. This is not serious for imaging simulated data sets in this paper, which takes about a few hours in total for each Stokes parameter, although it would be a relevant issue for imaging larger data sets or wider field-of-views. Recently, Obuchi & Kabashima (2016) have proposed an accurate approximation of the validating error for $n$-fold CV for LASSO, which can be derived from the image reconstruction of full data sets. A similar approximation for TV regularization has been also derived very recently (Obuchi et al., 2016). These approximations may allow validating errors to be estimated by imaging the full data set just once at each set of regularization parameters. We will implement these estimators for our algorithm, which will significantly reduce the whole computational time ($\sim$ 10 times shorter for 10-fold CV). We will also work on optimizing and accelerating the MFISTA algorithms by parallel computing such as GPGPU (General-Purpose computing on Graphics Processing Units). This will be helpful for extending our works to wider-FOV imaging or imaging of much larger data sets with many more stations than VLBI networks, such as ALMA.

K.A. thanks Dr. Michael D. Johnson, Dr. Lindy Blackburn, Katherine L. Bouman and Andrew Chael for many fruitful discussions and constructive suggestions on this work. K.A. and this work are financially supported by the program of Postdoctoral Fellowships for Research Abroad at the Japan Society for the Promotion of Science (JSPS). M.P. acknowledges support from the NASA Massachusetts Space Grant Consortium and the National Science Foundation’s (NSF) Research Experiences for Undergraduates program. M.M. acknowledges support from the ERC Synergy Grant (Grant 610058). Event Horizon Telescope work at MIT Haystack Observatory and the Harvard-Smithsonian Center for Astrophysics is supported by grants from the NSF (AST-1440254, AST-1614868) and through an award from the Gordon and Betty Moore Foundation (GMBF-3561). Work on sparse modeling and Event Horizon Telescope at the Mizusawa VLBI Observatory is financially supported by the MEXT/JSPS KAKENHI Grant Numbers 24540242, 25120007 and 25120008.

Appendix A Monotonic FISTA algorithm

We show the details of the algorithms which were used to solve equation (7). Our algorithms are variations of the monotonic FISTA (MFISTA) which was introduced in Beck & Teboulle (2009a, b). We first show the general framework of MFISTA and discuss how we applied it for our problem.

A.1 General framework of MFISTA

The general form of the problem is defined as follows,

[TABLE]

where, $\bm{x}\in\Re^{n}$ and $C\subseteq\Re^{n}$ is some closed subset of $\Re^{n}$. The properties assumed for $f(\bm{x})$ and $g(\bm{x})$ are summarized below.

• •

$f(\bm{x}):\Re^{n}\rightarrow\Re$, is a convex function of $\bm{x}$. It is continuously differentiable, and the gradient $\nabla f(\bm{x})$ is Lipschitz continuous, where $L(f)$ denotes the Lipschitz constant of $\nabla f(\bm{x})$ .

• •

$g(\bm{x}):\Re^{n}\rightarrow(-\infty,\infty]$, is a convex function of $\bm{x}$. It is not necessarily smooth.

The pseudo code of MFISTA is summarized in algorithm 1.

If the upper bound of the Lipschitz constant $L(f)$ is known, $L$ is set to the upper bound and ${\mbox{\sc Initialize}}(L_{k-1},\bm{y}_{k})$ can be omitted. Otherwise we need to find an appropriate value of $L_{k}$. Let us define a function $Q(\bm{x},\bm{y};L)$ as follows

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the inner product. The function $p_{C}(\bm{y};L)$ is the proximal map which is defined as follows,

[TABLE]

The practical form of the proximal map depends on the definition of $g(\bm{x})$.

The procedure $\mbox{\sc Initialize}(L;\bm{x})$ is defined as follows using $Q(\bm{x},\bm{y};L)$ and $p_{C}(\bm{y};L)$.

MFISTA only uses the gradient of $f(\bm{x})$ and is known to have a fast convergence rate. Let $\bm{x}^{\ast}$ be the optimal point of the problem in equation (7). The MFISTA algorithm has the following property (Beck & Teboulle, 2009a, b),

[TABLE]

A.2 Applying MFISTA for Polarimetric Imaging

We explain how we applied MFISTA to solve equation (7).

For the Stokes $I$ image, ${\bf S}={\bf I}$ and $I_{i}\geq 0$. Since $\|{\bf I}\|_{1}=\sum_{i}I_{i}$ holds, we can apply MFISTA by defining $f(\cdot)$, $g(\cdot)$ and $C$ as follows,

[TABLE]

The form of the proximal map $p_{C}(\bm{y};L)$ for the case $g(\bm{x})=\|\bm{x}\|_{\rm tv}$ has been discussed in Beck & Teboulle (2009a). We used their Fast Projected Gradient (FGP) method restricting $I_{i}\geq 0$.

For the Stokes $Q$, $U$ and $V$ image, each component can take negative value. MFISTA can be applied by defining $f(\cdot)$ and $g(\cdot)$ as follows,

[TABLE]

The proximal map $p_{C}(\bm{y};L)$ for this case can also be realized by a slight modification of FGP.