The Radio Interferometry Inverse Problem
========================================

What you'll learn
-----------------
This page explains the core mathematical and computational challenges in radio interferometric imaging, and how Pyralysis addresses them.

.. note::
   **Key Concept:** Radio interferometry reconstructs images from incomplete, irregularly sampled Fourier data (visibilities).

Radio interferometry allows astronomers to reconstruct images of the sky using an array of antennas :cite:`thompson_moran_swenson_2017`,
which act as a single telescope with a much larger effective aperture.

However, reconstructing an image from raw interferometric data is a **complex inverse problem**
that requires advanced computational techniques.

This document provides an overview of:

* **What a visibility is** and how it is formed.
* **Why interferometric datasets are irregularly sampled in Fourier space**.
* **The mathematics of gridding and its effects on image formation**.
* **Why the inverse problem is ill-posed** and requires deconvolution.
* **Different approaches to solving the inverse problem**.

---

What is a Visibility?
---------------------
Each pair of antennas in an interferometric array forms a **baseline**,
which measures a single **visibility** at a given time and frequency.

.. note::
   **Common Pitfall:** Visibilities are not images! They are samples in Fourier space.

A **visibility** is defined as the **temporal average of the correlation of the electric fields** received by two antennas:

.. math::

   V_{ij} = \langle E(r_i,t) E(r_j,t)^* \rangle_t

where:

* :math:`V_{ij}` is the measured visibility between antennas :math:`i` and :math:`j`.
* :math:`E(r_i,t)` and :math:`E(r_j,t)` are the electric fields received at antenna positions :math:`r_i` and :math:`r_j`.
* The brackets :math:`\langle \rangle_t` indicate a temporal average.

It was later discovered that visibilities are **related to the Fourier transform of the sky brightness distribution**
through the **Van Cittert-Zernike theorem** :cite:`thompson_moran_swenson_2017`:

.. math::

   V_{ij}(u,v,w) = \int \int \frac{I(l, m)}{\sqrt{1 - l^2 - m^2}} e^{-2\pi i (ul + vm + w(\sqrt{1 - l^2 - m^2} - 1))} \, dl \, dm

where:

* :math:`I(l,m)` is the true sky brightness distribution.
* :math:`(u,v,w)` are the spatial frequency coordinates sampled by the baseline.
* The **:math:`w`-term correction** accounts for non-coplanar baselines.

Thus, **radio interferometry does not directly measure an image**, but instead collects **samples in Fourier space**.

---

Gridding: From Irregular Data to a Regular Fourier Grid
-------------------------------------------------------

Gridding is the process of interpolating measured visibilities onto a **regular Fourier grid**
so that an **Inverse Fast Fourier Transform (IFFT2)** can be applied.

.. note::
   **Key Concept:** Gridding enables fast image reconstruction using FFTs.

**1️⃣ Count-in-Cell Gridding (Basic Approach)**
The simplest way to grid visibilities is **Count-in-Cell (CIC)**, where we assign each visibility
to the nearest Fourier grid cell:

.. math::

   V_{\text{gridded}}(u_g, v_g) = \sum_{k \mid (u_k, v_k) \in C_g} W_k V_k

where:

* :math:`(u_g, v_g)` is the **center of grid cell** :math:`C_g`, where the gridded visibility is computed.
* The summation includes **all visibilities \( V_k \) where \( (u_k, v_k) \in C_g \)**.
* :math:`W_k` are the **weights** assigned to each visibility.

However, CIC introduces **aliasing artifacts** because of the **sharp binning edges**.

**2️⃣ Convolutional Gridding (Aliasing Reduction)**
To reduce aliasing, we convolve visibilities with a **gridding convolution function** :math:`G(u,v)`:

.. math::

   V_{\text{gridded}}(u_g, v_g) = \sum_k W_k V_k G(u_g - u_k, v_g - v_k)

where:

* :math:`G(u,v)` is a **convolution kernel**, typically a **prolate spheroidal function**.
* This method smooths the gridding process and **reduces Fourier artifacts**.

Once visibilities are gridded, we apply **IFFT2** to obtain the **dirty image**.

