11institutetext: Commonwealth Scientific and Industrial Research Organization
(CSIRO), Australia

# The application of compressive sampling to radio astronomy

II: Faraday rotation measure synthesis

F. Li 11 S. Brown 11 T. J. Cornwell 11 F. de Hoog 11

(Received date; accepted date)

Faraday rotation measure (RM) synthesis is an important tool to study and
analyze galactic and extra-galactic magnetic fields. Since there is a Fourier
relation between the Faraday dispersion function and the polarized radio
emission, full reconstruction of the dispersion function requires knowledge of
the polarized radio emission at both positive and negative square wavelengths
$\lambda^{2}$. However, one can only make observations for $\lambda^{2}>0$.
Furthermore observations are possible only for a limited range of wavelengths.
Thus reconstructing the Faraday dispersion function from these limited
measurements is ill-conditioned. In this paper, we propose three new
reconstruction algorithms for RM synthesis based upon compressive
sensing/sampling (CS). These algorithms are designed to be appropriate for
Faraday thin sources only, thick sources only, and mixed sources respectively.
Both visual and numerical results show that the new RM synthesis methods
provide superior reconstructions of both magnitude and phase information than
RM-CLEAN.

††offprints: F. Li

##  1 Introduction

The intrinsic polarization of a synchrotron emitting source together with
knowledge of propagation effects through intervening media provide critical
diagnostics for magnetic field orientation and fluctuations in a wide range of
astrophysical contexts. Faraday rotation is a physical phenomenon where the
position angle of linearly polarized radiation propagating through a magneto-
ionic medium is rotated as a function of frequency. As introduced in Brentjens
& de Bruyn (2005) and Heald (2009), Faraday rotation measure synthesis is an
important tool for analysing radio polarization data where multiple emitting
regions are present along a single line of sight. Observations of
extragalactic sources, which by necessity must be viewed through the Faraday
rotating and emitting Galactic interstellar-medium (de Bruyn et al. 2006;
Brown & Rudnick 2009; Schnitzeler et al. 2007, 2009), are an obvious example
of this regime. Burn (1966) introduced the Faraday dispersion function
$F(\phi)$, which describes the intrinsic polarized flux per unit Faraday depth
$\phi$ (in rad m-2), and its relationship with the complex polarized emission
$P(\lambda^{2})$ as

$P(\lambda^{2})=\int_{-\infty}^{\infty}F(\phi)\mathrm{e}^{2\mathrm{i}\phi\lambda^{2}}\;\mathrm{d}\phi,$
(1)

where $\lambda$ is the wavelength. Note that $P$ can also be written as
$P=Q+\mathrm{i}U$, where $Q$ and $U$ represent the emission of Stokes $Q$ and
Stokes $U$, respectively.

To study multiple emitting and Faraday rotating regions along each line of
sight, we need to reconstruct the Faraday dispersion function, which is, in
general, a complex-valued function of the Faraday depth $\phi$. From Eq. (1) ,
we can invert the expression to yield:

$F(\phi)=\frac{1}{\pi}\int_{-\infty}^{\infty}P(\lambda^{2})\mathrm{e}^{-2\mathrm{i}\phi\lambda^{2}}\;\mathrm{d}\lambda^{2}.$
(2)

However, the problem is that we can not observe the polarized emission at
wavelengths where $\lambda^{2}<0$. Even for the wavelength range
$\lambda^{2}>0$, it is impossible to observe all wavelengths or frequencies.
Brentjens & de Bruyn (2005) propose a synthesis method by first introducing an
observing window function $M(\lambda^{2})$. The observed complex polarized
emission can then be described as

$\widetilde{P}(\lambda^{2})=M(\lambda^{2}){P}(\lambda^{2}).$ (3)

In this paper, the tilde denotes the observed quantities.

If the observing window function is $M(\lambda^{2})$ with $m$ channels, the RM
spread function (RMSF) is be defined by

$R(\phi)=K\sum_{i=1}^{m}M(\lambda_{i}^{2})\mathrm{e}^{-2\mathrm{i}\phi(\lambda^{2}_{i}-\lambda^{2}_{0})},$
(4)

where the parameter $\lambda^{2}_{0}$ is the mean of the sampled values
between $\lambda^{2}_{1}$ and $\lambda^{2}_{m}$ within the observation window
$M(\lambda^{2})$; $i$ is the $i^{\rm{th}}$ channel in the observation window,
and $K$ is a normalising constant of the window function $M(\lambda^{2})$. In
this paper, we assume as a simplification that all channels have uniform
weights for the $m$ channels in the observing window function.

In Brentjens & de Bruyn (2005), the reconstructed Faraday rotation measure
synthesis can be written in discrete form as

$\widetilde{F}(\phi)\approx
K\sum_{i=1}^{m}\widetilde{P}(\lambda^{2}_{i})\mathrm{e}^{-2\mathrm{i}\phi(\lambda^{2}_{i}-\lambda^{2}_{0})},$
(5)

where $\widetilde{F}(\phi)$ is the reconstructed Faraday dispersion function.
From Eq. (5), we can see that the Faraday dispersion function can be
reconstructed provided that the spectral coverage is sufficient.

However, the reconstructed results generally include some side lobes. Using
the terminology of radio interferometry, the result of Brentjens & de Bruyn’
method is a dirty version of the Faraday dispersion function and is
abbreviated as “the dirty curve”. It is the convolution of $F(\phi)$ and the
RMSF, and a deconvolution step may be used to clean it up. By borrowing the
cleaning procedure in the image deconvolution method of Högbom CLEAN (Högbom
1974). Heald (2009) proposes the RM-CLEAN method which deconvolves
$\widetilde{F}(\phi)$ with the RMSF to remove the sidelobe response.

