• Outline
• Glossary
• 8. Calibration
  • Previous: 8. Calibration
  • Next: 8.2 1GC calibration

Import standard modules:

Import section specific modules:

8.1 Calibration as a Least Squares Problem

In this section we discuss the procedure that is generally used in practice to perform calibration. We will use the unpolarized RIME in this section instead of the full-polarized RIME (see $\S$ 7 ➞). It provides us with a much simpler framework with which we can grasp the basics of calibration. Moreover, we assume for the sake of simplicity that the observed data are only corrupted by the instrument's antenna gains. This assumption results in a idealised treatment as there are many other factors that do in fact corrupt radio interferometric data (see $\S$ 7 ➞).

The unpolarized RIME is given by the following:

Unpolarized RIME
\begin{equation} d_{pq}(t) = g_p(t) g_q^*(t) \tilde{d}_{pq}(t) + \epsilon_{pq}(t), \end{equation}

where $d_{pq}(t)$ and $\tilde{d}_{pq}(t)$ denote the corrupted observed and uncorrupted visibility at time $t$ associated with baseline $pq$. Moreover, the factors $g_p$ and $g_q$ denote the complex gain of antenna $p$ and $q$. The term $\epsilon_{pq}$ is a zero mean (Gaussian) noise term, representing thermal noise.

We begin this section by generating the $uv$-tracks of a fictitious instrument in $\S$ 8.1.1 ⤵. In $\S$ 8.1.2 ⤵ we phrase the calibration problem (for the antenna gains) as a least squares minimization problem. Then in $\S$ 8.1.3 ⤵ we simulate "realistic" visibility data for the $uv$-tracks by including gain errors and adding noise to the resulting visibilities (similar to adding noise to a simple sinusoid as seen in $\S$ 2.11 ➞). We then vectorize the problem in $\S$ 8.1.4 ⤵, enabling us to use the built in scipy Levenberg-Marquardt algorithm to calibrate the data produced in $\S$ 8.1.3 ⤵. We implement the aforementioned steps via a wrapper ipython function called create_G_LM. We finish $\S$ 8.1.4 ⤵ by using create_G_LM to estimate the antenna gains corrupting the simulated data we produced in $\S$ 8.1.3 ⤵. The estimated antenna gains are then used to correct the corrupted data in $\S$ 8.1.5 ⤵.

8.1.1 Creating $uv$-Tracks: East-West Interferometer

We know from $\S$ 4.4.1.B.3 ➞ that when we work with an east-west interferometer things simplify to a large degree. Firstly: $XYZ = [0~|\mathbf{b}|~0]^T$, where $|\mathbf{b}|$ is the baseline length.

Moreover, we have that:

$uv$-Coverage of an EW-array (8.1)
\begin{eqnarray} \\ u &=&| \mathbf{b}|\cos H\\ v &=& |\mathbf{b}|\sin H \sin \delta, \end{eqnarray}

where $H$ is the hour angle of the field center and $\delta$ its declination. In this section we will be plotting the $uv$-coverage of a three element east-west interferometer. The ENU layout of a simple interferometer is given below. Note that $|\mathbf{b}|$ is measured in wavelengths.

Now consider an array made up of three antennas situated at 0, 100, 200 meters east of the array center as shown in the code fragment below.

We first need to set the hour angle range of our observation and the declination of our field center.

Our hour angle range is from -6h to 6h, and our declination is set to $60^{\circ}$.

As the earth rotates the antennas trace out $uv$-tracks (ellipses) as shown in the code fragment below, where the red tracks are due to baseline $pq$ and blue tracks are due to baseline $qp$. We can construct these $uv$-tracks by using Eq. 8.1 ⤵.

We can also pack the $uv$-coverage into a 2D-matrix. We denote the rows of this matrix with $p$ and the columns with $q$. The $pq$-th entry denotes the $uv$-track associated with baseline $pq$. The reason for packing the visibilities into a 2D structure will become apparent in Sec. 8.1.2 ⤵.

8.1.2. Unpolarized Calibration

As explained in $\S$ 7.2 ➞ the RIME assumes that our observed signal is polarized. For the sake of simplicity, however, we will now introduce the calibration problem with the underlying assumption that the observed signal is unpolarized. Unpolarized calibration is achieved by solving the following minimization problem:

Unpolarized Calibration
\begin{equation} \min_{\boldsymbol{\mathcal{G}}} \left \| \boldsymbol{\mathcal{D}} - \boldsymbol{\mathcal{G}}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{G}}^H \right \|, \end{equation}

where