**3️⃣ Gridding Correction Function (GCF)**
Since **convolution in Fourier space corresponds to multiplication in image space**,
the resulting image after IFFT2 is **multiplied by the IFFT2 of the gridding kernel**.
This introduces a **gridding function distortion** that must be corrected.

To obtain the correct sky brightness, we **divide the image by the Gridding Correction Function (GCF)**:

.. math::

   I_{\text{corrected}}(l, m) = \frac{I_{\text{dirty}}(l, m)}{\text{GCF}(l, m)}

where:

* :math:`I_{\text{corrected}}(l, m)` is the **final corrected image**.
* :math:`I_{\text{dirty}}(l, m)` is the **dirty image**.
* :math:`\text{GCF}(l,m)` is the **IFFT2 of the gridding kernel**.

---

Deconvolution: Cleaning the Dirty Image
---------------------------------------

Since gridding modifies the visibility function, the corresponding effect in **real space**
is a convolution with the **Point Spread Function (PSF)** or **Dirty Beam**:

.. math::

   I_{\text{dirty}}(l, m) = PSF(l, m) * I_{\text{true}}(l, m)

where:

* :math:`PSF(l,m)` is the **point spread function** (or dirty beam).
* :math:`*` represents convolution.

Thus, deconvolution is required to recover :math:`I_{\text{true}}(l,m)`.

---

Approaches to Solving the Inverse Problem
-----------------------------------------

1️⃣ **Hogbom CLEAN Algorithm (Traditional Method)**
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A simple iterative method that assumes the sky is composed of **point sources** :cite:`hogbom_1974,cornwell_2008_msclean`.

**Algorithm:**
1. Find the **brightest pixel** in the dirty image.
2. Subtract a scaled **PSF shifted to that position**.
3. Repeat until the residuals are below a threshold.
4. Convolve the recovered components with a **clean beam**.

**Limitations:**

* Works well for **sparse images** but fails for **extended emission**.
* Does **not solve a proper inverse problem**, relying on heuristics.

2️⃣ **Regularized Maximum Likelihood (RML) Imaging**
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Formulates image reconstruction as an **optimization problem** with multiple regularization terms :cite:`carrillo_2014_purify`:

.. math::

   \arg \min_{I} \quad || V_{\text{obs}} - V_{\text{model}} ||^2 + \sum_i \alpha_i \mathcal{R}_i(I)

where:

* :math:`\mathcal{R}_i(I)` are different **regularization terms** (entropy, sparsity, TV, etc.).
* :math:`\alpha_i` are **tuning parameters** that control their influence.

**Challenge:**
Choosing the **optimal** :math:`\alpha_i` is difficult and usually requires **cross-validation techniques**.

3️⃣ **Compressive Sensing (CS)**
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Expresses the image as a **linear combination of basis functions** :cite:`donoho_2006_cs,candes_romberg_tao_2006,wiaux_2009_cs_ri`:

.. math::

   I(l,m) = \sum_k c_k \phi_k(l,m)

where:

* :math:`\phi_k(l,m)` are basis functions (e.g., wavelets).
* :math:`c_k` are the **coefficients**, assumed to be **sparse**.

**Optimization problem:**
- Solve for :math:`c_k` with an :math:`\ell_1`-norm constraint:

.. math::

   \arg \min_{c_k} \quad || V_{\text{obs}} - V_{\text{model}} ||^2 + \lambda || c_k ||_1

4️⃣ **Neural Networks**
^^^^^^^^^^^^^^^^^^^^^^

* Uses deep learning to learn **complex priors**.
* Can reconstruct images **directly from visibilities**.

**Challenges:**

* Requires **large training datasets**.
* Needs **interpretability and uncertainty quantification**.

---

**Conclusion**
^^^^^^^^^^^^^^^^^^^^^^

The radio interferometric inverse problem is **ill-posed**, requiring:

* **Advanced mathematical models** for accurate image reconstruction.
* **Scalable computational techniques** to handle large-scale datasets.

Further Reading
---------------
- :doc:`gridding`
- :doc:`optimization`
- :doc:`autoapi/index`

References
----------

.. bibliography::
   :filter: docname in docnames

----

:doc:`gridding` | :doc:`optimization`