Recently, Frick et al. (2010) proposed a wavelet-based Faraday RM synthesis
method. In that approach, the authors assume specific magnetic field
symmetries in order to project the observed polarization emissions onto
$\lambda^{2}<0$.

Compressive sensing/sampling (CS) (Candès & Wakin 2008; Candès 2006; Wakin
2008) has been one of the most active areas in signal and image processing
over the last few years. Since CS was proposed, it has attracted very
substantial interest, and has been applied in many research areas (Wakin et
al. 2006; Lustig et al. 2007; Puy et al. 2010; Mishali et al. 2009; Bobin &
Starck 2009). In radio astronomy, CS has attracted attention as a tool for
image deconvolution. Wiaux et al. (2009a) compare the CS-based deconvolution
methods with the Högbom CLEAN method (Högbom 1974) on simulated uniform random
sensing matrices with different coverage rates. They apply compressive
sampling for deconvolution by assuming the target signal is sparse. Wiaux et
al. (2009b) proposed a new spread spectrum technique for radio interferometry
by using the non-negligible and constant component of the antenna separation
in the pointing direction. Recently, a new CS-based image deconvolution method
was introduced in Li et al. (2011) in which an isotropic undecimated wavelet
transform is adopted as a dictionary for sparse representation for sky images.

In this paper, we propose three new CS-based RM synthesis methods. In Section
2, the three CS-based RM synthesis methods are proposed. The implementation
details of the general experiment layout is given in Section 3. Simulation
results from the traditional methods are compared with those from CS-based
methods in Section 4. The final conclusions are given in Section 5.

##  2 CS-based RM synthesis

CS is primarily a sampling theory for sparse signals. A sensing matrix (Candès
et al. 2006b) is used to sample a signal with sparsity (few non-zero terms) or
a sparse representation with respect to a given basis function dictionary.
Given a limited number of measurements, generally less than the number of
unknowns in the target signal, the target signal can be reconstructed by
optimisation of an L1 norm. More information on the key concepts (such as
sparsity, incoherence, the restricted isometry property, and the L1 norm
reconstruction) and results can be found in Candès & Romberg (2007); Candès et
al. (2006a); Candès & Wakin (2008); Candès (2006).

CS includes two steps: sensing/sampling and reconstruction. This is in
contrast to Nyquist-Shannon theory which measures the target signal directly
without the reconstruction step. In this paper, we will focus on the
reconstruction step (calculating the Faraday dispersion function given an
observing window) rather than the sensing step (the selection of the observing
window), because the observing frequency range and the bandwidth for each
channel are usually fixed for a given telescope array.

To proceed with the CS approach, we rewrite the Fourier relationship as a
matrix equation. The projection of the Faraday dispersion function to the
polarized emission can be described as a matrix $Y$ of size $m\times N$

$\mathrm{Y}(j,N/2+k)=\mathrm{e}^{2\mathrm{i}\phi_{k}\lambda^{2}_{j}},j={1},\cdots,{m};k={1-N/2},\cdots,{N/2}.$
(6)

The inverse of the projection is the conjugate transpose of $Y$

$\mathrm{Y}^{\ast}(N/2+k,j)=\mathrm{e}^{-2\mathrm{i}\phi_{k}\lambda^{2}_{j}},j={1},\cdots,{m};k={1-N/2},\cdots,{N/2},$
(7)

where $\ast$ denotes the conjugate transpose. Suppose $f$ denotes the original
Faraday dispersion function $F(\phi)$ in a vector format of length $N$, then
the relationship between the Faraday dispersion function and the observed
radio emission is:

$\mathrm{Y}\mathbf{f}=\widetilde{\mathbf{p}},$ (8)

where $\widetilde{p}$ denotes the observed polarized emission in a vector
format of length $m$.

Because we can only measure a limited number of observations with the limited
number of channels, i.e. $m<<N$, there are many different potential Faraday
dispersion functions consistent with the measurements. To resolve these
ambiguities, the usual approach is to use some prior information to select a
solution. The prior information can be: the Faraday dispersion function is
real; the Faraday dispersion function has only point like signals which are
sparse in the Faraday depth domain or the Faraday dispersion function has a
sparse presentation with respect to a dictionary of basis functions, to name
just a few. Our three synthesis methods are based upon the last two structural
assumptions.

Before introducing our new RM synthesis methods, we need to review two
technical terms: Faraday thin and Faraday thick. A source can be either
Faraday thin if $\lambda^{2}\bigtriangleup\phi\ll 1$, or Faraday thick if
$\lambda^{2}\bigtriangleup\phi\gg 1$, where $\bigtriangleup\phi$ is the extent
of the source along the axis of Faraday depth $\phi$. Faraday thin sources can
be well described by Dirac $\delta$ function of $\phi$, while Faraday thick
sources have extensive support on the Faraday depth axis (Brentjens & de Bruyn
2005). Note that the definition of Faraday thin or thick is wavelength
dependent.

###  2.1 RM synthesis for Faraday thin sources

The relationship between the Faraday dispersion function and the observed
polarized radio emission is a Fourier pair if $\lambda^{2}=\pi u$, where $u$
is a wavelength related parameter. Since the space and Fourier domain are
perfectly incoherent (Candès & Romberg 2007), we can apply CS for RM synthesis
in a straightforward manner provided there are Faraday thin sources along the
line of sight since the screen is necessarily sparse.

In this context, CS recommends solving for the Faraday dispersion function by
minimising the L1 norm (summed absolute value) of the dispersion function as
inimising the L1 norm optimises the sparsity of the reconstruction. There
remains one further obstacle - the dispersion function is complex. We handle
this by summing the L1 norm of the real and imaginary parts:

${\rm{min}}\;\\{{\|\mathrm{Re}(\mathbf{f})\|_{l_{1}}+\|\mathrm{Im}(\mathbf{f})\|_{l_{1}}}\\}\;\;\;s.t.\;\mathrm{Y}\mathbf{f}=\widetilde{\mathbf{p}},$
(9)

where Re$(\bullet)$ and Im$(\bullet)$ denote the real and the imaginary parts,
respectively. By forming a real-valued vector of double length (comprising of
the real part and the imaginary part) of the complex-valued vector, almost all
L1 norm optimization solvers can be used for Eq. 9. This CS-based rotation
measure synthesis for Faraday thin sources is abbreviated as CS-RM-Thin. This
is similar in concept to RM-CLEAN because the assumption for RM-CLEAN is that
the Faraday dispersion function comprising of spike like signals. However,
results in Section 4 show that CS-RM-Thin can provides superior results to RM-
CLEAN.

###  2.2 RM synthesis for Faraday thick sources

CS-RM-Thin can work effectively when the Faraday dispersion function includes
Faraday thin sources along the line of sight. This limits its application for
the case when there are some Faraday thick sources along the line of sight.
However, CS can still reconstruct the Faraday dispersion function efficiently
provided that we can find a suitable dictionary of basis functions that can
decompose the extended sources into a sparse representation as described in
Candès (2006). In this paper, we adopt the Daubechies D8 wavelet transforms
(Daubechies 1992) as the dictionary. Other wavelet transforms can also be
adopted, the selection depends on the property of the Faraday dispersion
function. We choose the D8 wavelet transform, because we assume that the
Faraday dispersion function with thick sources is a sinc-like signal.

We can rewrite Eq. (8) as

$\mathrm{Y}\mathrm{W}^{-1}\alpha=\widetilde{\mathbf{p}},$ (10)

where $\mathrm{W}^{-1}$ is the inverse wavelet transform matrix of size
$N\times N$; $\alpha$ is the wavelet coefficient of the Faraday dispersion
function $\mathbf{f}$. The wavelet transform matrix is denoted as
$\mathrm{W}$, therefore, $\alpha=\mathrm{W}\mathbf{f}$. Other symbols follow
the definitions in Eq. (8). Under the condition that
$\mathrm{Y}\mathbf{f}=\widetilde{\mathbf{p}}$, we adopt the following
assumption: both the real part and the imaginary part of the Faraday thick
sources will have a sparse representation in the wavelet domain independently.
Then the wavelet based CS RM synthesis method for Faraday thick sources can be
written as

${\rm{min}}\;\\{{\|\mathrm{W}\cdot{\rm{Re}}(\mathbf{f})\|_{l_{1}}+\|\mathrm{W}\cdot{\rm{Im}}(\mathbf{f})\|_{l_{1}}}\\}\;\;\;s.t.\;\mathrm{Y}\mathbf{f}=\widetilde{\mathbf{p}}.$
(11)

This CS-based rotation measure synthesis for Faraday thick sources is
abbreviated as CS-RM-Thick.

###  2.3 RM synthesis for Faraday mixed sources

So far, we have proposed two RM synthesis methods: CS-RM-Thin and CS-RM-Thick
for solving Faraday thin sources and thick sources, respectively. However,
this begs the question: which method should be selected if there are both
Faraday thin sources and thick sources along the line of sight? Moreover, how
can we make a selection if we have no prior information about the Faraday
dispersion function, i.e. we are not sure what it looks like? Clearly, neither
of them is suitable, we therefore need another solution for solving the above
problems. Let us assume that there are both Faraday thin sources and Faraday
thick sources in $F(\phi)$. Suppose $\mathbf{f}_{\rm{thin}}$ denotes the
Faraday thin sources in $F(\phi)$ in a vector format of length $N$;
$\mathbf{f}_{\rm{thick}}$ denotes the Faraday thick sources in $F(\phi)$ in a
vector format of length $N$, then
$\mathbf{f}_{\rm{thin}}+\mathbf{f}_{\rm{thick}}=\mathbf{f}$. Eq.(8) can be
rewritten as

$\mathrm{Y}\mathbf{f}_{\rm{thin}}+\mathrm{Y}\mathbf{f}_{\rm{thick}}=\widetilde{\mathbf{p}}.$
(12)

Since we know that the L1 norm can preserve sparsity; Faraday thin sources
show sparsity in the Faraday depth domain; Faraday thick sources show sparsity
in the wavelet domain, we propose the following solution for the mixed
circumstance

$\displaystyle{\rm{min}}\;\\{{\|\mathrm{Re}(\mathbf{f}_{\rm{thin}})\|_{l_{1}}+\|\mathrm{Im}(\mathbf{f}_{\rm{thin}})}\|_{l_{1}}+\|{\mathrm{W}\cdot{\rm{Re}}(\mathbf{f}_{\rm{thick}})\|_{l_{1}}}$
$\displaystyle+\|\mathrm{W}\cdot{\rm{Im}}(\mathbf{f}_{\rm{thick}})\|_{l_{1}}\\}\;\;\;\;\;s.t.\;\mathrm{Y}\mathbf{f}=\widetilde{\mathbf{p}},$
(13)

where the definition of $\mathrm{W}$ is the same as the above subsection. The
above solution for Faraday mixed sources is still based on the spirit of CS by
preserving the sparsity in the Faraday depth domain for Faraday thin sources
and in the wavelet domain for Faraday thick sources, simultaneously. This CS-
based rotation measure synthesis for Faraday mixed sources is abbreviated as
CS-RM-Mix.

##  3 Implementation details

In this section, the implementation details of the above three proposed CS-
based RM synthesis methods will be given.

###  3.1 Preparation

To create a general experiment layout, we borrow some definitions and
conclusions from Brentjens & de Bruyn (2005). Both the diagram of the
wavelength square $\lambda^{2}$ and $\phi$ are displayed in figure 1. The
maximum observable Faraday depth $\|\phi_{\mathrm{max}}\|$ is given in
Brentjens & de Bruyn (2005)