• $\boldsymbol{\mathcal{D}}$ is the observed visibility matrix. Each entry, which we denote by $d_{pq}$, represents the visibility measured by the baseline formed by antennas $p$ and $q$.
• $\boldsymbol{\mathcal{M}}$ is the model visibility matrix. The entry $m_{pq}$ of $\boldsymbol{\mathcal{M}}$ denotes a true or model visibility which was created with the calibration sky model and a $uv$-point on the $uv$-track associated with baseline $pq$.
• $\boldsymbol{\mathcal{G}} = \textrm{diag}(\mathbf{g})$ is the antenna gain matrix, where $\mathbf{g}=[g_1,g_2,\cdots,g_N]^T$ denotes the antenna gain vector. The operator $\textrm{diag}(\cdot)$ forms a diagonal matrix from a vector by putting the elements of the vector on the main diagonal. The vector $\mathbf{g}$ represents the instrumental response of the antennas, i.e. the complex antenna gains. These antenna gains are chosen in such a way that they minimize the difference between the observed and model visibilities.
• $\boldsymbol{\mathcal{G}}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{G}}^H$ is the predicted visibility matrix. This matrix contains the model visibilities after the antenna gains have been applied to them.

The superscript $(\cdot)^H$ denotes the Hermitian or conjugate transpose and $\left \| \cdot \right \|$ denotes the norm used. Most calibration algorithms use the Frobenius norm for matrices and the 2-norm or Euclidean norm for vectors, thus treating calibration as a least squares problem.

8.1.3. Creating an Unpolarized Visibility Matrix (create_vis_mat)

In this section we present a function that allows us to create the observed visibility matrix $\boldsymbol{\mathcal{D}}$ and the model visibility matrix $\boldsymbol{\mathcal{M}}$. The function employs three seperate steps to produce a visibility matrix, namely:

• We first take the Fourier transform of the skymodel and then sample the result using the sampling function (i.e. $uv$-coverage). The sky model can only consist of point sources. Mathematically we may represent our sky model as $I(l,m) = \sum_k A_k\delta(l-l_k,m-m_k)$, where $A_k$ denotes the flux of the $k$-th source and $(l_k,m_k)$ denotes the direction cosine position vector that is associated with the $k$-th source. We then have that $V(u,v) = \mathscr{F}\{I(l,m)\} = A_k e^{-2\pi \imath l_k u + m_k v}$, where $\mathscr{F}\{\cdot\}$ denotes the Fourier transform of its operand. This result stems from the fact that the Fourier transform of a delta function is a complex exponential. If we now apply the sampling function we finally obtain $V_{pq}(u_{pq},v_{pq}) = A_k e^{-2\pi \imath l_k u_{pq} + m_k v_{pq}}$. We now use $V_{pq}$ to construct a 2D model visibility matrix. The skymodel is passed to the function via the variable point_sources. The sampling function is passed in via u_m and v_m.
• We then corrupt the visibilities with the antenna gains that were passed into the function via g. We use g to construct $\boldsymbol{\mathcal{G}}$. We corrupt our visibilites by multiplying by $\boldsymbol{\mathcal{G}}$ on the left of the model visibility matrix and on the right by $\boldsymbol{\mathcal{G}}^H$.
• The last step is to add some noise to our visibilities. The standard deviation of the noise is passed in via sig.

It should now be obvious how we can use the same function to produce both $\boldsymbol{\mathcal{M}}$ and $\boldsymbol{\mathcal{D}}$. In the case of $\boldsymbol{\mathcal{M}}$, we do not corrupt our visibilities, nor add any noise. See the function create_vis_mat below.

We now use create_vis_mat to create an example $\boldsymbol{\mathcal{M}}$ and $\boldsymbol{\mathcal{D}}$. Note that there are two sources in our sky model.

We now plot the baseline entries of $\boldsymbol{\mathcal{M}}$ and $\boldsymbol{\mathcal{D}}$.

The images above contain the real part of the corrupted (green) and uncorrupted (blue) visibilities as a function of timeslots for baseline 01, 02 and 12 respectively.

8.1.4 Levenberg-Marquardt (create_G_LM)

We are now ready to use least squares to calibrate $\boldsymbol{\mathcal{D}}$ (see GPU accelerated nonlinear optimization in radio interferometric calibration ⤴).

We first present a brief review of least squares minimization. Suppose we wish to fit a model $\mathbf{f}\left( \mathbf{m},\breve{\mathbf{g}}\right)$, where $\mathbf{m}$ and $\breve{\mathbf{g}}$ denote the model input values and a vector of unknown variables respectively, to some data $\left\{\mathbf{d}_{i},\mathbf{m}_{i}\right\}$. The vector of unknown variables $\breve{\mathbf{g}}$ parametrize the model. A standard method for determining which parameter vector $\breve{\mathbf{g}}$ best fits the data is to minimize the sum of the squared residuals. This technique is referred to as least squares minimization. The residual vector is denoted by $\mathbf{r}(\mathbf{m},\mathbf{d},\breve{\mathbf{g}}) = \mathbf{d} - \mathbf{f}\left( \mathbf{m},\breve{\mathbf{g}}\right)$. The objective function (the function we wish to minimize) associated with least squares is: $\sum_i \mathbf{r}_i^2$. The function optimize.leastsq is scipy's built-in least squares solver and employs the Levenberg-Marquardt algorithm in the background. The Levenberg-Marquardt algorithm is discussed in more detail in $\S$ 2.11 ➞. To use optimize.leastsq one needs a function, here called err_func, that calculates the residual vector $\mathbf{r}$. An initial guess of the parameter vector $\breve{\mathbf{g}}$ is also required.

For calibration the above variables become:

Vectorizing

• $\mathbf{d} = [\textrm{vec}(\Re\{\boldsymbol{\mathcal{D}}\}),\textrm{vec}(\Im\{\boldsymbol{\mathcal{D}}\})]$

• $\mathbf{m} = [\textrm{vec}(\Re\{\boldsymbol{\mathcal{M}}\}),\textrm{vec}(\Im\{\boldsymbol{\mathcal{M}}\})]$

• $\breve{\mathbf{g}} = [\Re\{\mathbf{g}\},\Im\{\mathbf{g}\}]$

• $\mathbf{f}\left(\mathbf{m},\breve{\mathbf{g}}\right) = [\textrm{vec}(\Re\{\boldsymbol{\mathcal{G}}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{G}}^H\}),\textrm{vec}(\Im\{\boldsymbol{\mathcal{G}}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{G}}^H\})]$, where $\boldsymbol{\mathcal{M}} = \textrm{vec}^{-1}(\mathbf{m}_U)+\imath\textrm{vec}^{-1}(\mathbf{m}_L)$ and $\boldsymbol{\mathcal{G}} = \textrm{diag}(\breve{\mathbf{g}}_U)+\imath\textrm{diag}(\breve{\mathbf{g}}_L)$

In the above bullets $\textrm{vec}(\cdot)$, $\textrm{vec}^{-1}(\cdot)$, $(\cdot)_U$, and $(\cdot)_L$ denote vectorization, matrization, the upper half of a vector and the lower half of a vector respectively. Moreover, $\Re\{\cdot\}$ and $\Im\{\cdot\}$ denote the real and imaginary part of their operands. The first thing we need to define in order to perform calibration by using optimize.leastsq is the function err_func, which we do below.

We are now able to define a wrapper function create_G_LM that in turn calls optimize.leastsq. The wrapper function translates the calibration problem into a format that optimize.leastsq can interpret. The input of create_G_LM is $\boldsymbol{\mathcal{D}}$ and $\boldsymbol{\mathcal{M}}$, while the output is $\mathbf{g}$ and $\boldsymbol{\mathscr{G}}=\mathbf{g}\mathbf{g}^H$.

We may now calibrate $\boldsymbol{\mathcal{D}}$ by using create_G_LM.

The above function works by vectorizing the real and imaginary part of $\boldsymbol{\mathcal{D}}$ and storing the result in $\mathbf{d}$. The vector $\mathbf{m}$ is generated in a similar manner. The error vector $\mathbf{r}$ is calculated by err_func. We initialize $\breve{\mathbf{g}}$ with $\breve{\mathbf{g}}_0=[\mathbf{1},\mathbf{0}]$. We can then call

optimize.leastsq(self.err_func, g_0, args=(d, m)).

We can now calculate $\mathbf{g} = \breve{\mathbf{g}}_U+\imath\breve{\mathbf{g}}_L$ and $\boldsymbol{\mathscr{G}}=\mathbf{g}\mathbf{g}^H$. This is repeated for each observational time-slot.

8.1.5 Corrected Visibilites

Before imaging, we have to correct our observed visibilities by removing the effect that the antenna gains had on the observed visibilities. This can be accomplished by using

Correcting Visibilities
\begin{equation} \boldsymbol{\mathcal{D}}^\mathrm{(c)} = \boldsymbol{\mathcal{G}}^{-1}\boldsymbol{\mathcal{D}}\boldsymbol{\mathcal{G}}^{-H} = \boldsymbol{\boldsymbol{\mathscr{G}}}^{\odot-1}\odot\boldsymbol{\mathcal{D}}, \end{equation}

• $\boldsymbol{\mathcal{D}}^\mathrm{(c)}$ is the corrected visibility matrix.
• $\boldsymbol{\mathscr{G}}^{\odot-1}$ denotes the visibility calibration matrix, which is computed by taking the Hadamard (element-wise) inverse of $\boldsymbol{\mathscr{G}}$.

The superscript $(\cdot)^{-1}$ denotes matrix inversion, while $(\cdot)^{-H}$ denotes the inverse of the Hermitian transpose. The operator $\odot$ denotes the Hadamard product. We calculate the corrected visibilities below.

We plot the corrected visibilities below. Note that the model and corrected visibilities align well, implying that calibration was successfull.

The images above contain the real part of the corrected (green) and uncorrupted (blue) visibilities as a function of timeslot for baseline 01, 02 and 12 respectively.

Next: 8.2 1GC calibration