$\|\phi_{\mathrm{max}}\|\approx\frac{\sqrt{3}}{\delta\lambda^{2}},$ (14)

where $\delta\lambda^{2}$ is the width of an observing channel. The full width
at half maximum (FWHM) of the main peak of the RMSF can be estimated by

$\delta\phi\approx\frac{2\sqrt{3}}{\bigtriangleup\lambda^{2}},$ (15)

where $\bigtriangleup\lambda^{2}$ is the width of the total $\lambda^{2}$
distribution.

Before using CS-based RM synthesis, the following two steps are needed:

  1. 1.

Select the resolution of the Faraday depth $\phi_{R}$. Since we know the
maximum observable Faraday depth $\phi_{\rm{max}}$ from Eq. (14) and the FWHM
of the main peak of the rotation measure spread function from Eq. (15), we can
select a grid resolution parameter $\phi_{R}$ in $\phi$ space, which should be
four or five times less than $\delta\phi$, to achieve Nyquist sampling.
However, for some observational window functions $M(\lambda^{2})$, the maximum
scale (Faraday thickness) that one is sensitive to, estimated by
$\frac{\pi}{\lambda^{2}_{\mathrm{min}}}$, is actually smaller than
$\delta\phi$. In these cases, it might be practical to Nyquist sample this
smaller scale in order to calculate $\phi_{R}$. Based on $\phi_{\rm{max}}$ and
$\bigtriangleup\phi$, we can calculate the number of grid points $N$ as

$N=\rm{floor}(\frac{2\phi_{\rm{max}}}{\phi_{R}}).$ (16)

  2. 2.

Constructing the two matrices $\mathrm{Y}$ and $\mathrm{Y}^{\ast}$.

Figure 1: This diagram shows the relationship between parameters in
$\lambda^{2}$ domain and $\phi$ domain, respectively.

The selection of these CS-based RM synthesis methods depends on the prior
knowledge about the Faraday dispersion function. If we assume it includes
Faraday thin sources only along the line of sight, we should select CS-RM-
Thin. On the other hand, CS-RM-Thick should be used if we know that there are
Faraday thick sources only. When we know that there are both Faraday thin
sources and thick sources along the line of sight, CS-RM-Mix should be used.
In most circumstances, we have no prior information about the Faraday
dispersion function, CS-RM-Mix can always be used to reconstruct a reliable
result as a compromise.

###  3.2 L1 norm solvers for CS-Based RM synthesis methods

For CS-RM-Thin and CS-RM-Thick, many optimization methods (Beck & Teboulle
2009; Becker et al. 2011; Boyd & Vandenberghe 2004) can be used to solve the
L1 norm minimization problem in a straightforward manner. There are many
solvers or toolboxes, for example, L1-Magic Matlab toolbox which can be
download from http://www.acm.caltech.edu/l1magic/. In this paper, L1-Magic is
adopted for solving equations (9) and (11). Fast Iterative Shrinkage-
Thresholding Algorithm (FISTA) (Beck & Teboulle 2009) can also be used for
solving this problem if we rewrite Eq. (9) or (11) in a Lagrangian form.

As far as CS-RM-Mix is concerned, the solvers or toolboxes introduced above
can be used for solving Eq. (2.3) indirectly. Suppose
$\alpha_{\mathrm{thick}}$ denotes the wavelet coefficients of the thick
sources $\mathbf{f_{\mathrm{thick}}}$ in the Faraday dispersion function in a
vector format, and $\mathrm{W}^{-1}$ is the inverse wavelet transform matrix,
then we have

$\mathbf{f_{\mathrm{{thick}}}}=\mathrm{W}^{-1}\alpha_{\mathrm{thick}}.$ (17)

Bring the above equation into Eq. (12), we have

$\mathrm{Y}\mathbf{f}_{\rm{thin}}+\mathrm{Y}\mathrm{W}^{-1}\alpha_{\mathrm{thick}}=\widetilde{\mathbf{p}}.$
(18)

Furthermore, Eq. (18) can be rewritten as

$\left[\mathrm{Y}\;\mathrm{Y}\right]\left[\begin{array}[]{cc}\mathrm{I}&\mathrm{O}\\\
\mathrm{O}&\mathrm{W}^{-1}\end{array}\right]\left[\begin{array}[]{c}\mathbf{f}_{\mathrm{thin}}\\\
\alpha_{\mathrm{thick}}\end{array}\right]=\widetilde{\mathbf{p}}_{m\times 1},$
(19)

where $\mathrm{I}$ denotes the identity matrix of size $N\times N$, and
$\mathrm{O}$ is the matrix of all zeros with the size of $N\times N$. If we
denotes $\mathrm{Y_{mix}}=[\mathrm{Y}\;\mathrm{Y}]_{m\times 2N}$,
$\mathrm{T}$=$\left[\begin{array}[]{cc}\mathrm{I}&\mathrm{O}\\\
\mathrm{O}&\mathrm{W}^{-1}\end{array}\right]_{2N\times 2N}$ and
$\mathbf{c}=\left[\begin{array}[]{c}\mathbf{f}_{\mathrm{thin}}\\\
\alpha_{\mathrm{thick}}\end{array}\right]_{2N\times 1}$, almost all L1 norm
minimization solvers can be used to solve Eq. (2.3) with:

(20)

To help readers who are unfamiliar with L1 norm minimization to use or
implement our proposed CS-based RM synthesis methods, we have developed a
simple algorithm based on the iterative soft-thresholding algorithm (ISTA)
(Beck & Teboulle 2009) for CS-RM-Mix as an example. The algorithm is as
follows:

  1. 1.

Initialization:

    1. (a)

Choose parameters: the soft-threshold $\tau$ (this can be set by 1 for most
circumstances) , the stopping-threshold $\delta$ (this can be set by the noise
level)

    2. (b)

Set the number of iteration $l=\rm{floor}(\tau/\delta)$

    3. (c)

$\mathbf{f}_{\rm{thin}}=\mathrm{Y}^{\ast}\widetilde{\mathbf{p}}$;
$\mathbf{f}_{\rm{thick}}=0$

  2. 2.

Within $l$ iterations:

    1. (a)

Reconstructing the Faraday thin sources

      1. i.

Calculate the residual
$\mathbf{r}=\widetilde{\mathbf{p}}-\mathrm{Y}\mathbf{f}_{\rm{thin}}-\mathrm{Y}\mathbf{f}_{\rm{thick}}$

      2. ii.

Calculate the gradient $\mathbf{d}=\mathrm{Y}^{\ast}\mathbf{r}$

      3. iii.

Update $\mathbf{f}_{\rm{thin}}=\mathbf{f}_{\rm{thin}}+\mathbf{d}$

      4. iv.

Soft threshold ${\rm{Re}}(\mathbf{f}_{\rm{thin}})$ and
${\rm{Im}}(\mathbf{f}_{\rm{thin}})$, respectively. Set any values below $\tau$
to zero and update $\mathbf{f}_{\rm{thin}}$

    2. (b)

Reconstructing the Faraday thick sources

      1. i.

Calculate the residual
$\mathbf{r}=\widetilde{\mathbf{p}}-\mathrm{Y}\mathbf{f}_{\rm{thin}}-\mathrm{Y}\mathbf{f}_{\rm{thick}}$

      2. ii.

Calculate the gradient $\mathbf{d}=\mathrm{Y}^{\ast}\mathbf{r}$

      3. iii.

Update $\mathbf{f}_{\rm{thick}}=\mathbf{f}_{\rm{thick}}+\mathbf{d}$

      4. iv.

Calculate the wavelet coefficients for both the real part and imagery part of
$\mathbf{f}_{\rm{thick}}$, i.e.
$\mathrm{W}\cdot{\rm{Re}}(\mathbf{f}_{\rm{thick}})$ and
$\mathrm{W}\cdot{\rm{Im}}(\mathbf{f}_{\rm{thick}})$

      5. v.

Soft threshold the wavelet coefficients of
$\mathrm{W}\cdot{\rm{Re}}(\mathbf{f}_{\rm{thick}})$ and
$\mathrm{W}\cdot{\rm{Im}}(\mathbf{f}_{\rm{thick}})$. Set any values below
$\tau$ to zero and update $\mathrm{W}\cdot{\rm{Re}}(\mathbf{f}_{\rm{thick}})$
and $\mathrm{W}\cdot{\rm{Im}}(\mathbf{f}_{\rm{thick}})$

      6. vi.

Calculate the inverse wavelet transform for both the real part and imaginary
part, respectively, then update $\mathbf{f}_{\rm{thick}}$

    3. (c)

$\tau=\tau-\delta$

  3. 3.

Reconstructed Faraday dispersion function
$\widetilde{\mathbf{f}}=\mathbf{f}_{\rm{thin}}+\mathbf{f}_{\rm{thick}}$

Note the Faraday thin sources and thick sources are reconstructed separately.
This can be helpful when astronomers focus on either Faraday thin sources or
thick sources only.

If we ignore the step 2b and set $\mathbf{f}_{\rm{thick}}=0$, the above
algorithm will degenerate to CS-RM-Thin. On the contrary, if we ignore the
step 2a and set $\mathbf{f}_{\rm{thin}}=0$, the algorithm will degenerate to
CS-RM-Thick. In this paper, these CS-based rotation measure synthesis methods
are implemented in MATLAB. Our code may be found at
http://code.google.com/p/csra/downloads11footnotetext: Download the file
“CS_RM.zip” which includes both CS-RM-Thin and CS-RM-Thick algorithms.

##  4 Experimental results

We have adopted the standard test platform in Brentjens & de Bruyn (2005). In
this platform there 126 observing channels within Window 1 (0.036 to 0.5m)
evenly distributed in $\lambda^{2}$. Three different Faraday dispersion
functions are simulated to test these CS-Based RM synthesis methods for
Faraday thin sources, Faraday thick sources and mixed sources cases.

###  4.1 Simulation results for Faraday thin sources

We simulate a Faraday dispersion function containing four Faraday thin
sources. See figure 2 for the function. From left to right, these sources are:
$F(-10)=10-4\mathrm{i}$ Jy m2 rad-1, $F(-17)=-7+5\mathrm{i}$ Jy m2 rad-1,
$F(40)=9-7\mathrm{i}$ Jy m2 rad-1 and $F(88)=-4+3\mathrm{i}$ Jy m2 rad-1. The
true Faraday dispersion function is shown in the top left corner of figure 2.
The thin solid line shows the real value, the dashed line the imaginary part,
and the thick solid line the amplitude. The Faraday dispersion function in
this test is complex valued, i.e. the intrinsic polarization angles are non-
zero. We have selected this simulated dispersion function rather than the
standard test (real valued Faraday dispersion function) in Brentjens & de
Bruyn (2005) to investigate the behaviour when the intrinsic polarization
angles are non-zero. From Eq. 15, we can calculate that the FWHM of the RMSF
is around $14$ rad m2. The dirty curve which is calculated by assuming all
unmeasured emission to be zeros, is shown in figure 2 (b). The result of RM-
CLEAN is shown in figure 2 (c). Note that RM-CLEAN cannot correctly
reconstruct the magnitude or phase of the Faraday dispersion function in this
case. It has been demonstrated previously that RM-CLEAN has difficulty with
the separation of sources below the FWHM (Farnsworth et al. 2011). In figure 2
(a), from left to right, the distance between the first and the second sources
is smaller than the FWHM. In figure 2 (c), the first source and second source
are merged to form an unrealistic source. Results of CS-RM-Thin, CS-RM-Thick
and CS-RM-Mix are shown in the figure 2 (d), (e) and (f), respectively. As one
might expect, the result of CS-RM-Thick is poor, because the Faraday
dispersion function has no Faraday thick sources. In figure 2, though CS-RM-
Mix does a better job than CS-RM-Thick and RM-CLEAN in terms of magnitude and
phase of the reconstructed disperse functions, the result is far from good
enough. The result of CS-RM-Thin is shown in figure 2) (d). We can see that
CS-RM-Thin gives the best result, reconstructing the Faraday dispersion
function without any error. This is consistent with CS theory which says that
we can reconstruct the sparse signal exactly with “overwhelming probability”
(Candès & Wakin 2008; Candès 2006; Wakin 2008). In this test, we select
$\phi_{R}=3.6$ rad m2 which is around one quarter of FWHM, so $N=480$.

To carry out a numerical comparison, we use the root mean square (RMS) error
to characterise the difference between the reconstructed
$\widetilde{\mathbf{f}}$ and the Faraday dispersion function $\mathbf{f}$:

$RMS=\sqrt{\frac{\sum_{-N/2+1}^{N/2}(\mathbf{f}-\widetilde{\mathbf{f}})^{2}}{N}}.$
(21)

The RMS error is calculated for all the candidate methods, and the results are
listed in table 1. Results for this test can be found in the first row of the
table. CS-RM-Thin gives the best result.

Figure 2: We have tested our methods on a Faraday dispersion function with four Faraday thin sources. From left to right in the first row are: (a) Original $F(\phi)$, (b) Dirty curve, (c) RM-CLEAN. From left to right in the second row are: (d) CS-RM-Thin, (e) CS-RM-Thick, (f) CS-RM-Mix. The thin solid line shows the real value, the dashed line the imaginary part, and the thick solid line the amplitude. All horizontal axis units are rad m-2, i.e. $\phi$, and all vertical axis units are Jy m2 rad-1.  Figure 3: Reconstructed results of a Faraday dispersion function with two Faraday thick sources. From left to right in the first row are: (a) Original $F(\phi)$, (b) Dirty curve, (c) RM-CLEAN. From left to right in the second row are: (d) CS-RM-Thin, (e) CS-RM-Thick, (f) CS-RM-Mix. Table 1: Numerical comparison results (RMS error) | Dirty Faraday dispersion function | RM-CLEAN | CS-RM-Thin | CS-RM-Thick | CS-RM-Mix
---|---|---|---|---|---
Test with Faraday thin sources | 1.40 | 0.78 | 0.00 | 0.84 | 0.76
Test with Faraday thick sources | 2.18 | 0.91 | 1.07 | 0.72 | 0.77
Test with Faraday mix sources | 2.45 | 1.03 | 0.95 | 0.81 | 0.80

###  4.2 Simulation results for Faraday thick sources

We now test our CS-based methods for a Faraday dispersion function with
Faraday thick sources. Here we assume that the Faraday dispersion function
includes two sources $F(\phi)=2-2\mathrm{i}$ Jy m2 rad-1 where
$-120\leq\phi\leq 40$ and $F(\phi)=-6-3\mathrm{i}$ Jy m2 rad-1 where
$30\leq\phi\leq 70$. The simulated Faraday dispersion function is shown in
figure 3 (a). We adopt the previous observing window for this test and the
dirty Faraday dispersion function is shown in figure 3 (b). Note that this
only provides us with the approximate shape of the Faraday dispersion
function. As mentioned in Frick et al. (2010), the magnitude of $F(\phi)$
indicates the polarized emission of the region with Faraday depth $\phi$ and
its phase defines the intrinsic position angle. For the study of polarized
emission of galaxies, the magnitude of $F(\phi)$ is very important, and for
the study of orientation of the magnetic field perpendicular to the line of
sight, the phase information of $F(\phi)$ is crucial. Unfortunately, Brentjens
& de Bruyn’ method can not reconstruct reliable phase information for
$F(\phi)$. The cleaned version is shown in figure 3 (c). RM-CLEAN also failed
to reconstruct the phase information. CS-RM-Thin does not work well for this
case, because $F(\phi)$ is not sparse. The result of CS-RM-Thin is shown in
figure 3 (d). CS-RM-Mix is also used for this test by assuming that there are
both Faraday thin sources and thick sources. The result of CS-RM-Mix is shown
in figure 3 (f). The reconstructed result from CS-RM-Thick is shown in figure
3 (e). Even though CS-RM-Mix gives a much better result than CS-RM-Thin and
RM-CLEAN, the result is not as good as that of CS-RM-Thick. CS-RM-Thick
provides the best approximation to the original Faraday dispersion function in
terms of both magnitude and phase. This is also supported by the numerical
comparison in table 1. From the second row of the table, we can see that CS-
RM-Thick gives the smallest RMS error $0.72$ which is slightly better than
that of CS-RM-Mix $0.77$.

###  4.3 Simulation results for Faraday mixed sources

So far we have tested our methods for both Faraday thin sources and Faraday
thick sources, then we will test our CS-based methods for the mixed
circumstances - both Faraday thin sources and thick sources. The simulated
Faraday dispersion function for this test is shown in the top left corner of
figure 4. There are three sources along the line of sight. From left to right,
these sources are: $F(-58)=-4+3\mathrm{i}$ Jy m2 rad-1,
$F(-30)=10-4\mathrm{i}$ Jy m2 rad-1, $F(\phi)=2-6\mathrm{i}$ Jy m2 rad-1 where
$41\leq\phi\leq 100$. The original Faraday dispersion function is shown in the
top left corner of figure 2. We use the previous observing window for this
test. Based on the observing window (with 126 observing channels distributed
between 0.036 m to 0.5 m), the RMSF is calculated and shown in figure 4 (b).
the dirty curve is shown in figure 4 (c). The cleaned version by RM-CLEAN is
shown in figure 4 (d). RM-CLEAN performs badly in this test, because there is
a Faraday thick source. In general, RM-CLEAN can only work well when there are
Faraday thin sources along the line of sight. The results of CS-RM-Thin and
CS-RM-Thick are shown in figure 4 (e) and (f), respectively. We can see that
CS-RM-Thin reconstructs the two Faraday thin sources nicely, but failed to
reconstruct the Faraday thick source. On the contrary, CS-RM-Thick can
properly reconstruct the Faraday thick source, but expends the two Faraday
thin sources by mistake. As introduced above, CS-RM-Mix can separate the
Faraday thin components $\mathbf{f}_{\rm{thin}}$ and the Faraday thick
components $\mathbf{f}_{\rm{thick}}$ during the reconstruction. In this test,
the soft-threshold $\tau=1$ and $\delta=0.001$ in the proposed algorithm for
CS-RM-Mix. The results of CS-RM-Mix are shown in the last row of figure 4. The
separated Faraday components $\mathbf{f}_{\rm{thin}}$ and
$\mathbf{f}_{\rm{thick}}$ are shown in figure 4 (g) and (h), respectively. The
sum of the separated components is the reconstructed Faraday dispersion
function which is shown in figure 4 (i). We can see that the result of CS-RM-
Mix takes advantage of the results of both CS-RM-Thin and CS-RM-Thick, and it
gives the closest approximation to the original $F(\phi)$. From objective
evaluation point of view, CS-RM-Mix gives the minimum RMS error $0.80$ which
can be seen from the third row of table 1.

Figure 4: Reconstructed results of a Faraday dispersion function with two
Faraday thin sources and a thick source. From left to right in the first row
are: (a) Original $F(\phi)$, (b) RMSF of the observing window with 126
observing channels distributed between 0.036 m to 0.5 m, (c) Dirty curve. From
left to right in the second row are: (d) RM-CLEAN, (e) CS-RM-Thin, (f) CS-RM-
Thick. From left to right in the third row are: (g) Thin components
$\mathbf{f}_{\rm{thin}}$ by using CS-RM-Mix, (h) Thick components
$\mathbf{f}_{\rm{thick}}$ by using CS-RM-Mix, (i) CS-RM-Mix i.e.
$\mathbf{f}_{\rm{thin}}+\mathbf{f}_{\rm{thick}}$. All horizontal axis units
are rad m-2, and all vertical axis units are Jy m2 rad-1.

###  4.4 Discussion

In the figures above, we show the reconstructed dispersion function without
smoothing and addition of the residuals - a step commonly know as restoring.
For real applications, restoration is an option if the robustness is
insufficient.

From the above three tests, we can see that there is no single CS-based RM
synthesis method with an outstanding performance for any circumstances. The
best reconstruction can only be achieved when we have some prior knowledge
about the Faraday dispersion function and select the relevant CS-based RM
synthesis method. If such information is not available, CS-RM-Mix can always
be used as a compromise. Another option is that we can either select CS-RM-
Thin with a large $\phi_{R}$ or CS-RM-Thick with a small $\phi_{R}$. If a
large $\phi_{R}$ is selected, $\mathbf{f}$ is likely to be a sparse vector, so
CS-RM-Thin should be selected for the reconstruction. On the other hand, a
small $\phi_{R}$ can expand compact sources into extended sources, so CS-RM-
Thick will become suitable for the reconstruction. It does not mean that CS-
RM-Thick can be used for any cases with a small $\phi_{R}$, because a smaller
$\phi_{R}$ (which means a larger $N$ from Eq. 16) brings more unknowns in
$\mathbf{f}$ and more uncertainty. We have to balance these factors.

For RM synthesis, the observing window is quite similar to the frequency
filters. For example, the previously introduced observing Window 1 is like a
low pass filter in the wavelength squared domain. Therefore, we should bear in
mind, to observe radio sources with Faraday thick sources under the same
restriction of $m$ and $\delta\lambda^{2}$, the higher frequency observing
band the better.

These CS-based RM synthesis methods are not limited to the optimisation
methods L1-Magic solver, FISTA and ISTA. Other L1 norm optimization solvers
can also be adopted for solving Eqs. (9), (11) and (2.3). Though the wavelet
transform is used as the sparse representations dictionary in this paper,
there are some other potential basis functions can be used to achieve sparsity
for the Faraday thick sources.

In summary, the performance of CS-based RM synthesis methods depend on the
observing window, the resolution of $\phi$, the number of measurements, and
the sparsity of the Faraday dispersion function. The reconstruction of these
CS-based RM synthesis methods take less time than RM-CLEAN in general. For
example, CS-RM-Thin takes 3 seconds for the above tests in contrast with 5
seconds of RM-CLEAN. The calculation time really depends on the construction
of the matrix $\mathrm{Y}$, the larger $N$ and $m$ ($N=480$, $m=126$ for the
above tests), the more time it takes. The computer is a 2.53-GHz Core 2 Duo
MacBook Pro with 4GB RAM.

##  5 Conclusions

Faraday rotation measure synthesis is a very useful tool to study
astrophysical magnetic fields. The problem in RM synthesis is to reconstruct
the Faraday dispersion function given incomplete observations. From CS, we
know that a signal with sparsity can be well reconstructed based on few
measurements. We propose three CS-based RM synthesis methods by finding sparse
representations of the Faraday dispersion functions $F(\phi)$ for different
circumstances. CS-RM-Thin, CS-RM-Thick and CS-RM-Mix can be used for Faraday
thin sources only, thick sources only and mixed sources, respectively. In
general, Faraday thin sources show sparsity in the Faraday depth domain
$\phi$, therefore, we apply the CS reconstruction methods (L1 norm
optimization solvers) in a straightforward manner i.e. CS-RM-Thin. Although
Faraday thick sources are not sparse in the Faraday depth domain, they are
sparse in the wavelet domain for a suitably chosen basis wavelet. Therefore,
we apply the L1 norm optimization solvers in the wavelet domain i.e. CS-RM-
Thick. When there are Faraday mixed sources along the line of sight, we
preserve the sparsity by using L1 norm in the Faraday depth domain and the
wavelet domain simultaneously i.e. CS-RM-Mix.

As shown in the experimental results, the performance of these CS-based
methods is markedly superior to the traditional RM synthesis methods
(Brentjens & de Bruyn 2005; Heald 2009) in terms of magnitude and angle of the
reconstructed Faraday dispersion function. Exemplified by figure 2, both
Brentjens & de Bruyn’ method and RM-CLEAN do not work well in disentangling
two closely spaced sources. In contrast, CS-RM-Thin can separate the sources.

##  6 Acknowledgements

We thank Jean-Luc Starck for early discussions on Compressive Sampling. In
addition, comments, suggestions and derivations made by Dr Brentjens during
the review process are very much appreciated.

## References

  * Beck & Teboulle (2009) Beck, A. & Teboulle, M. 2009, SIAM Journal on Imaging Sciences, 2, 183
  * Becker et al. (2011) Becker, S., Bobin, J., & Candès, E. 2011, SIAM Journal on Imaging Sciences, 4, 1
  * Bobin & Starck (2009) Bobin, J. & Starck, J.-L. 2009, Wavelets XIII, 7446, 74460I
  * Boyd & Vandenberghe (2004) Boyd, S. & Vandenberghe, L. 2004, Convex optimization (Cambridge University Press), 716
  * Brentjens & de Bruyn (2005) Brentjens, M. & de Bruyn, A. 2005, Astronomy and Astrophysics, 441, 1217
  * Brown & Rudnick (2009) Brown, S. & Rudnick, L. 2009, The Astronomical Journal, 137, 3158
  * Burn (1966) Burn, B. 1966, Monthly Notices of the Royal Astronomical Society, 133, 67
  * Candès (2006) Candès, E. 2006, Int. Congress of Mathematics, Madrid, Spain, 2006, 3, 1433
  * Candès & Romberg (2007) Candès, E. & Romberg, J. 2007, Inverse Problems, 23, 969
  * Candès et al. (2006a) Candès, E., Romberg, J., & Tao, T. 2006a, IEEE Transactions on information Theory, 52, 489
  * Candès et al. (2006b) Candès, E., Romberg, J., & Tao, T. 2006b, Communications on Pure and Applied Mathematics, 59, 1207
  * Candès & Wakin (2008) Candès, E. & Wakin, M. 2008, IEEE Signal Processing Magazine, 25, 21 30
  * Daubechies (1992) Daubechies, I. 1992, Ten lectures on wavelets (Society for Industrial and Applied Mathematics Philadelphia, PA, USA)
  * de Bruyn et al. (2006) de Bruyn, A., Katgert, P., Haverkorn, M., & Schnitzeler, D. 2006, Astronomische Nachrichten, 327, 487
  * Farnsworth et al. (2011) Farnsworth, D., Rudnick, L., & Brown, S. 2011, Arxiv preprint: 1103.4149
  * Frick et al. (2010) Frick, P., Sokoloff, D., Stepanov, R., & Beck, R. 2010, Monthly Notices of the Royal Astronomical Society: Letters, 401, L24
  * Heald (2009) Heald, G. 2009, Cosmic Magnetic Fields: From Planets, 259, 591
  * Högbom (1974) Högbom, J. 1974, Astronomy and Astrophysics Supplement, 15, 417
  * Li et al. (2011) Li, F., Cornwell, T., & de Hoog, F. 2011, A&A, 528, A31
  * Lustig et al. (2007) Lustig, M., Donoho, D., & Pauly, J. 2007, Magnetic Resonance in Medicine, 58, 1182
  * Mishali et al. (2009) Mishali, M., Eldar, Y. C., Dounaevsky, O., & Shoshan, E. 2009, CCIT Repor No. 751
  * Puy et al. (2010) Puy, G., Wiaux, Y., Gruetter, R., Thiran, J., & De, D. V. 2010, IEEE International Symp. on Biomedical Imaging: From Nano to macro
  * Schnitzeler et al. (2007) Schnitzeler, D., Katgert, P., & de Bruyn, A. 2007, Astronomy and Astrophysics, 471, L21
  * Schnitzeler et al. (2009) Schnitzeler, D., Katgert, P., & de Bruyn, A. 2009, Astronomy and Astrophysics, 494, 611
  * Wakin (2008) Wakin, M. 2008, IEEE Signal Processing Magazine, 1623
  * Wakin et al. (2006) Wakin, M., Laska, J., Duarte, M., & Baron, D. 2006, International Conference on Image Processing, ICIP 2006, Atlanta George, Oct 2006, 1437
  * Wiaux et al. (2009a) Wiaux, Y., Jacques, L., Puy, G., Scaife, A., & Vandergheynst, P. 2009a, Monthly Notices of the Royal Astronomical Society, 395, 1733
  * Wiaux et al. (2009b) Wiaux, Y., Jacques, L., Puy, G., Scaife, A. M. M., & Vandergheynst, P. 2009b, Monthly Notices of the Royal Astronomical Society, 395, 1733